Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Strike-slip fault structure and fault-system evolution : a numerical study applying damage rheology Finzi, Yaron 2010

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2010_spring_finzi_yaron.pdf [ 4.23MB ]
Metadata
JSON: 24-1.0052947.json
JSON-LD: 24-1.0052947-ld.json
RDF/XML (Pretty): 24-1.0052947-rdf.xml
RDF/JSON: 24-1.0052947-rdf.json
Turtle: 24-1.0052947-turtle.txt
N-Triples: 24-1.0052947-rdf-ntriples.txt
Original Record: 24-1.0052947-source.json
Full Text
24-1.0052947-fulltext.txt
Citation
24-1.0052947.ris

Full Text

STRIKE-SLIP FAULT STRUCTURE AND FAULT-SYSTEM EVOLUTION: A NUMERICAL STUDY APPLYING DAMAGE RHEOLOGY by  YARON FINZI  B.Sc., The Hebrew University of Jerusalem, Israel, 2001 M.Sc., The Hebrew University of Jerusalem, Israel, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Geophysics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  January 2010 © Yaron Finzi, 2010  ABSTRACT  In seismically active regions, faults nucleate, propagate, and form networks that evolve over time. Progressive strain localization and periodic fault pattern re-configuration induce the accumulation and healing of fault zone damage. The damage zones are characterized by distributed fractures, veins, and secondary faults, and may act as barriers for propagating earthquake ruptures, or as nucleation sites for earthquakes. They interact with seismic waves, promoting strong surface motions during earthquakes, and can focus fluid flow and enhance mineralization. In spite of their great scientific, social, and economic significance, interactions between these evolving damage zones and crustal deformation remain unresolved. Indeed, geodynamic models generally treat active faults as surfaces embedded in a medium with nonevolving material properties. For my dissertation projects, I have simulated fault system evolution over thousands of years, applying a rheological model which incorporates concepts of damage mechanics. This model accounts for crack nucleation, growth and concentration (i.e., material degradation), macroscopic failure, and material healing. My Simulations show that strike-slip faults form as segmented structures before evolving into contiguous, simpler structures. Flower structures rapidly form along fault segments (before a total offset of 0.05 km), and stepovers display extensive, permanent damage and ongoing seismicity throughout the seismogenic crust. My models also indicate that the “effectiveness” of material healing strongly affects the spatial extent of damage zones and long-term fault complexity. Effective healing promotes rapid evolution of segmented faults to a simpler through-going fault, and ineffective healing preserves fault complexities, resulting in long-lasting, distributed deformation. I also find that lateral contrasts in lithosphere viscosity structure (or effective plate thickness) attract evolving faults and cause damage and strain to concentrate on the “weaker” side. Realistic contrasts in crustal rigidity, however, have only a minor effect on the symmetry of damage, deformation, or the propagation of faults. In addition, lower crust and mantle viscosity ii  contrasts induce the formation of broad shear zones with relatively high strain-rate in the “weaker” side of the interface. These results demonstrate that reasonable, lateral contrasts in viscosity (rather than extreme, unrealistic contrasts in elasticity) can explain GPS observations of highly asymmetric, interseismic deformation around strike-slip faults.  iii  TABLE OF CONTENTS Abstract ..................................................................................................................................... ii Table of contents ...................................................................................................................... iv List of tables ........................................................................................................................... viii List of figures ........................................................................................................................... ix Acknowledgments.................................................................................................................... xi Dedication ............................................................................................................................... xii Co-authorship statement ........................................................................................................ xiii  1 Introduction ......................................................................................................................... 1 1.1 Background ................................................................................................................... 1 1.2 Thesis themes ................................................................................................................ 3 1.3 Format and outline of thesis .......................................................................................... 7 Bibliography ....................................................................................................................... 9  2 Structural properties and deformation patterns of evolving strike-slip faults: Numerical simulations incorporating damage rheology ..................................................................... 14 Summary ........................................................................................................................... 14 2.1 Introduction ................................................................................................................. 15 2.2 Geological and geophysical observations of fault zone structure ............................... 16 2.3 Damage rheology framework ...................................................................................... 19 2.3.1 Theoretical background .................................................................................. 19 2.3.2 Damage model parameters .............................................................................. 22 2.4 Geophysical constraints on healing rate parameters ................................................... 24 iv  2.4.1 Parameters C1 and C2, and healing as a function of time................................ 24 2.4.2 In-Situ geophysical constraints for healing parameters C1 and C2.................. 28 2.5 Damage and strain distribution across active strike-slip faults ................................... 31 2.5.1 Model setup ..................................................................................................... 31 2.5.2 Model output: Examples and interpretation .................................................... 33 2.5.3 Damage, rigidity and strain distribution across strike-slip fault segments ..... 40 2.5.4 Fault stepovers ................................................................................................ 45 2.5.5 Fault system complexity as a function of time ............................................... 46 2.6 Discussion and conclusions ........................................................................................ 50 Acknowledgments............................................................................................................. 52 Bibliography ..................................................................................................................... 53  3 Damage-zone healing and strike slip fault evolution: Numerical results and geophysical constraints ......................................................................................................................... 61 Summary ........................................................................................................................... 61 3.1 Introduction ................................................................................................................. 62 3.2 Damage rheology and healing formulation................................................................. 63 3.2.1 Damage rheology: Theory .............................................................................. 64 3.2.2 Healing of damaged material: Constraints from fault zone observations ....... 68 3.3 Results ......................................................................................................................... 70 3.3.1 Fault zone structure as a function of healing .................................................. 70 3.3.2 Fault zone evolution patterns and temporal stability of fault stepovers ......... 75 3.4 Discussion and conclusions ........................................................................................ 80 Bibliography ..................................................................................................................... 84 v  4 Effect of lateral rheological contrasts on damage zone structure and surface deformation around strike-slip faults..................................................................................................... 88 Summary ........................................................................................................................... 88 4.1 Introduction ................................................................................................................. 89 4.2 Methods and model set-up .......................................................................................... 90 4.2.1 Numerical model and damage rheology ......................................................... 90 4.2.2 Model set-up ................................................................................................... 92 4.3 Results ......................................................................................................................... 97 4.3.1 Damage zone asymmetry and fault location relative to the material interface ................................................................................................................... 97 4.3.2 Fault propagation toward and along bi-material interfaces .......................... 102 4.3.3 Asymmetric deformation patterns and surface velocities ............................. 106 4.4 Discussion and conclusions ...................................................................................... 113 4.4.1 Crustal interfaces, mantle interfaces and crust-mantle interactions.............. 113 4.4.2 Interpretation of surface velocity data along material interfaces .................. 116 4.4.3 Propagation of faults to material interfaces .................................................. 120 4.4.4 Conclusions ................................................................................................... 121 Bibliography ................................................................................................................... 123  5 Concluding chapter ......................................................................................................... 128 5.1 Main contributions .................................................................................................... 128 5.1.1 Establishing damage rheology applicability to studies of natural fault systems ................................................................................................................. 129  vi  5.1.2 Damage-zone structure and deformation patterns along segmented strike-slip faults ........................................................................................................ 130 5.1.3 Revealing the important role healing efficiency plays in fault system evolution ................................................................................................................. 131 5.1.4 Asymmetric surface deformation due to lateral contrasts in rheology, and interpretation of geodetic data ................................................................ 132 5.2 Updates on working hypotheses and on current knowledge in the field .................. 132 5.3 Weakness of reported research and outline of potential future research .................. 135 5.3.1 Loose constraints on material parameters and sparse observation of in situ damage evolution .................................................................................... 135 5.3.2 Computational limitations ............................................................................. 136 Bibliography ................................................................................................................... 139  Appendix A  Numerical procedures ................................................................................... 141  Bibliography ................................................................................................................... 143 Appendix B  Boundary conditions ................................................................................... 144  Bibliography ................................................................................................................... 146 Appendix C  Code verification: Viscoelastic component .................................................. 147  Bibliography ................................................................................................................... 149  vii  LIST OF TABLES 2.1 Damage rheology parameters and their constraints……………………………….. 23 2.2 Model parameters used in our fault evolution study……………………...............  32  4.1 Material parameters used in simulations with lateral material contrasts ………....  95  4.2 Asymmetry of fault damage zones along material interfaces ………………….....  99  4.3 Fault-parallel surface deformation patterns along material interfaces …………… 111  viii  LIST OF FIGURES 2.1. Healing parameter space and expected damage level as a function of time……… 27 2.2. Geophysical, analytical and experimental constraints on healing parameters……. 30 2.3. Typical lithospheric structure used in the numerical simulations……………....... 33 2.4. Plan views of a segmented strike slip fault at several depths…………………….. 35 2.5. Cross-sections of a typical “flower type” damage zone………………………….. 36 2.6. DOF damage zone width and depth plotted against cumulative offset………….. 38 2.7. Sensitivity of model results to element size……………………………………… 39 2.8. DOF damage zone width and depth as a function of material parameters………. 41 2.9. LAF damage zone width and depth as a function of healing parameter C2……… 42 2.10. Correlation between damage and strain across simulated fault-zones………….. 44 2.11. Cross-sections of damage patterns within fault stepovers………………………. 46 2.12. Fault stepover evolutionary stages………………………………………………. 48 2.13. Fault evolution snap-shots………………………………………………………. 49 3.1. Damage as a function of time for various sets of parameters C1 and C2…………. 67 3.2. Geophysical, analytical and experimental constraints on healing parameters…… 69 3.3. DOF and LAF damage along a simulated strike-slip fault……………………..... 72 3.4. Damage zones along faults with effective and ineffective healing conditions......  74  3.5. Sensitivity of damage zone dimensions to healing effectiveness ……………….. 75 3.6. Slip-rate and damage zone evolution with time………………………………...... 78 3.7. Stepover evolution under effective and ineffective healing conditions.................  80  4.1. Effective viscosity profiles of simulated lithologies …………………………...... 96 4.2. Damage zone asymmetry measures ……………………………..........................  98  4.3. Damage zone structures along various types of material interfaces …………...... 101 ix  4.4. Damage maps of simulated faults along various types of material interfaces …… 104 4.5. Strain rate and velocity cross-sections of faults along material interfaces ………. 108 4.6. Velocity profile asymmetry measures ……………………………………………. 111 4.7. Velocity and effective viscosity across a fault zone that evolved along a mantle viscosity interface ……………………………………………………………...... 115 4.8. Modeled and observed (GPS) velocity profiles for the NAFZ ………………...... 118 B.1. Schematic illustration of the variable force boundary condition ……………….. 145 C.1. Surface velocities from simulations of viscoelastic post-seismic deformation …. 148  x  ACKNOWLEDGMENTS First, I would like to express my gratitude to my thesis advisor, Elizabeth Hearn. Elizabeth has been a wonderful supervisor providing overwhelming support and great scientific insight. Liz guided and inspired me in my scientific and academic endeavor. Her reassuring and kind personality helped me through the tough times and provided me a great example to follow in my future career. I‟ve had the pleasure to work closely with Vladimir Lyakhovsky, whose creative approach, patience and good will helped me from my very first steps in numerical modeling, through various frustrating problems and dull moments, and on to great satisfaction from my scientific discoveries. I‟m very fortunate to have had a chance to collaborate with Yehuda BenZion. Invigorating and overwhelming discussions with Yehuda provided me important insights into geology, science and life in academia. I would also like to thank the friends, peers, staff and professors of the Earth and Ocean Department. Specifically I would like to thank my friends and office mates Ali and Jounada for many discussions and for a fresh perspective on life and mankind. I would like to thank Andrew, David, and Rebecca-Ellen for their friendship, and Cornell for his ongoing support. And I would like to thank Michael Bostock, Jim Mortensen, Garry Clarke and Erik Eberhardt for their advice, patience and support. Finally, and most importantly, I wish to thank my family. I am forever grateful for their love and support. I hope my graduation will bring my parents, siblings and grandfather pride and pleasure, and will help us all overcome the sacrifice of living so far away from each other. For sharing the love of nature, for encouraging my scientific curiosity, and for motivating me to contribute to a better future on planet Earth, I thank my two amazing sons, Ma‟ayan Raz and Dolev Bar. For being the love of my life and for making our life so amazing, I thank my wife Michelle. You are my core, you are my sun, you are my magnitude 10. Much appreciated financial support was provided by the Southern California Earthquake Center (SCEC), the Faculty of Graduate Students of UBC, and the Lorentzen Scholarship.  xi  DEDICATION To my grandfather Alberto Finzi, and parents Yair and Arna Finzi who got me started on my pursuit of happiness.  xii  CO-AUTHORSHIP STATEMENT Chapters 2 and 3 were published with my supervisor, Elizabeth Hearn, and with collaborators Vladimir Lyakhovsky and Yehuda Ben-Zion (from the Geological Survey of Israel, and University of Southern California, respectively). I identified the research program and performed the research. Vladimir Lyakhovsky assisted with algorithm development and theoretical analysis of damage rheology parameters. Elizabeth Hearn assisted in developing research ideas and in interpretation of geodetic data. Yehuda Ben-Zion assisted in many discussions about earthquakes, faults and geophysical data analysis. I wrote the manuscripts with inputs from Elizabeth Hearn. Vladimir Lyakhovsky and Yehuda Ben-Zion reviewed the papers at various stages of the work. Chapter 4 was prepared for publication with Elizabeth Hearn. The research goals were defined together, and the numerical analysis was performed by me. Elizabeth Hearn assisted in interpretation of surface velocity results. I wrote the manuscripts, and Elizabeth Hearn reviewed and commented on the manuscript.  xiii  1  INTRODUCTION  1.1  Background Many aspects of fault evolution and plate-boundary deformation remain an enigma in  spite of their great scientific, social, and economic significance. In actively deforming regions, we know that faults nucleate, propagate, and form systems that mature over time, but the interactions between evolving structural properties of the crust (including the fault zones) and patterns of stress and deformation remain unresolved. A decade-long surge of activity in the geological and geophysical communities has been aiming to resolve fundamental questions related to the spatio-temporal evolution of fault systems and their associated seismicity (see reviews by Ben-Zion and Sammis, 2003; Rundle et al., 2003). Various studies of active and exhumed strike-slip faults have yielded high resolution descriptions of their internal structure and slip localization patterns. These studies indicate that most of the slip along a fault is accommodated along a single „„principal fracture surface‟‟ within a narrow (10–20 cm) layer of ultracataclasite. The ultracataclasite layer referred to as the „„core‟‟ of the fault zone, is parallel to the macroscopic slip vector and is surrounded by a damage zone which is a relic structure of the progressive coalescence and localization of the active fault zone (e.g., Sibson, 1986; Chester et al., 1993). The damage zone consists of various geologically observable features such as open and healed fractures, veins, and secondary faults that are associated with deformation along the main fault zone (Chester and Logan, 1986; Chester and Chester, 1998; Schulz and Evans, 2000; Rockwell and Ben-Zion, 2007). While damage along fault segments is typically confined to the top 3-5 km of the crust (Peng and Ben-Zion, 2004; Lewis et al., 2005; Dor et al., 2006, 2008), persisting fault complexities that exhibit strain hardening at greater depths (e.g., fault offsets, kinks, and bends) sustain larger and deeper damage zones consisting of secondary fractures at different scales (e.g., Segall and Pollard, 1980; Ben-Zion and Sammis, 2003; Kim et al., 2004; De Paola et al., 2007). It is currently understood that the internal structure of fault systems and their damage evolution patterns have important implications for fault propagation (e.g., Vermilye and Scholz, 1  1998, 1999) and earthquake initiation and arrest (Sibson, 1985; King, 1986). Persisting rigidity loss and local modification of the stress state within fault stepovers and other fault complexities may cause a local change in attenuation and ground shaking patterns, and could yield a barrier for propagating earthquakes or a nucleation zone for earthquakes (Harris, 2004; Oglesby, 2005; Wesnousky, 2006). Increased permeability and enhanced fluid flow through faults and cracks within stepover zones has an additional effect on rupture dynamics and implications for mineral exploration (Sibson, 1996, Martel and Boger, 1998; Micklethwaite and Cox, 2004, 2006, Sheldon and Micklethwaite, 2007). While fluid flow and pore-pressure have direct effects on the stress state and overall strength of fault segments and stepovers (e.g., Sibson, 1985), we are currently unable to simulate the coupled evolution of fluid-flow and rock deformation. Regardless of the overwhelming evidence that natural fault zones are evolving complex 3-D structures, most studies of fault systems, rupture dynamics and seismic hazards represent faults as idealized planar features. Until recently, numerical models of coupled fault-structure and deformation-pattern evolution have been confined to simulations of seismic activity along a single fault or a simple fault system (e.g., Robinson and Benites, 1995; Rice and Ben-Zion, 1996; Ziv and Rubin, 2003; Lapusta and Rice, 2003). Some models applied 2-D calculations of elastostatic deformation in a thin plate to simulate the development of regional faults (e.g., Sornette et al., 1994; Ward, 1996). However, such models typically neglected the time-dependent inelastic deformation of the ductile lower crust and its important role in loading faults in the brittle upper crust (e.g., Thatcher, 1983; Li and Rice, 1987; Reches et al., 1994), and they often neglected the fact that the geometry and rheological properties of fault systems evolve with ongoing deformation (e.g., King, 1986; Andrews, 1989; Oglesby et al., 2003). A few recent models that impose fault surfaces governed by friction employ either a Coulomb stress analysis or a plasticity analysis to calculate off fault strain or to explain distributed deformation (e.g., Ben-Zion and Shi, 2005; Templeton and Rice, 2008; Ma, 2008). While these models can account for evolving rheological properties of fault zones and may be used to study off-fault deformation, they do not enable fault-system evolution in the form of generation of new faults. To simulate brittle rock deformation including fault nucleation and propagation (i.e. creation of new frictional surfaces) and the generation of distributed deformation, a rheological model must incorporate concepts of fracture mechanics. Therefore, realistic models of the fault 2  system evolution should employ concepts such as damage rheology that account for all stages of the seismic cycle including subcritical crack growth during early stages of the loading, material degradation due to increasing crack concentration, macroscopic failure, off-fault deformation, and healing (e.g., Rundle et al., 2003; Turcotte and Glasscoe, 2004; Bercovici and Ricard, 2003; Regenauer-Lieb and Yuen, 2003).  1.2  Thesis themes In my doctoral research I have applied a recently developed damage rheology model  (Lyakhovsky et al. 1997a,b, 2001) in numerical simulations of strike-slip fault systems to study their evolving structure and deformation patterns. The damage rheology model used in my work is based on the interdisciplinary view that the brittle portion of the Earth's lithosphere is “damaged” in the sense that it contains an internal distribution of microcracks, joints and faults (e.g., Brace and Kohlstedt, 1980). During continuing deformation, the material around such flaws is subjected to large stress concentrations which necessarily lead to increased cracking (damage accumulation) and a reduction of the elastic strength (Cowie and Scholz, 1992). Subsequently, when shear stresses drop (e.g., after a seismic failure of a fault), the internal damage could undergo healing, which increases the elastic strength of the material (e.g., Rice, 1983). In damage rheology models, an evolving damage variable quantifies the extent of internal damage and modifies the elastic properties of material in the crust, simulating the creation and healing of fault systems (Malvern, 1969; Krajnovic, 1996). Various damage rheology frameworks are currently used in a wide range of deformation studies, and are considered to be most appropriate for analyzing large scale geologic deformation in the upper crust (e.g., Turcotte and Glasscoe, 2004; Bercovici and Ricard, 2003; Regenauer-Lieb and Yuen, 2003; Allix and Hild, 2002; Auth et al., 2003).    Preparatory theme: Improving our damage model and constraining its parameters. The damage rheology model used in my work was previously applied by Ben-Zion et al.  (1999), Lyakhovsky et al. (2001) and Ben-Zion and Lyakhovsky (2006) to study deformation patterns associated with large earthquake cycles, frequency-size statistics of earthquakes, 3  accelerated seismic release and more. In recent years, Lyakhovsky et al. (2005) and Lyakhovsky and Ben-Zion (2008) demonstrated that simulations with damage rheology can reproduce widely observed phenomena of rate- and state-dependent friction, and they improved the coseismic strain calculation in their framework so it is analogous to Drucker-Prager plasticity. However, the model parameters were previously constrained based on analytic considerations and laboratory results (Hamiel et al., 2004, 2006; Lyakhovsky et al., 1997a,b) rather than fault zone observations. Furthermore, though the code incorporates viscoelasticity, calculations without damage evolution have been benchmarked against just a few, simple analytical solutions. To improve our damage rheology model and to further establish its utility to the Earth science community, I benchmark our calculation of mantle viscoelastic relaxation and I demonstrate the applicability and advantages of a non-standard variable force boundary condition (Chapter 2). One very significant contribution to the damage rheology model is an analysis of the healing formulation and a refinement of the constraints on healing parameters based on seismic and geodetic data from various strike-slip faults (Chapters 2, 3). It is pertinent to realize that the damage rheology applied in this work does not explicitly describe any micromechanism of healing, and rather it adopts a well-accepted healing rate form based on friction experiments (Dieterich 1978, 1979). While the healing formulation is the most suitable for crustal processes, the calibration of healing parameters is challenging due to the scarcity of insitu measurements of damage levels and healing rates during the various stages of the earthquake cycle. The great importance of our constraints on these parameters is twofold: As these parameters have a controlling effect on damage evolution, constraining them based on fault-zone observations establishes the applicability of the model to studies of natural fault evolution (Chapter 2). In addition, the improved analysis of healing parameters reveals the important role that healing effectiveness plays in the evolution of fault structures and fault systems (Chapter 3).    First theme: Fault-zone properties and structures, and strain distribution patterns The first objective of my thesis is to numerically study the evolution of strike-slip fault  structures and associated deformation patterns, and to determine their dependence on lithosphere characteristics. In chapter 2, I present simulations of evolving strike-slip faults that explore the characteristics of commonly observed fault zone structures (flower structures, distributed 4  damage, fault segments and stepovers), strain localization along active fault segments, and interseismic seismicity and damage accumulation within fault stepovers. The flower structures in our simulations consist of a broad damage zones in the top few kilometers of the crust and highly localized damage at depth. The flower structures form during an early evolutionary stage of the fault system, and they persist as continued deformation localizes further along narrow slip zones. While tectonic strain is concentrated along the highly damaged fault cores at depth, near the surface a small portion of the strain is accommodated over a broader region that correlates with the distributed shallow damage. Chapter 2 also shows that simulated fault stepover zones are locations of ongoing interseismic deformation and significant damage (reduced rigidity and shear wave velocity) to depths of 10 to 15 km. Results presented in chapter 2 indicate that the damage zone structures are robust features that show little or no dependence on crustal properties and model parameters. However, an improved healing analysis based on additional fault zone observations (Chapter 3) reveals that the damage zone dimensions show significant variations as a function of the model parameters that control healing effectiveness.    Second theme: Fault-system evolution and temporal stability of large-scale fault complexities The second objective of my thesis is to study the structural evolution of fault structures  from the single-fault scale to the fault-system (or plate boundary) scale. In the seismogenic zone, which is governed by a strain weakening rheology, faults are expected to evolve toward geometrically simple structures (Ben-Zion and Sammis, 2003). This is because strain weakening produces zones of localized deformation and strength reduction, leading to further strain localization and strength reduction. Along natural faults, however, heterogeneities of material properties, stress state and strain history, and dynamic processes produce complications that prevent complete localization. Observational data (e.g., Evans et al., 2000; Chester et al., 1993; de Joussineau and Aydin., 2007) indicate that during a short initial phase, fault system complexity increases due to strain hardening and formation of fault structures of various length scales. With increasing deformation under the same applied loads, this process is replaced by strain weakening and localization to tabular zones that become the main carriers of subsequent strain. At that stage, most of the complex initial structure becomes mechanically passive and the 5  fault zones evolve with continuing deformation toward planar (curvilinear) geometry, progressive simplicity and regularization. In addition, Wesnousky (1994) and Stirling et al. (1996) observed a transition from the Gutenberg-Richter power-law distribution of earthquake magnitudes on relatively disordered immature faults, to the characteristic earthquake distribution on relatively regular mature faults. These observed seismicity and structural evolution patterns are compatible with various computer simulations and analytical results from fault system models (Fisher et al., 1997; Dahmen et al., 1998; Ben-Zion et al., 1999; Lyakhovsky and BenZion, 2009). In chapter 3, I present simulations that explore various aspects of fault system evolution including the structural evolution of individual segments, the evolution of fault slip-rate, and the evolution of segmented fault systems towards a simpler, through-going fault. In addition I show how model parameters that control healing effectiveness play an important role in the evolution of fault structures and fault systems. Results presented in chapter 3 indicate that highly effective healing leads to a wide and shallow damage zone along fault segments, and short-lived, shallow damage zones at fault stepovers which do not impede the evolution of the segmented fault into a through-going fault. Inefficient healing in our simulations leads to a narrow and deep damage zone along fault segments, and long-lived, deep and extensively-damaged fault stepovers which do not heal and therefore impede the evolution of the segmented fault into a through-going fault. Evolution patterns of slip rate along simulated strike-slip faults (presented in chapter 3) complement and support my results on damage zone formation (Chapter 2). In agreement with field observations cited above, my numerical results (Chapters 2 and 3) suggest that the structural maturation of a fault segment is rapid and occurs before large offsets accumulate.    Third theme: strike-slip fault deformation patterns along bi-material interfaces The third objective of my thesis is to use numerical simulations to study the effects of  large-scale lithospheric heterogeneities on the evolution of fault systems and their associated deformation patterns. Large faults juxtapose materials that generally have different physical properties. Thus, in the simplest strike-slip case, asymmetry of deformation across the fault reflects asymmetry in material properties. In some cases large-scale continental strike-slip faults juxtapose oceanic lithosphere to continental lithosphere. Such extreme cases exist for the San 6  Andreas, Alpine and Philippine faults (Le Pichon et al., 2005). Large damage zones generated by repeated earthquakes, and narrow slip zones that tend to localize on one side of these damage zones, constitute another form of large-scale bi-material interface (e.g., Ben-Zion and Andrews, 1998; Ben-Zion and Sammis, 2003). Various studies suggest that tectonic settings such as oceanic-continental interfaces and embedded weak zones may control the structure and evolution of strike-slip fault systems, which in turn control the seismicity patterns and spatial distribution of deformation. Brietzke and BenZion, (2006) show that 2-D in-plane ruptures are attracted to and tend to migrate toward material interfaces such as embedded (damaged) weak zones. Le Pichon et al. (2003), Le Pichon and Kreemer (2005), and Wdowinski et al. (2007) use geodetic measurements to analyze significant asymmetry in interseismic and coseismic strain along large continental strike-slip faults. In their work they indicate that intensive damage along fault zones may lead to asymmetric deformation patterns. Schmalzle et al. (2006) use finite element models to investigate the effect that largescale material heterogeneities have on surface deformation patterns. They conclude that imposing homogeneous mechanical properties in a model may bias the conclusions about fault behavior and cause erroneous interpretations of geodetic data. In chapter 4 I test three hypotheses regarding how bi-material interfaces affect the structure and evolution patterns of strike-slip faults. Simulations with various configurations of bi-material interfaces relative to model boundaries test whether evolving faults tend to migrate toward and propagate along such interfaces. Other simulations focus on the damage zone structure and surface velocity profiles across strike-slip faults that coincide with bi-material interfaces. These results address the growing consensus that asymmetric surface deformation patterns reflect a structural asymmetry of fault zones and large scale lithospheric heterogeneities, and suggest some considerations in the interpretation of geodetic observations of such deformation patterns.  1.3 Format and outline of thesis This thesis is an amalgamation of published and submitted manuscripts resulting from my PhD research at the University of British Columbia and collaboration with co-authors from 7  different institutions. Following the regulations of UBC Faculty of Graduate Studies, minor format changes were made. Appendices from all of the manuscripts were combined and placed after the concluding chapter of the thesis. As the content of manuscripts was unaltered to conform to journal copyrights, some repetitions of introductory background and theory are expected. The thesis introduction (Chapter 1) outlines the field of study and reviews the literature and concepts required to develop a basic understanding of the research. The Introduction also describes the themes of my research and the relations between the manuscript chapters (Chapters 2-4). The concluding chapter of the thesis discusses the results and conclusions of the manuscript chapters and relates them to each other and to the overall field of study. The concluding chapter also comments on the significance of contributions and the limitations of the thesis research, and it proposes several ideas for future work to improve and complement the presented research.  8  Bibliography Allix, O. Hild, F. (2002), Continuum Damage Mechanics of Materials and Structures, Elsevier, p. 396. Amsterdam Holland. Andrews, D.J., 1989. Mechanics of fault junctions, J. geophys. Res., 94, 9389–9397. Auth, C., D. Bercovici, and Christensen, U. (2003), Two-dimensional convection with a self lubricating, simple-damage rheology, Geophys. J. Int., 154, 783–800. Ben-Zion, Y. and Andrews, D. J. (1998), Properties and implications of dynamic rupture along a material interface, Bull. Seism. Soc. Am., 88(4), 1085–1094. Ben-Zion, Y., Dahmen, K., Lyakhovsky, V., Ertas, D. Agnon, A. (1999). Self-driven mode switching of earthquake activity on a fault system, Earth planet. Sci. Lett., 172, 11–21. Ben-Zion, Y. Lyakhovsky, V. (2002), Accelerating seismic release and related aspects of seismicity patterns on earthquake faults, Pure and appl. Geophys., 159, 2385–2412. Ben-Zion, Y. and Sammis, C. (2003), Characterization of Fault Zones, Pure and Appl. Geophys. 160, 677–715. Ben-Zion, Y. and Shi, Z. (2005), Dynamic rupture on a material interface with spontaneous generation of plastic strain in the bulk, Earth Planet. Sci. Lett., 236, 486-496, DOI: 10.1016/j.epsl.2005.03.025. Ben-Zion, Y. and Lyakhovsky, V. (2006), Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology, Geophys. J. Int. 165, 197-210. Bercovici, D. and Ricard, Y. (2003), Energetics of a two-phase model of lithospheric damage, shear localization and plate-boundary formation, Geophys. J. Int. 152 (3), 581–596 doi:10.1046/j.1365-246X.2003.01854.x Brace, W. F. and Kohlstedt, D. L. (1980), Limits on the lithospheric stresses imposed by laboratory experiments. J. Geophys. Res. 85, 6248-6252. Brietzke, G. B. and Y. Ben-Zion, (2006), Examining tendencies of in-plane rupture to migrate to material interfaces, Geophys. J. Int., 167, 807–819. Chester, F.M., Logan, J.M. (1986), Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California. Pure Appl. Geophys. 124, 79–106. Chester, F. M., J. P. Evans, and Biegel, R. L. (1993), Internal structure and weakening mechanisms of the San Andreas Fault, J. Geophys. Res. 98, 771–786. Chester, F. M. and Chester, J. S. (1998), Ultracataclasite structure and friction processes of the Punchbowl Fault, San Andreas System, California, Tectonophysics 295 199–221. Cowie, P.A. and Scholz, C.H. (1992), Physical explanation for the displacement–length relationships of faults using a post-yield fracture mechanics model. J. Struct. Geol. 14, 1133–1148. Dahmen, K., D. Ertas, and Y. Ben-Zion, (1998), Gutenberg Richter and Characteristic Earthquake Behavior in Simple Mean-Field Models of Heterogeneous Faults, Phys. Rev. E 58, 1494–1501. 9  de Joussineau, G., and A. Aydin (2007), The evolution of the damage zone with fault growth in sandstone and its multiscale characteristics, J. Geophys. Res., 112, B12401, doi:10.1029/2006JB004711. De Paola, N., Holdsworth, R.E., Collettini, C., McCaffrey, K. W. and Barchi, M. R. (2007), The structural evolution of dilational stepovers in regional transtensional zones, Cunningham W.D. and Mann, P. (eds), Tectonics of Strike-Slip Restraining and Releasing Bends. Geological Society, London. Special Publication, 290, pp. 433-445. Dieterich, J. H. (1978), Time-dependent friction and the mechanics of stick-slip, Pure and Appl. Geophys. 116, 790-805. Dieterich, J. H. (1979), Modeling of rock friction 1. Experimental results and constitutive equations, J. Geophys. Res. 84, 2161-2168. Dor, O., Rockwell, T. K., and Ben -Zion, Y. (2006), Geological observations of damage asymmetry in the structure of the San Jacinto, San Andreas and Punchbowl Faults in Southern California: a possible indicator for preferred rupture propagation direction, Pure appl. Geophys. 163, 301–349, DOI 10.1007/s00024-005-0023-9. Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., Duman, T. Y. (2008), Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction, Geophys. J. Int., doi: 10.1111/j.1365-246X.2008.03709.x. Evans, J. P., Shipton, Z. K., Pachell, M. A., Lim, S. J., and Robeson, K. (2000), The Structure and Composition of Exhumed Faults, and their Implication for Seismic Processes, In Proc. of the 3rd Conf. on Tect. problems of the San Andreas system, Stanford University. Fisher, D. S., K. Dahmen, S. Ramanathan, and Y. Ben-Zion, 1997. Statistics of earthquakes in simple models of heterogeneous faults, Phys. Rev. Lett., 78, 4885–4888. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y. and Lockner, D. (2004), A Visco-elastic damage model with applications to stable and unstable fracturing, Geophys. J. Int. 159, 1155–1165. Hamiel, Y., Katz, O., Lyakhovsky., Reches, Z. and Fialko, Y. (2006), Stable and unstable damage evolution in rocks with implications to fracturing of granite, Geophys. J. Int. 167, 1005-1016. Harris, R. A., 2004. Numerical Simulations of Large Earthquakes: Dynamic Rupture Propagation on Heterogeneous Faults. Pure appl. Geophys. 161, 2171–2181. Kim, Y. S., Peacock, D. C. P., and Sanderson, D. J. (2004), Fault damage zones, J. of Struct. Geology, v. 26, p. 503-517. King, G. (1986), Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure, Pure Appl. Geophys. 124, 567–585. Krajcinovic, D. (1996), Damage Mechanics, Elsevier, Amsterdam. Lapusta, N. and J. R. Rice, (2003), Nucleation and early seismic propagation of small and large events in a crustal earthquake model. J. Geophys. Res. 108 (B4): Art. No. 2205. 10  Le Pichon, X., C. Kreemer, and N. Chamot-Rooke, 2005. Asymmetry in elastic properties and the evolution of large continental strike-slip faults, J. Geophys. Res., 110, B03405, doi:10.1029/2004JB003343. Le Pichon, X., N. Chamot-Rooke, C. Rangin, and A. M. C. Sengor, 2003. The North Anatolian fault in the Sea of Marmara, J. Geophys. Res., 108(B4), 2179, doi:10.1029/2002JB001862. Lewis, M. A., Peng, Z., Ben-Zion, Y., and Vernon, F. (2005), Shallow seismic trapping structure in the San Jacinto fault zone, Geophys. J. Int. 162, 867–881, doi:10.1111/j.1365246X.2005.02684.x, 2005. Li, Y. G., P. Leary, K. Aki, and P. Malin (1990), Seismic trapped modes in the Oroville and San Andreas fault zones, Science, 249, 763-766. Li, V. C., and J. R. Rice, 1987. Crustal deformation in great California earthquake cycles, J. Geophys. Res., 92: 11,533-11,551. Lyakhovsky, V., Ben-Zion, Y. and Agnon, A. (1997a) Distributed damage, faulting, and friction, J. Geophys. Res. 102, 27,635-27,649. Lyakhovsky, V., Reches, Z., Weinberger, R. and Scott, T.E. (1997b), Non-linear elastic behavior of damaged rocks, Geophys. J. Int. 130, 157-166. Lyakhovsky, V. Ben-Zion, Y. and Agnon, A. (2001), Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere, J. Geophys. Res. 106: 4103-4120. Lyakhovsky, V., Ben-Zion, Y. and Agnon, A. (2005), A viscoelastic damage rheology and rateand state-dependent friction, Geophys. J. Int. 161, 179-190. Lyakhovsky, V., and Ben-Zion, Y. (2008), Scaling relations of earthquakes and aseismic deformation in a damage rheology model, Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365-246X.2007.03652.x. Lyakhovsky V., and Ben-Zion, Y. (2009), Evolving fault zone structures in a damage rheology model, in review, Geochem., Geophys., Geosys. Ma, S., 2008. A physical model for widespread near-surface and fault zone damage induced by earthquakes. Geochem. Geophys. Geosyst., 9, Q11009, doi:10.1029/2008GC002231. Malvern, L.E. (1969), Introduction to the Mechanics of a Continuum Medium, New Jersey, Prentice-Hall, Inc., 713 p. Martel, Y., S. J. Boger, 1998. Geometry and mechanics of secondary fracturing around small three-dimensional faults in granitic rock. Journal of Geophysical Research 103, 21299– 21314. Micklethwaite, S. and Cox, S. F. (2004), Fault-segment rupture, aftershock-zone fluid flow and mineralization, Geology, 32, 813–816. Micklethwaite, S., Cox, S.F., 2006. Progressive fault triggering and fluid flow in aftershock domains: Examples from mineralized Archaean fault systems. Earth and Planet. Sci. Lett., Vol. 250, Issues 1-2, 318-330.  11  Oglesby, D. D., 2005. The dynamics of strike-slip step-overs with linking dip-slip faults. Bull. Seismol. Soc. Am. 95, 1604–1622. Peng, Z. and Ben-Zion, Y. (2004), Systematic analysis of crustal anisotropy along the KaradereDuzce branch of the north Anatolian fault, Geophys. J. Int. 159, 253-274, doi:10.1111/j.1365- 46X.2004.02379.x. Reches, Z., G. Schubert, and C. Anderson, (1994), Modeling of periodic great earthquakes on the San Andreas fault: Effects of nonlinear crustal rheology, J. Geophys. Res., 99: 21,98322,000. Regenauer-Lieb, K., Yuen, D.A. (2003), 'Modeling shear zones in geological and planetary sciences: solid- and fluid-thermal-mechanical approaches', Earth - Science Reviews, 63, pp. 295-349. Rice, J.R., (1983), Constitutive relations for fault slip and earthquake instabilities. Pure Appl. Geophys. 121, 443–475. Rice, J.R. and Ben-Zion, Y. (1996), Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci. USA 93: 3811-1818. Rice, J. R., C. G. Sammis, and R. Parsons, (2005), Off-fault secondary failure induced by a dynamic Slip Pulse, Bull. Seis. Soc. Am., 95: 109 - 134. Robinson, R. and R. Benites, (1995), Synthetic seismicity models of multiple interacting faults, J. Geophys. Res., 100: 18,229-18,238. Rockwell, T.K., Ben-Zion Y., 2007. High localization of primary slip zones in large earthquakes from paleoseismic trenches: Observations and implications for earthquake physics. J. Geophys. Res. 112, B10304, doi:10.1029/2006JB004764. Rundle, J. B., D. L. Turcotte, R. Shcherbakov, W. Klein, and C. Sammis, (2003), Statistical physics approach to understanding the multiscale dynamics of earthquake fault systems. Rev. Geophys. 41 (4): Art. No. 1019, Dec. 18. Schmalzle, G., T. Dixon, R. Malservisi, and R. Govers, (2006), Strain accumulation across the Carrizo segment of the San Andreas Fault, California: Impact of laterally varying crustal properties, J. Geophys. Res., 111, B05403, doi:10.1029/2005JB003843. Schulz, S. E. and Evans, J. P. (2000), Mesoscopic Structure of the Punchbowl Fault, Southern California and the Geologic and Geophysical Structure of Active Strike-slip Faults, J. of Struct. Geol. 22, 913–930. Segall, P., and Pollard, D. D. (1980), Mechanics of discontinuous faults, J. Geophys. Res. 85, 4337-4350, 1980. Sheldon H. A., and Micklethwaite S. (2007), Damage and permeability around faults: Implications for mineralization, Geology: Vol. 35, No. 10 pp. 903–906. Sibson, R.H. (1985), Stopping of earthquake ruptures at dilational fault jogs, Nature 316, 248– 251. Sibson, R.H. (1986), Brecciation processes in fault zone: Inferences from earthquake rupturing. Pure Appl. Geophys. 124, 159–175. 12  Sibson, R.H., (1996), Structural permeability of fluid-driven fault-fracture meshes. Journal of Structural Geology 18, 1031–1042. Sibson, R.H. (2003), Thickness of the Seismic Slip Zone, Bull. Seis. Soc. Am. 93, 3, 1169-1178. Sornette, D., P. Miltenberger, and C. Vanneste, 1994. Statistical physics of fault patterns selforganized by repeated earthquakes, Pure Appl. Geophys., 142: 491-527. Stirling, M.W., Wesnousky, S.G. and Shimazaki, K. (1996), Fault trace complexity, cumulative slip, and the shape of the magnitude-frequency distribution for strike slip faults: a global survey, Geophys. J. Int. 124, 833–868. Sylvester, A. G., and Smith, R, (1976), Tectonic transpretions and basement controlled deformation in the San Andreas Fault zone, Salton trough, California, AAPG Bull. 60, 2081-2102. Sylvester, A. G. (1988), Strike-slip faults, Geol. Soc. Am. 100, 1666-1703. Tchalenko, J.S. (1970), Similarities between shear zones of different magnitudes, Geological Society of America Bulletin 81, 1625–1640. Templeton, E. L., and Rice, J. R. (2008), Off-fault plasticity and earthquake rupture dynamics, 1. Dry materials or neglect of fluid pressure changes, J. Geophys. Res. Thatcher, W., 1983. Nonlinear strain build-up and earthquake cycle on the San Andreas fault, J. Geophys. Res., 88: 5893-5902. Turcotte, D. L., and Glasscoe, M.T. (2004), A damage model for the continuum rheology of the upper continental crust, Tectonophysics, 383, 71-80. Vermilye, J. M., C. H. Scholz, (1998), The process zone: a microstructural view of fault growth. Journal of Geophysical Research 103, 12223–12237. Vermilye, J.M., C.H. Scholz, (1999), Fault propagation and segmentation: insight from the microstructural examination of a small fault. J. Structural Geology 21, 1623–1636. Ward, S., (1996), A synthetic seismicity model for southern California: cycles, probabilities, hazard, J. Geophys. Res., 101: 22,393-22,418. Wdowinski, S., B. Smith-Konter, Y. Bock, and D. Sandwell, (2007), Diffuse interseismic deformation across the Pacific–North America plate boundary. Geology, 35,4; 311–314; doi: 10.1130/G22938A.1 Wesnousky, S. (1994), The Gutenberg-Richter or characteristic earthquake distribution, which is it?, Bull. Seismol. Soc. Am. 84, 1940–1959. Wesnousky, S. G. (2006), Predicting the endpoints of earthquake ruptures, Nature, 444, 358-360. Wilcox, R.E., Harding, T.P.,, Seely, D.R. (1973), Basic wrench tectonics, AAPG Bull., 57, 7496. Ziv, A. and A. M. Rubin, 2003. Implications of rate-and-state friction for properties of aftershock sequence: quasi-static inherently discrete simulations. J. Geophys. Res. 108 (B1): 2051.  13  2  STRUCTURAL PROPERTIES AND DEFORMATION PATTERNS OF EVOLVING STRIKE-SLIP FAULTS: NUMERICAL SIMULATIONS INCORPORATING DAMAGE RHEOLOGY1  Summary We present results on evolving geometrical and material properties of large strike-slip fault zones and associated deformation fields, using 3-D numerical simulations in a rheologically-layered model with a seismogenic upper crust governed by a continuum damage framework over a viscoelastic substrate. The damage healing parameters we employ are constrained using results of test models and geophysical observations of healing along active faults. The simulated fault zones form initially as complex segmented structures and evolve overall with continuing deformation toward contiguous, simpler structures. Along relativelystraight mature segments, the models produce flower structures with depth consisting of a broad damage zone in the top few kilometers of the crust and highly localized damage at depth. The flower structures form during an early evolutionary stage of the fault system (before a total offset of about 0.05 to 0.1 km has accumulated), and persist as continued deformation localizes further along narrow slip zones. The tectonic strain at seismogenic depths is concentrated along the highly damaged fault cores, although at shallow depths a small portion of the strain is accommodated over a broader region. This broader domain corresponds to shallow damage zones which have been identified in several seismic and geodetic studies of active faults. The models produce releasing stepovers between fault zone segments that are locations of ongoing interseismic deformation. Material within the fault stepovers remains damaged during the entire earthquake cycle (with significantly reduced rigidity and shear wave velocity) to depths of 10 to 15 km. These persistent damage zones should be detectable by geophysical imaging studies and could have important implications for earthquake dynamics and seismic hazard.  1  A version of this chapter has been published. Finzi, Y., Hearn, E.H., Ben-Zion, Y. and Lyakhovsky V. (2009) Structural properties and deformation patterns of evolving strike-slip faults: Numerical simulations incorporating damage rheology. Pure Appl. Geophys., 166, 1537-1573 doi: 10.1007/s00024-009-0522-1.  14  2.1  Introduction Understanding the geometrical and mechanical properties of fault zones is important for  many geo-science fields, including earthquake mechanics, crustal hydrology and mineral exploration. Since crustal faults generally grow and evolve through repeated earthquake ruptures, there are fundamental feedback mechanisms between the earthquakes sustained by a fault and its structural evolution. Observational and theoretical studies indicate that the temporal and frequency-size statistics of earthquakes change as faults evolve with cumulative slip, from disordered structures to more regular mature fault zones (e.g., Wesnousky, 1994; Ben-Zion, 1996; Stirling et al., 1996; Lyakhovsky et al., 2001). The evolution of the permeability structures around large faults affects the fluid flow properties of the crust and deposition of minerals (e.g., Micklethwaite and Cox, 2004). The evolving fluid flow regime influences in turn the mechanics of earthquakes and faults (e.g., Hickman et al., 1995 and references therein). Fault stepovers and other geometrical heterogeneities affect the initiation, propagation and termination of earthquakes (e.g., Sibson, 1985; King, 1986; Harris and Day, 1999; Oglesby et al., 2003; Wesnousky, 2006). Contrasts of elastic and permeability properties across faults affect the mode and properties of dynamic ruptures, seismic radiation and aseismic slip (e.g., Ben-Zion and Andrews, 1998; Rudnicki and Rice, 2006; Yamashita, 2007; Ampuero and Ben-Zion, 2008; Dunham and Rice, 2008). Detailed mapping of several exhumed fault zones (e.g., Chester et al., 1993; Evans et al., 2000; Sibson, 2003) and additional observations summarized by Ben-Zion and Sammis (2003) indicate that the internal structure of fault zones evolves from an early stage associated with distributed deformation and band-limited fractal structures at several hierarchies, through localization to principal slip zones, to a mature stage characterized by large-scale faults with tabular damage zones and narrow cores of ultra-cataclasite. However, the ranges of conditions over which such evolution takes place, and the coupling between the evolving structures and distributions of crustal stress and strain, are not well-understood. Various studies attempted to model changes of fault properties with ongoing deformation. For example, Olson and Pollard (1989) modeled the evolution of joints based on linear elastic fracture mechanics. Cowie et al. (1993) simulated evolving geometrical properties of fault networks using a scalar elastic field on 15  a lattice model with spring-like elements. Andrews (2005), Ben-Zion and Shi (2005) and Templeton and Rice (2008) simulated the generation of off-fault plastic strain during propagation of dynamic ruptures on frictional faults surrounded by a solid governed by Coulomb plastic yielding. While these studies provide important insights for various topics, they do not account for the evolution of elastic properties that accompanies the generation of cracking and inelastic strain, and they are also typically done within 2D “plane strain” frameworks. In the present work we attempt to understand some general aspects of the evolution of large strike-slip fault zone structures. The study is based on three-dimensional numerical simulations with a regional lithospheric model consisting of a seismogenic crust governed by damage rheology over a viscoelastic substrate (Ben-Zion and Lyakhovsky, 2006). Using this framework with parameters constrained by laboratory and geophysical observations, we examine the evolving geometrical and elastic properties of fault zones and the associated deformation patterns. In the next section we summarize observational results on fault zone structures that are relevant to our study. In section 2.3 we review the damage rheology framework and key aspects of the numerical model employed in this work. In section 2.4 we use geophysical observations of strength degradation and recovery within active fault zones to narrow the range of admissible damage rheology parameters. Section 2.5 contains the results of our parameter-space study on structural evolution of large strike-slip fault zones. Our simulations produce for ranges of realistic conditions flower structures with depth and secondary faulting within stepovers comparable to those documented in geological and seismic studies. The results support the view that fault zones display highly localized slip embedded within a wider shallow damage zone. The implications of the results to plate boundary dynamics and suggestions for continuing studies are discussed in section 2.6.  2.2  Geological and geophysical observations of fault zone structure The geometrical properties of fault structures and earthquake slip zones have been  documented in many geological and geophysical studies. In general, strike-slip fault zones display a nested hierarchy of damage zones and slip surfaces that form “flower structures” with 16  depth (e.g., Wilcox et al., 1973; Sylvester and Smith, 1976; Sylvester, 1988). In a typical fault structure, the principal slip zone is surrounded by gouge and embedded within a tabular or wedge-shaped damage zone (Ben-Zion and Sammis, 2003, and references therein). The extent of the damage zone may be defined as the region in which the density of deformation features exceeds the average density in the surrounding host rock (Chester, 1995). Studies of earthquake slip in exhumed faults and paleoseismic trenches indicate that within the top few kilometers of the crust the majority of coseismic slip is accommodated along very narrow slip zones (e.g., Sibson, 2003; Rockwell and Ben-Zion, 2007). The highly localized slip zone and surrounding ultracataclasite layer are referred to as the „„core‟‟ of the fault zone. This fault core is typically parallel to the macroscopic slip vector and is surrounded by a cataclasite layer which is a few meters thick (e.g., Chester and Chester, 1998; Schulz and Evans, 2000). The damage zone (DZ) around the fault core typically consist of a zone of intense damage, and possibly pulverized rocks, with width of a few hundred meters (Dor et al., 2006, 2008), which is surrounded by a broader, several kilometers wide zone, of distributed damage. The latter is probably a relic structure of the progressive coalescence and localization of the active fault zone over time (Ambraseys, 1970; Kim et al., 2004; Sibson, 2003). The near surface observations of distributed DZ are supported and complemented by a variety of geophysical studies that associate the DZ with a negative gravity anomaly and low seismic velocities (e.g., Stierman, 1984; Mooney and Ginzburg, 1986), along with anisotropic seismic waves (e.g., Cochran et al., 2003; Liu et al. 2004; Peng and Ben-Zion, 2004, 2005) and elevated seismic scattering (e.g., Revenaugh, 2000). Within the fault zone, seismic waves may be trapped in a narrow zone of intense coherent damage that is significantly distinct from the wider distributed damage zone (e.g., Ben-Zion and Aki, 1990; Li et al., 1990). Such fault zone trapped waves have been observed along large faults of the North Anatolian Fault System (NAFS), the Eastern California Shear Zone (ECSZ) and the San Andreas Fault System (e.g., Li et al., 1994; Ben-Zion et al., 2003). Systematic inversions of trapped waves indicate ~100 m wide tabular zones that extend typically to ~3-4 km depth and are characterized by strong attenuation and ~30-50% shear wave velocity (Vs) reduction relative to their surroundings (Peng et al., 2003; Korneev et al., 2003; Lewis et al., 2005).  17  Damage zones also show up in geodetic measurements that detect amplification of deformation signals along fault zones. Ben-Zion et al. (1990) observed amplified strain and water level signals at several locations along the Mojave segment of the San Andreas fault. Fialko et al. (2002) and Hamiel and Fialko (2007) interpreted measurements of Interferometeric Synthetic Aperture Radar (InSAR) along several large strike-slip faults in terms of compliant zones that are 1-2 kilometers wide, 3 km to over 10 km deep, and having rigidity () reduction of 50%-70% relative to the host rock. The discrepancies in dimensions (particularly, width) of the geodetically determined “compliant zones” and seismically determined “trapping structures” probably reflect differences between the broader long-term, quasi-passive, damage structure and the narrower active zone associated with recent earthquake ruptures. Other effects may feed into these discrepancies. For example, distributed microfractures affect the geodetically observed static strength more than they affect the seismically observed dynamic strength of rock (e.g., Ide, 1936; Eissa and Kazi, 1988). Other seismic observations such as seismic anisotropy (e.g., Cochran et al., 2003; Liu et al. 2004; Peng and Ben-Zion, 2004, 2005) and elevated scattering (e.g., Revenaugh, 2000) near large faults conform with the geodetically determined wide damage zones (1-2 km wide at the top few kilometers), indicating perhaps that the “trapping structures” are much smaller than the entire damage zone. To date, the best evidence of high localization of seismic slip at depths larger than 3-5 km is associated with the general tendency of seismicity to localize along relatively-straight fault segments to zones with width that is comparable to the smallest dimension that is resolvable by the data analysis. In places with good network coverage and relocated seismicity, the width of such zones is only a few tens of meters (e.g., Poupinet et al., 1984; Nadeau et al., 1994; Schaff et al., 2002; McGuire and Ben-Zion, 2005; Thurber et al., 2006). A significant deviation from the relatively simple DZ structure described above occurs at fault stepover zones. Fault zones generally display higher geometrical complexity and broader damage zones within stepovers than along relatively straight segments. While the major fault segments reflect a positive feedback of strain weakening and strain localization along the fault cores, persisting geometrical features such as fault offsets, kinks, and bends, can produce strain hardening that leads to local complexity and secondary fractures at different scales (Ben-Zion and Sammis, 2003). Many studies have characterized macroscopic structural features within 18  enlarged damage zones at geometrical irregularities (e.g., Segall and Pollard, 1980; Kim et al., 2004). Our study attempts to clarify the evolution of structural properties of fault zones along relatively straight segments as well as near large persisting stepovers.  2.3  Damage rheology framework  2.3.1 Theoretical background In the past decade continuum damage mechanics models have been successfully applied (e.g., Bercovici and Ricard, 2003; Turcotte and Glasscoe, 2004) in various studies of long-term tectonic deformation. Lyakhovsky et al. (1997a, b), Hamiel et al. (2004) and references therein developed a thermodynamically-based continuum damage model for evolving elastic properties of rocks sustaining irreversible brittle deformation. The employed damage rheology is applicable to volumes with a sufficiently large number of cracks that allow quantitative description through properties of the crack distribution rather than those of the individual cracks (Lyakhovsky and Myasnikov, 1984, 1985). The model generalizes the strain energy function of a solid to account for first-order macroscopic effects of distributed cracks (i.e., damage), and makes the elastic moduli functions of an evolving damage state variable  representing the local crack density. An undamaged material with =0 is the ideal solid governed by 3-D linear elasticity, while a material with =1 is completely destroyed. Using the balance equations of energy and entropy, the damage model quantifies the effective elastic behavior of a cracked solid for all intermediate values of the damage variable (0<<1). Detailed reviews and recent developments of the model can be found in Ben-Zion and Lyakhovsky (2006) and Lyakhovsky and Ben-Zion (2008). Here we only summarize the main ingredients of the model that are relevant for our work. The effects of distributed cracks (i.e., existing damage) on the elastic properties of a solid are accounted for in the damage model by generalizing the strain energy function to the form:  U  1  2   I1  I 2  I1 I 2  2   (2.1)  where I1=kk and I2=ijij are the first and second invariants of the elastic strain tensor ij,  is the mass density,  and  are the Lamé parameters, and is a third modulus of a damaged solid. The 19  first two terms of equation (2.1) give the classical strain potential of linear elasticity (e.g., Malvern, 1969). The third term may be derived using the effective medium theory of Budiansky and O'Connell (1976) for non-interacting cracks that dilate and contract in response to tension and compression (Lyakhovsky et al., 1997b), or by expanding the strain energy potential as a general second-order function of I1 and I2 and eliminating non-physical terms (Ben-Zion and Lyakhovsky, 2006). The kinetic aspects of the damage rheology involve making the elastic moduli functions of the damage state variable, and deriving an equation for the evolution of damage. Lyakhovsky et al. (1997a) showed that the leading term of the damage evolution equation, satisfying energy conservation and non-negative entropy production, can be written as   C d I 2    0 , d   dt    C1  exp  I 2    0 ,   C2   where I1  for    0  (2.2) for    0  I2 is referred to as the strain-invariants ratio, and Cd, C1, C2 are damage rate  parameters further described in sections 2.3.2 and 2.4.1. The parameter 0 is a yielding threshold separating states of deformation involving material degradation (d/dt > 0) when    0 , and material healing (d/dt < 0) when  < 0. Agnon and Lyakhovsky (1995) and Lyakhovsky et al. (1997a) related this parameter to the angle of internal friction by considering the critical shear stress for Mohr-Coulomb sliding. They obtained 0 = −0.8 for rock with internal friction coefficient of f = 0.6 and Poisson‟s ratio = 0.25, and noted that this value varies only slightly (0.7 to -0.9) for rocks with Poisson ratios between 0.2 and 0.3. Equation (2) was derived assuming for simplicity   = constant;    m   (2.3)  = m ; 20  where m is the maximum value of the third elastic modulus defined by normalization of the damage variable. The dependencies of elastic moduli on the damage variable produce the following changes during loading: as the damage variable  increases, the shear modulus  decreases, Poisson ratio  increases, and the modulus  increases from 0 (damage free) to m. Following the onset of positive damage evolution above the elastic limit at    0 and before the final macroscopic failure, the model incorporates a gradual accumulation of inelastic strain, ijv, described in Appendix A. When the damage variable reaches a critical value cr, there is brittle instability leading to rapid conversion of deviatoric elastic strain to permanent plastic strain. The reduced deviatoric stress at the end of the brittle failure episode typically leads to a state    0 that is associated (equation 2.2) with healing. The exponential dependency of the damage recovery (healing) is motivated by the logarithmic healing with time that is observed for rocks and other materials (e.g., Dieterich, 1978; 1979). Lyakhovsky et al. (2005) showed that the above damage model reproduces the main observed features of rate- and state-dependent friction, and constrained the healing parameters C1, C2 by comparing model calculations with lab frictional data. The main components of the numerical procedure, utilizing the Fast Lagrangian Analysis of Continua (FLAC) algorithm (e.g., Cundall and Board, 1988; Poliakov et al., 1993), are presented in Appendix A. To simulate long-term deformation processes with appropriate boundary conditions at the edges of our model domain, we use boundary conditions with variable forces (Lyakhovsky and Ben-Zion, 2008). These boundary conditions account for the stress buildup and abrupt drop during each seismic cycle, and for the evolution of elastic properties and cumulative plastic strain within the model domain. Appendix B presents an overview of the variable-force boundary conditions, and Appendix C describes the viscoelastic rheology applied in our models and applies the variable force boundary conditions in a test study verifying the viscoelastic component of our code. Additional details on the employed damage model and comparisons of results with laboratory fracture and friction data are given by Lyakhovsky et al. (1997a,b, 2005), Hamiel et al. (2004, 2006), Ben-Zion and Lyakhovsky (2006) and Lyakhovsky and Ben-Zion (2008).  21  2.3.2 Damage model parameters A fundamental set of results of previous damage-based models is that rheological damage parameters have significant impact on the evolving geometrical properties of fault zones, seismicity patterns and spatial distribution of deformation. Lyakhovsky et al. (2001) and BenZion and Lyakhovsky (2006) suggested that damage zone structure is primarily controlled by (1) the ratio between loading rate and healing rate, (2) the overall degree of “brittleness” of crustal deformation which may be parameterized by the seismic coupling coefficient , and (3) the susceptibility to propagation of rupture associated with dynamic weakening and related dynamic time scale r. Below we review the main material parameters and outline their effect on fault zone evolution. Table 2.1 presents a synthesis of the plausible range of values for each parameter, and justifications for these values. As discussed in the context of equation (2.2), the employed damage rheology includes three material parameters that affect the rate of damage evolution with time: C1 and C2 are both healing rate parameters, and Cd is the damage accumulation rate parameter. In our simulations, the loading rate is closely linked to the specified tectonic strain rate and the value of Cd is fixed. Therefore the ratio of loading rate to healing rate is governed by the healing rate parameters C1 and C2. Simulations with high healing rates compared to the loading rate result in rapid nearcomplete healing of fault damage. In such settings, ruptured fault zones quickly regain their strength, enabling larger interseismic stress accumulation and coseismic stress drops, and the evolving fault zones have more complex geometries than cases with low healing rates. The degree to which crustal deformation is brittle is controlled by the material parameter Cv, the coefficient of damage-related inelastic deformation (Appendix A). This material parameter determines the ratio of aseismic to seismic components of deformation, or the seismic coupling coefficient, given (Ben-Zion and Lyakhovsky, 2006) by  = 1/(1+R) through the nondimensional R-value with R=0Cv and 0 being the initial rigidity of the material. Lyakhovsky and Ben-Zion (2008) showed that higher crustal Cv values induce larger components of aseismic deformation in the seismogenic zone, and therefore lower coseismic stress-drops. The susceptibility to rupture propagation determines the degree of dynamic weakening and dynamic 22  time scale r during the occurrence of brittle instability. Lyakhovsky et al. (2001) found that higher r values induce larger ruptures and lead to relatively simple failure histories consisting of system-sized events occurring in a single fault zone.  Parameter  Preferred range  Cd  Damage accumulation rate  C1 , C2  Healing rate parameters  0.5 - 5 s-1  Comments  Lyakhovsky et al. (1997a): Fracture experiments with granite at relatively high confining pressure (100 MPa)  Lyakhovsky et al. (2005) show that at shallow depths (z<-5 km) Cd increases with the decrease in confining pressure (Cdsurface >10 s-1)  Lyakhovsky et al. (2005):  C1= 10-24-10-4 s-1 Analysis of 1D damage and See discussion and new constraints C2=0.1- 0.01 comparison with rate and state in section 2.4 dependent friction parameters  Agnon & Lyakhovsky (1995): Analytic analysis indicates that a Lyakhovsky et al. (1997a) also friction coefficient of f~0.6-0.7 indicate a week dependency of 0 on  o  Critical strain invariant ratio (generalized internal friction)  Justification for value, and references  -0.7 to -1  corresponds to 0=-0.8 Poisson ratio, , but concluded that Lyakhovsky et al. (1997a): 0=-0.8 for various rocks with 3D faulting experiments yield =0.2-0.3 0= -0.7 to -1  Cv  Yang & Ben-Zion (2008): Analysis of aftershock sequence in southern California, and comparison to damage rheology predictions  r  Ben-Zion & Lyakhovsky (2006): Largest aftershock magnitude analysis (comparing numeric results and expected log-linear relation)  damage-related 10-4-5·10-6MPa-1 inelastic strain accumulation  Characteristic time scale for seismic wave damping  3·102-3·104s  1. Cv· d/dt is the damage related compliance 2. Based on these Cv values, the fraction  of elastic strain released during a seismic cycle as brittle deformation is estimated to be 30%85%  Table 2.1 Damage rheology parameters and their constraints.  23  The above material parameters are currently constrained mainly by analytical considerations and by fracture and friction experiments (Table 2.1). In the next section, we use geophysical data to narrow the admissible range of healing rate parameters, and to better relate these parameters to natural deformation processes observed along active fault systems.  2.4  Geophysical constraints on healing rate parameters  2.4.1 Parameters C1 and C2, and healing as a function of time As follows from equation (2.2), the healing rate depends on the strain magnitude (I2), strain invariant ratio (), level of material damage ( and material properties including critical strain invariant ratio ( and two healing rate parameters C1 and C2. The critical strain invariant ratio is well constrained and fairly constant (  -0.8). The strain invariants ratio () varies  from  3 to  3 . The post failure shear strain (I2) is mostly controlled by the lithostatic pressure and may vary with depth by two orders of magnitude or less. Therefore, the rate and the overall effectiveness of the healing process are primarily determined by the rate parameters C1 and C2. While C2 determines the rate dependence on the damage state  and varies within one order of magnitude, C1 may vary by many orders of magnitude (Lyakhovsky et al., 2005). Depending on the combination of these parameters the healing process may be fast or slow, and may yield insignificant or near-complete healing of the damaged material over long timescales. Equation (2.2) indicates that very large C1 corresponds essentially to zero memory, in which case damaged material heals rapidly and almost completely. If damage accumulation is also rapid (Cd is very large) the model will display ideal elastoplastic behavior. Very small C1 yields insignificant healing except for the special case of C2  0 in which healing is near instantaneous  and complete. Extremely large C2 values lead to a healing rate proportional to C1 (i.e., d/dt = C1I2(0)). In such cases, for C1 < 1∙10-10 s-1 the healing is slow and insignificant while for C1 > 1∙10-05 s-1 the healing is fast and almost complete (Fig. 2.1). To better understand the healing process and constrain the parameters C1 and C2, we define a time scale for healing (h) during which the relative change of the elastic moduli () is above an arbitrarily chosen rate of 0.1% yr-1 (d/dt =3 ∙10-11sec-1). Simplifying equation (2.2) 24  for a uniform strain invariant ratio suitable for healing, assuming ( 0) to be of the order of one, and setting the healing rate to this chosen threshold, the expected damage level is given by         f = C2 ln ( 3∙10-11/ C1I2),  (2.4)  where f is the damage level as the healing becomes slower than 0.1% yr-1. This value is referred to hereafter as the final damage level, even though slow but possibly significant healing continues on after f is reached (at geological time scales). Substituting f in equation (8) of Lyakhovsky et al. (2005) for the damage level as a function of time, the time scale for healing (i.e., the duration required to reduce the damage level from = 1 to f) can be estimated by  1  f    1 exp C 2   h  C  1  I 2 1 exp  C2  C2   (2.5)  Within the healing parameter-space that represents materials that undergo significant healing (i.e., damaged materials with f < 0.8) the time scale for healing (h) falls in the range 10-110 yr (Fig. 2.1a). The time scale for healing indicates how fast a damaged material heals to a near constant damage level, but it does not reveal the final damage level. Therefore, materials with similar h may exhibit a wide range of final damage levels (Fig. 2.1). Figure 2.1b displays six healing processes determined by six sets of C1 and C2 parameters (marked and labeled in Fig. 2.1a). Healing processes with short time scales (e.g., 1a and 1b in Fig. 2.1b) display higher initial healing rate and faster decay of the healing rate compared to processes with longer time scales (e.g., 3a and 3b in Fig. 2.1b). As evident from Fig. 2.1b, in order to fully determine the healing parameters of a material one would either need to establish the damage level at two distinct times after failure or to determine the current damage level, the duration since failure and the healing time scale relevant to the specific healing process. Lyakhovsky et al. (2005) suggested that the parameter C2 is closely related to the parameter b of rate and state friction (b~10-1; C2~10-2-10-1), and obtained the following relation between parameters C1 and C2:  25     2 C1  BC2 exp  0  /  cmp  C2   (2.6)  where B (~1-2 s-1) is a laboratory-determined time scale for the evolution of static friction with hold time (Dieterich, 1972, 1978) and cmp is the compaction strain estimated by the ratio between lithostatic stress and the bulk modulus (K). Lyakhovsky et al. (2005) estimated cmp ~102  for crustal rocks at seismogenic depths, but noted that this strain level may vary significantly  for various lithologies and depths. Since the lithospheric stress within the seismogenic zone (depth 1-20 km) ranges between 20 MPa and 400 MPa, and the bulk modulus of typical crustal rocks varies by an order of magnitude (Christensen and Mooney, 1995), it is safe to assume that the compaction strain could vary by 2-3 orders of magnitude. Given that the parameter B may differ from the well-constrained lab-based values and that compaction strain may vary significantly, we consider (Fig. 2.2) a wide range of C1 values per C2 value (six orders of magnitude).  26  Figure 2.1 (A) Damage levels of typical crustal material at 10 km depth after 50 years of healing, for the damage model parameter space suggested by Lyakhovsky et al., 2005. Black and red lines show contours of equal damage and final damage (f), respectively. White lines are contours of equal healing time scale (h). Black diamond symbols indicate sets of healing parameters referred to in Figure 2.1B. (B) Damage as a function of time after failure for six sets of C1 and C2 parameters (indicated in Figure 2.1A). The healing time scale (h) and final damage level (f) of each healing process are indicated by diamonds (red – 1a,b; black- 2a,b; blue – 3a,b). Note that some healing occurs after f is attained (see section 2.4.1).  27  2.4.2 In-Situ geophysical constraints for healing parameters C1 and C2 We use data from seismic surveys along large fault zones with significant fault-related damage to better constrain the healing parameters C1 and C2. Simplifying equation (2.3) for uniform shear deformation = 0(1-), and using the relation between rigidity, density and shear wave velocity (Vs2), we convert reported seismic velocity and rigidity reductions (section 2.2) into damage level estimates. The seismic and geodetic studies indicate that major strike-slip fault zones rapidly heal in the top few km to   0.75 during the early postseismic  stage (Ben-Zion et al., 2003; Peng et al., 2003; Lewis et al., 2005), and thereafter display damage levels of  > 0.5 (Hamiel and Fialko, 2007; Fialko et al., 2002; Fialko, 2004). While the healing at greater depth is expected to be higher, these observations may indicate that damage zones of large active faults do not completely heal over time scales of typical earthquake cycles. This argument is supported by the abundance of ancient fault zones that remain weaker than the surrounding rock (Tchalenko, 1970; Sengor et al., 2005; Armijo et al., 1996; Powell and Weldon, 1992; Evans et al., 2000). Additional support comes from previous numerical models (Lyakhovsky et al., 2001; Finzi et al., 2006) and experimental work (Tenthorey et al., 2003) showing that damage zones do not heal completely during the earthquake cycle. The above is expected to be valid for low-porosity crystalline rocks. In contrast, deformation bands in sandstones and other high porosity rocks are frequently denser and stronger than their host rock (Aydin and Johnson, 1983; Shipton and Cowie, 2003). The limited healing argument implies that long-term interseismic healing in low porosity rocks is typically minor and that damage generation and healing in such rocks occurs predominantly in the seismogenic crust during the co- and early post-seismic interval (e.g., over weeks to months). This is supported by various postseismic healing rate studies (e.g., Li et al., 2006; Karabulut and Bouchon, 2007; Peng and Ben-Zion, 2006; Schaff and Beroza, 2004; Rubenstein and Beroza, 2004; Baisch and Bokelmann, 2001; Wu et al., 2009). We note that the similar interseismic damage levels for the different fault zones mentioned above may reflect resolution limitations of the seismic and geodetic methods (i.e., perhaps materials with  < 0.5 are not reliably detected by these techniques). If this is the case, the argument that long term healing is minor may not be valid. Furthermore, the healing parameters may be pressure and/or temperature dependent (as the damage accumulation rate 28  parameter Cd). Further work should be done to better constrain the healing parameters at seismogenic depths (e.g. > 3 km). Finally, the healing computations in our parameter space study do not take into account ongoing deformation (and damage accumulation) due to continuous tectonic loading or nearby earthquakes, and therefore they may underestimate damage levels in natural systems. Based on the above considerations we suggest two general constraints for the healing parameters suitable for models of natural processes: (a) the minimum damage level expected in shallow crustal fault zones during the interseismic stage should be above   0.4, and (b) the  healing time scale representative of natural damage zones should be shorter than h  40 yrs  (yielding healing rates of 5-10% yr-1 after 4-5 months of healing, and a very low rate of approximately 0.1 % yr-1 after 40 years of healing). The resulting healing parameter-space is outlined in Fig. 2.2 by a hatched pattern. The admissible values of C1 and C2 in that subspace are 10-24 s-1 to 10-10 s-1 and 0.015 to 0.035, respectively. A final analytical constraint for healing parameters can be derived from the convexity condition for macroscopic failure used in our damage rheology framework (Lyakhovsky et al., 1997a). This stability condition indicates that near the surface, where normal stress is low compared to shear stress and the strain invariants ratio is approximately zero, the maximum sustainable (stable) damage level is approximately  = 0.75. Therefore, healing parameters that yield lasting damage levels greater than  = 0.75 at shallow depths (z < 3 km) are assumed to be non-realistic (Fig. 2.2).  29  Figure 2.2 (A) Geophysical, analytical and laboratory-based constraints on healing parameters. Damage levels at a depth of 1 km after 50 years of healing are shown for C1 values within three orders of magnitude of their expected values (Eq. 2.6). The hatched region bracketed by the maximum healing time scale (h=40 yr) and the minimum and maximum admissible damage levels (0.75> f >0.4) represents the healing parameters relevant to models of natural damage zones (see text). Diamond symbols indicate two healing parameter sets plotted in B for illustration. (B) Estimated damage versus depth after 0.1 yr (gray lines), 10 yr (dark gray lines) and 100 years of healing (thick black lines) under lithospheric stress conditions, for healing parameters representative of natural fault zones (solid lines: C1=1∙10-18s-1, C2=0.02; dashed lines: C1=1∙10-12s-1, C2=0.03). 30  2.5  Damage and strain distribution across active strike-slip faults To investigate the structure of damage zones associated with evolving strike-slip fault  systems, we use several realizations of a three-dimensional model of transform plate boundary incorporating damage rheology in the upper crust.  2.5.1 Model setup A typical model setup (Fig. 2.3) consists of a layered seismogenic crust governed by damage rheology, underlain by viscoelastic lower crust and upper mantle layers. The modeled region is 100-250 km in the along-strike direction, 100 km wide and 50 km deep. Detailed description of such a model setup is given by Ben-Zion and Lyakhovsky (2006). Here we only summarize the main ingredients. A diabase flow law is used to represent the rheology of the crystalline crust (Carter and Tsenn, 1987) and a dunite/olivine flow law is used for the upper mantle (Kirby and Kronenberg, 1987). We assume a geothermal gradient of 20 C km-1. The assumed flow laws and geotherm are kept fixed in our simulations. A range of damage model parameters, chosen from the values given in Table 2.1, is used in the models (Table 2.2). Since we are not trying to characterize surface damage structures, and as frequent failure of surface elements due to low confining stress would be computationally time consuming, we suppress damage accumulation in the simulated surface layer (typically top 3 km of the crust) by setting there Cd = 0. However, based on test models with Cd > 0 within the surface layer (not presented here), we expect that the surface damage zone is slightly wider and consists of locally higher damage levels than the underlying (simulated) damage zone. A variable force boundary condition (Lyakhovsky and Ben Zion, 2008, 2009) is applied to the sides and bottom of the model domain, simulating a constant far field fault-parallel velocity with relative rate of 32 mm/yr (corresponding to the San Andreas Fault). These boundary conditions are further discussed in Appendix B. The boundary driving forces are set to represent imposed fault zones outside the model domain, and they induce faulting near the centers of the north and south edges of the model (see for illustration Fig. B1). The top model boundary is stress free.  31  Model name  Grid spacing (km)  Healing C1 (s-1)  Healing C2  Dynamic weakening r  Seismic ratio   Initial damage heterogeneity  Initial r heterogeneity  NB_0 NB_2 NB_6 NB_7 Nb_1_lap Nb_3_lap Nb_5_lap Nb_7_lap Nb_8_lap Nb_9_lap Nb_11_lap Nb_geos_4 Nb_geos_6 Prop_2 Prop_7 Prop_9 Prop_lap_2 Prop_lap_6 Prop_run_a9 Prop_geos_5 Prop_lin_4.5 Stepover_8 Stepover_9 Stepover_10 Stepover_11 long-term_1  3.2 3.2 2.2 2.2 3.2 3 2.2 2.2 2.2 2.2 2.2 1.6 1.6 3.2 3.2 3.2 3.2 3.2 2.2 2.2 4.5 0.25 0.25 0.25 0.25 0.6  6.0E-11 1.0E-13 1.0E-13 1.0E-13 1.0E-13 1.0E-13 1.0E-13 1.0E-13 1.0E-13 1.0E-13 5.0E-11 1.0E-13 1.0E-13 1.0E-13 2.0E-09 6.0E-11 1.0E-13 1.0E-13 1.0E-12 1.0E-12 1.0E-13 1.0E-10 1.0E-10 1.0E-10 1.0E-10 1.0E-12  0.07 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.05 0.02 0.02 0.03 0.15 0.07 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.025 0.025 0.03  1.0E+4 1.0E+4 5.0E+3 5.0E+3 1.0E+4 1.0E+4 1.0E+4 3.0E+3 3.0E+3 5.0E+3 4.0E+3 1.0E+4 2.0E+3 8.0E+3 1.0E+4 1.0E+4 9.0E+3 1.0E+4 5.0E+3 7.0E+3 1.0E+4 3.0E+1 1.0E+2 4.0E+2 6.0E+2 6.0E+4  80% 80% 61% 99% 80% 80% 80% 67% 97% 57% 61% 80% 67% 80% 80% 80% 80% 80% 61% 72% 80% 40% 40% 40% 40% 80%  0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 30% 30% 25% 15% 25% 25% 10% 10% 10% 0% 100% 100% 50% 100% 25%  30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30% 30%  Table 2.2 Model parameters used in our fault evolution study. Other material parameters were set to represent crustal materials and were not varied in our models (these parameters include:  0=-0.8. Cd= 5 s-1; Sedimentary layer: density =2.4∙10-3 Kg m-3, Newtonian viscosity =1019 Pa S; Crustal rheology: Young's modulus E=80 GPa, Poisson's ratio = 0.3, =2.8∙10-3 Kg m-3, viscosity flow law coefficients: A=6.3∙10-20 Pa-n S-1, n=3.05, Q=276 Kj mol-1; Mantle rheology: E=150 GPa, =0.3, =3.3∙10-3 Kg m-3, A=7.0∙10-14 Pa-n S-1, n=3.0, Q=520 Kj mol-1). 32  Figure 2.3 A block diagram of typical 3-D lithospheric structure used in the numerical simulations. The fault parallel extent of the model domain varies from 100 to 250 km in different simulations. Imposed damage (not shown in Figure 2.3) is applied only in a few simulations of long term fault stepover evolution (i.e. initial conditions of high-resolution fault stepover models included damage zones representing a segmented fault).  2.5.2 Model output: examples and interpretation The model outputs include the level of damage  and strain . We calculate surface velocities, rigidity  and other related quantities from these variables. Since ,  and  are computed throughout the model domain, both plan views (including depth slices) and cross sectional views of these parameters may be plotted at any time step. Figures 2.4 and 2.5 show examples of model outputs, and illustrate some features which correspond to observed geological structures such as fault segments, stepovers, and flower structures. 33  Contiguous sets of elements that fail repeatedly, resulting in a higher level of  (and a lower ) than their surroundings, are interpreted as fault segments (Fig. 2.4 and 2.5). Because of their relative weakness, these fault segments are also the centers of high velocity gradients and high strain rate (Fig. 2.4 and 2.5). Cross-sectional profiles through modeled fault segments (Fig. 2.5) display “flower structures” with depth, which comprise of localized damage along the active fault core with a superimposed, broader zone of distributed damage in the top 3-10 kilometers of the crust (Fig. 2.4 and 2.5). Based on these observations, we define two damage sub-zones that are distinct in their evolution patterns, damage level, and spatial distribution: (1) Localized Active-Fault (LAF) damage which represents the highly localized damage along the active fault cores (Fig. 2.4 and 2.5). The LAF damage is coseismically very high along the primary slip zone, but it rapidly heals at depth. (2) Distributed Off-Fault (DOF) damage which is sustained cumulative damage resulting from many earthquakes. The DOF damage develops during the early stages of fault-system evolution, and thereafter its spatial extent is stable and the degree of damage within it evolves locally (Fig. 2.4, 2.5 and 2.6). Descriptive analyses of damage structures along simulated strike-slip fault segments are given in section 2.5.3. Plan views of the model domain show several examples of stepover zones where segments are offset from one another (Fig. 2.4). Our models produce just extensional stepovers; that is, during propagation, new offset faults form in areas where end-effects from existing faults contribute a tensional mean normal stress. Stepovers are characterized by a wide DOF damage zone, and by high levels of damage, high strains, and low rigidity within the stepover. Descriptive analyses of damage structures within simulated stepover-zones are given in section 2.5.4.  34  Figure 2.4 Plan views of a segmented strike slip fault at several depths, showing damage level (), rigidity () and the second invariant of the deviatoric strain, Sd (Sd= sqrt(eij ∙eij) where eij =  ij-ij kk /3 and ij is the Kronecker delta). Shallow (z= 3 km) damage (top left panel) is distributed within the stepover and around the fault segments. At depth (z= 5-10 km) damage is highly localized along the fault segments and is distributed within the fault stepover. At the lower part of the seismogenic crust, damage within stepover may persist long after the localized damage along fault segments heals. 35  (b)  (a)  (c)  Figure 2.5 Cross-sections of a typical “flower type” damage zone along a strike-slip fault, displaying damage levels and rigidity (A and B), and fault-parallel velocity (C). Annotations in (A) show the dimensions of the Localized Active Fault (LAF) damage (black lines) and Distributed Off Fault (DOF) damage (white lines) as they were measured in this study. Dashed black and white lines show the location of the LAF and DOF damage zones (respectively) in the smoothed contour plots (B and C).  36  Before we discuss the results of our models, we need to (1) introduce observable quantities which may be systematically measured and then used to compare between model results, (2) define a threshold criterion for the maturity of our modeled faults to ensure we base our analysis only on models that were run long enough to form mature damage zones, and (3) confirm that the model results presented here are fairly insensitive to numerical element dimensions. This is necessary because the models shown on Table 2.2 were run for different durations, and for a range of element dimensions. To systematically describe the spatial extent of damage zones around faults we have chosen a threshold damage level of = 0.35 (presumably above any expected background damage level). The four quantities used to define the extent of damage are the widths and depths of both the LAF and DOF damage zones. Measurements of these quantities are performed on plan-view and cross-section plots of the simulated fault zones (without smoothing). The width of the DOF damage represents the maximal spatial extent of the damage zone, and its depth is the average depth extent of the shallow distributed damage away from the active fault core (Fig. 2.5a). The width of the LAF damage represents the maximum width of the localized damage along the fault core, and its depth is the maximum depth extent of the damage zone (Fig. 2.5a). To facilitate comparison with geodetic studies we define also the Fault Compliant Zone (FCZ) in our models as the volume in which the average material rigidity is reduced by 50% relative to the host rock, consistent with the compliant zones of Fialko et al. (2002). The simulated FCZ typically consists of the entire LAF damage and most of the DOF damage. Strike-slip fault systems evolve over time, first becoming complex and then gradually simplifying to a more continuous configuration with fewer fault segments (e.g., Ben-Zion and Sammis, 2003 and references therein; Lyakhovsky and Ben-Zion, 2009). This complicates directly comparing model runs that may have been cut off at different evolutionary stages. Figure 2.6 illustrates that the width and depth of the DOF damage initially grows rapidly, starting to stabilize after a total relative displacement of approximately 0.05 km. This simulated stage of rapid DOF damage growth represents the initial stage of fault growth and complexity increase after which the fault configuration starts to stabilize and the strain localizes along the fault. Hence to analyze mature damage zone structures we can only use fault simulations with a total displacement larger than 0.05 km. 37  All results displayed on subsequent plots are for models with mature damage zones in which at least 0.05 km of displacement has accrued (corresponding to a modeled time interval of about 1600 years).  Figure 2.6 DOF damage zone width and depth plotted against cumulative strike-slip offset. After an initial stage with relatively fast damage zone growth, the DOF damage zone dimensions remain fairly constant (at offsets exceeding 0.05 km).  Since brittle failure in our simulations is associated with an abrupt transition from initial (static) to final (dynamic) stress levels, the numerical models are inherently discrete (e.g., Rice and Ben-Zion, 1996) and some aspects of the results are expected to be grid-dependent. An analysis of our entire set of damage zone simulations indicates that while the spatial extent of damage is somewhat grid dependent, the average level of damage within voluminous damage zones (e.g., stepover zones) is probably not grid-size dependent (Fig. 2.7). Furthermore, the spatial extent of damage zones in our models with element dimensions between 0.6 and 4 km (Fig. 2.7) are not significantly sensitive to element size. Results of models within this range of element dimensions will be shown together on subsequent plots. Models with finer and coarser elements were also run. For models with 4.5 km elements, results were smeared forming 38  apparently wider and deeper damage structures. Models with 0.25 km elements ran very slowly and were numerically unstable. These simulations were terminated after 10-20 days of CPU time, during which the simulated damage zones did not reach a stable width. Based on Figure 2.6 we interpret the narrow and shallow damage zones that formed in these simulations as being immature, and we therefore do not incorporate them in our analysis.  Figure 2.7 Sensitivity of model results to element size. (A) For models with element dimensions between 0.6 and 3.2 km, DOF damage zone dimensions are fairly insensitive to grid size. (B) Damage levels in fault stepovers at depth of 5-6 km and 10-12 km. These damage levels are insensitive to element size within the range we model. The scatter in results on both panels is due to different values of damage model parameters and variations in cumulative slip. 39  2.5.3 Damage, rigidity and strain distribution across strike-slip fault segments Simulated damage zones along strike-slip faults form flower structures consisting of a shallow DOF damage zone 6 to 14 kilometers wide and up to 7 kilometers deep, and a more intense LAF damage zone around the fault core. The DOF damage is the result of cumulative damage from past earthquakes (Fig. 2.5 and 2.6). The LAF damage is narrow and steeply dipping, and it may extend to the bottom of the seismogenic zone (Fig. 2.4 and 2.5) coseismically and in the early postseismic interval. The deeper parts of the LAF damage zone are only one model element wide and appear discontinuous (Fig. 2.4 and 2.5). This probably indicates a tendency to evolve to a much narrower extent than our model element dimensions, in agreement with observed extreme localization of active slip zones (e.g., Chester and Chester, 1998; Sibson, 2003; Rockwell and Ben-Zion, 2007). The highest damage levels along fault segments are usually found at the shallowest part of the LAF damage zone. The parts of the LAF damage in the shallow crust that are more than one element thick may correspond to the observed distributed DZ along the surface trace of large strike-slip faults. The flower structure in our simulations is a robust feature that shows little dependence on damage parameters (C1, C2, , r), and limited sensitivity to the presence of various degrees of material heterogeneities (Table 2.2). The insensitivity of the DOF damage to healing parameters (Fig. 2.8) may imply that the shallow crust is readily damaged and it experiences limited healing regardless of material parameters. This apparent insensitivity to the healing parameters may also indicate that the current range of modeled healing parameters is insufficient and further analysis is needed to understand the role of healing in damage zone evolution. The deeper sections of the LAF damage are more sensitive to the healing rate parameters (Fig. 2.9a). This reflects the fact that healing processes are sensitive functions of the normal stress. Our results indicate that the depth extent of the LAF damage (after the early postseismic interval) ranges from the entire seismogenic zone for materials with extremely slow healing (C2 > 0.05; h > 55 yrs; Fig. 2.9b) to a few kilometers for the faster healing materials (C2 < 0.03; h < 25 yrs; Fig. 2.9c). Model realizations with healing parameters representative of natural processes (see section 2.4) indicate that during most of the seismic cycle the contiguous well-developed LAF damage is limited to the top section of the crust (e.g., Fig. 2.9c and Fig. 2.10a). These results are consistent with numerical simulations of plastic strain generation during dynamic rupture (Ben-Zion and Shi, 40  2005), and analyses of large seismic data sets recorded around active faults (e.g., Ben-Zion et al., 2003; Peng et al., 2003; Korneev et al., 2003; Cochran et al., 2003; Lewis et al., 2005; Graymer et al., 2007).  Figure 2.8 Width and depth of the distributed part of the damage zone (DOF damage) for models with a wide range of material parameters. DOF damage zone dimensions are insensitive to the healing parameters and the seismic coupling ratio. 41  Figure 2.9 (A) Width and depth of the LAF damage zone as a function of healing parameter C2 (width measured at 5-8 km depth, just below the DOF damage). The hatched region indicates admissible values of C2 for modeling natural processes (see section 2.4). (B,C) Examples of deep and shallow fault core damage zones in models with long (B) and short (C) healing time scales. The geometry of these zones is stable through most of the interseismic interval.  The elastic strength of simulated fault zones is reduced as the damage level increases. The relation between rigidity and damage level (equation 2.2) implies that the effective rigidity is primarily a function of . However, the effective rigidity is further reduced near the surface and along the fault core because the strain invariants ratio  is higher at these localities. Our 42  damage zone models indicate that significant rigidity reduction (50-70% reduction) occurs within the shallow DOF damage and along the upper part of the fault-core (e.g., top 5 km of the LAF damage). Along the LAF at depths exceeding 12 km healing is rapid and the long-lasting (interseismic) reduction on  is relatively small (Fig. 2.4 and 2.5). Our fault evolution models indicate that the bulk of inelastic strain is concentrated in the highly damaged cores of the fault zones (i.e., along the LAF damage). However, our models also exhibit low-gradient strain beyond the fault cores in the uppermost crust (top 3-5 km). The modeled off-fault strain is typically confined to a shallow layer approximately 10-15 km wide that exhibits significant DOF damage ( > 0.4: see Fig. 2.10). Where the shallow DOF damage is of lower intensity (< 0.4) our models do not display significant strain, suggesting that relatively thin layers of slightly damaged rock may not modify regional surface deformation patterns and may not be easily detectable by means of geodesy. At distances greater than 10 km from fully formed fault-zones, the total deviatoric strain is negligible and the long-term average strain-rate is uniform, indicating that the undamaged upper crust behaves kinematically as a rigid block. The correlation between damage and strain distribution confirms recent interpretations of observed surface deformation patterns above fault-related compliant zones (e.g., Fialko et al., 2002; Fialko, 2004). This correlation indicates that geodetically observed compliant zones are related to relatively high damage levels ( > 0.4), and suggests that parts of the damage zones ( < 0.4) may not be geodetically observable.  43  Figure 2.10. Correlation between damage and strain across simulated fault-zones. (A) Damage at depths 1, 5, 10 km and superimposed fault-parallel velocity. (B) Plots of damage level (solid curve, right ordinate), long-term (multi-cycle) fault parallel velocity (dashed curve, outer-left ordinate) and deviatoric strain (gray dotted curve, inner-left ordinate) for the same depths, extracted along the central section shown in (A) (at y= 65 km).  44  2.5.4 Fault stepovers While major fault segments display a positive feedback of strain weakening and localization along highly damaged fault cores, persisting geometrical features such as fault offsets, kinks, and bends, may display strain hardening and produce local complexity and new fractures at different scales (Ben-Zion and Sammis, 2003). In our models, segmented fault zones display continuous distributed seismic and aseismic deformation within fault stepover zones. This consists of aseismic strain and small earthquakes; however, moderate earthquakes (M L< 5.5) also occur occasionally. While much of the damage along fault segments heals during the interseismic stage, the damage level at stepovers remains persistently high (Fig. 2.11a, compare with Fig. 2.9c). Fault stepovers and segment termination zones undergo significant damage accumulation during the interseismic stage, typically displaying average  between 0.5 and 0.75 (measured at a depth of 5-10 km), depending on the healing rate parameters. The damage within stepovers extends to greater depths than along fault segments and in many simulations it reaches the bottom of the seismogenic zone (Figs. 2.11 and 2.4). Models with realistic healing parameters typically exhibit significant damage ( > 0.5) down to depths of 10-15 km (Fig. 2.11). An important implication of the permanently elevated damage level within fault stepovers is that these regions of reduced  (Fig. 2.4) affect rupture propagation and strong ground motion patterns.  45  Figure 2.11. Cross-sections of damage patterns within fault stepovers. These zones of fault complexity are interseismically active, retaining high damage levels throughout the entire seismic cycle. Significant differences in damage level at depth (i.e., at Z >10 km) arise from varying healing parameters: models with fast and efficient healing (A) yield lower damage levels than those with slow and inefficient healing (B). Based on our healing rate analysis (section 2.4), the faster healing models (e.g., A, and Fig. 2.9c) are better representatives of natural DZ processes. The healing rates in these models are identical to those in Figure 2.9b and 2.9c.  2.5.5 Fault system complexity as a function of time According to our results, the DOF damage dimensions depend on the maturity of the fault system. Simulations with a wide range of material properties and numerical characteristics (e.g., element dimensions and boundary conditions) indicate that the damage zone grows until the fault accumulates an offset of about 0.05-0.1 km (Fig. 2.6). During this early evolutionary stage, the fault system‟s complexity increases as additional segments nucleate and propagate, forming new damage zones. This stage culminates as deformation localizes along narrow slip-zones. As mentioned, an important exception to this evolution stage is the widening of the DZ at bends and stepovers and other local fault complexities (Fig. 2.12a-c and Fig. 2.13). The fault system evolution continues as through-going faults bridge sites of fault complexity such as stepovers 46  (Fig. 2.12d and Fig. 2.13). When through-going faults fully form and bridge the entire depth of the stepovers, these structures may become inactive. After this stage the width and depth of distributed damage remain fairly constant, until the existing fault configuration cannot accommodate the evolving regional stress, at which point new faults form and migration of faulting may occur. Parallel faults may sustain simultaneous damage accumulation (Fig. 2.13) or exhibit alternating deformation periods before the initial faults become inactive and the new fault configuration stabilizes. These results are consistent with existing multi-disciplinary knowledge of fault systems (Ben-Zion and Sammis, 2003 and references therein; Dolan et al., 2007) and the related numerical results of Lyakhovsky and Ben-Zion (2009).  47  D  Figure 2.12 Fault stepover evolutionary stages displayed in four snap-shots of damage levels around a stepover zone (at 3 km depth). (A) Segmented fault. (B, C) The stepover zone displaying extensive damage. In (C) distinct lateral (subsidiary) “faults” (regions of high damage) link between the two fault segments. (D) Formation of a through-going fault through the entire stepover zone. These results are from high-resolution, small-domain models focusing on stepovers.  48  Figure 2.13 Fault evolution snap-shots showing damage levels at 10 km depth. The originally segmented fault (left panel; notice large stepovers) is smoothened with time (right, mature fault 1; note that some of the apparent small stepovers are actually numerical artifacts formed because the fault is not parallel to the grid). An additional fault formed after approximately 5000 yrs to better accommodate regional stress (right, fault 2). Annotations of times in the images indicate the age of each fault in the simulation. Finely dashed line long fault 1 in both panels indicates its original segmented outline (as shown in left panel).  49  2.6  Discussion and conclusions We perform a large numerical parameter-space study on general aspects of structural  evolution of large strike-slip faults and related deformation fields, using a layered lithospheric model with an upper crust governed by a continuum damage rheology. One barrier to the widespread incorporation of damage rheology in crustal deformation models has been the numerous damage model parameters, whose relationships with observable phenomena are sometimes unclear. Much work in recent years has gone into relating these parameters to results of brittle deformation experiments and thermodynamic theory (e.g., Lyakhovsky et al., 2001, 2005; Hamiel et al., 2004, 2006). We take this further by using the observed shear modulus reduction in damaged fault zones to constrain the ranges of likely values for healing parameters C1 and C2. We find that admissible values of C1 and C2 are 10-24 to 10-10 s-1 and 0.015 to 0.035, respectively (Fig. 2.2). This range is significantly reduced relative to previous studies (e.g., Lyakhovsky et al., 2005). Our models with reasonable damage and viscoelastic parameters yield general deformation patterns that are comparable to those seen in natural strike-slip fault systems. Flower structures, stepovers, localized strain around fault segments and permanent damage in the shallow upper crust and within stepovers are all reproduced. Due to model simplifications (e.g., no damage accumulation in the surface layer) and element size limitations, our models cannot conclusively predict the details of surface damage patterns, the width of the fault core, or the geometry of small faults and fractures within flower structures and stepover regions. Our simulations would probably yield narrower damage zones if we incorporated depth-dependent damage-rate parameters as suggested by Lyakhovsky et al. (2005). The fault stepover zones in our models exhibit extensive damage and elasticity degradation sustained during many earthquake cycles. The simulated tensional stepovers show damage patterns consistent with intense tensile fracturing and dilation, and therefore are expected to exhibit long-lived enhanced permeability. Such damage patterns are consistent with recent structural evolution models for dilational stepovers (De Paola et al., 2007), and with mineral exploration studies that relate hydrothermal ore deposits to long-lasting extensive damage and increased permeability within fault stepovers (Micklethwaite and Cox, 2004;  50  Sheldon and Micklethwaite, 2007). The permanent damage zones our models predict should be detectable with detailed seismic and geodetic imaging studies. An important implication of the predicted damage zones at stepovers concerns the interaction between damaged material and propagating earthquakes. During the interseismic stage, weakened stepover zones experience continuous earthquakes and proportionally more inelastic strain than the surrounding crust. This reduces the interseismic stress accumulation in the stepover region, which could aid earthquake rupture arrest. Various studies address rupture propagation across stepovers using quasi-static models (e.g., Segall and Pollard, 1980), dynamic models (Harris and Day, 1999; Harris et al., 1991) and field observations (e.g., Wesnousky, 2006). Such studies show that it is easier for a rupture to jump across dilational stepovers (such as those our model produces) than compressional stepovers, because dilational normal stress brings nearby faults closer to failure. Sibson (1985) and Hamiel et al. (2005) argue, however, that during an earthquake, the normal stress change in a dilational stepover could lead to a sudden opening of fluid-filled cracks, reducing pore pressure and causing material hardening. An analogous effect is also seen in some of the models of Harris and Day (1993), where fluid within a dilational stepover inhibits the ability of rupture to propagate across it. As long-term damage accumulation and coseismic pore pressure decrease have competing effects on the elastic strength of a stepover, simulations with coupled evolution of deformation, permeability and fluid flow will be required to clarify the ramifications of damage for rupture propagation and arrest. While recent developments of our numerical framework (Hamiel et al., 2005) enable the evaluation of permeability evolution, the coupling of fluid flow and deformation is not yet incorporated in our simulations. We also note that the extensive rock damage near stepovers should produce amplified ground motion, and hence higher seismic hazard, at those regions.  While the presented model simulations have several limitations, the following features appear to be robust:   Distributed fault zone damage develops early in the evolution of a fault system (approaching steady state damage zone dimensions after .05 km of total slip in our models).    The strain generally localizes to narrow segments which are wider near the surface than at depth (representing flower structures). 51    Along fault segments, the damage heals postseismically at depths exceeding 5-10 km and is permanent at shallower depths.    The off-fault strain distribution correlates with the permanent shallow damage exhibited along fault-zones.    Stepover regions may be permanently damaged to mid-crustal depths, and such damage zones should be detectable with focused seismic and geodetic studies.  Our numerical simulations of fault zone evolution suggest that the overall aspects of fault zone deformation along large faults at seismogenic depths can be modeled effectively over time intervals of several large earthquake cycles by a collection of narrow segmented zones or planar surfaces (i.e. without incorporating off-fault damage evolution). However, the evolving fault zone structures and off-fault damage play important roles in the evolution of geometrical properties of fault sections, the formation of fault systems, and in the deformation patterns along plate boundaries. Although regional scale static stress transfer in the crust is probably not affected much by the presence of shallow weakened damage zones along fault segments, longlived volumes of extensively damaged material within fault stepovers and near other persistent geometrical irregularities can significantly affect earthquake propagation, strong ground motion, (locally) crustal stress state and surface deformation patterns.  Acknowledgments We thank Yann Klinger and Harsha Bhat for useful comments that helped to improve the readability of the paper. The study was supported by the Southern California Earthquake Center (based on NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008) and the US–Israel Binational Science Foundation, Jerusalem Israel (2004046).  52  Bibliography Agnon, A., and Lyakhovsky, V. (1995), Damage distribution and localization during dyke intrusion, In: G. Baer and A. Heimann (eds), The Physics and Chemistry of Dykes, Rotterdam, Balkema, 65-78. Ambraseys, N. N. (1970), Some characteristic features of the North Anatolian fault zone, Tectonophysics 9, 143–165. Ampuero, J.-P. and Y. Ben-Zion (2008), Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction, Geophys. J. Int., 173, 674–692, doi: 10.1111/j.1365-246X.2008.03736.x. Andrews, D. J. (2005), Rupture Dynamics with Energy Loss outside the Slip Zone, J. Geophys. Res. 110, B01307, doi:10.1029/2004JB003191. Aydin, A. and Johnson, A.M. (1983), Analysis of faulting in porous sandstones, J. Struct. Geol. 5, 19– 31. Armijo, R., Meyer, B., King, G. C. P., Rigo, A., and Papanastassiou, D. (1996), Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean, Geophys. J. Int. 126, 11–53. Baisch. S, and Bokelmann, G. H. R. (2001), Seismic waveform attributes before and after the Loma Prieta earthquake: scattering change near the earthquake and temporal recovery, J. Geophys. Res. 106, 16,323-16,337. Ben-Zion, Y. and K. Aki (1990), Seismic radiation from an SH line source in a laterally heterogeneous planar fault zone, Bull. Seism. Soc. Am., 80, 971-994. Ben-Zion, Y., T. Henyey, P. Leary and Lund, S. (1990), Observations and implications of water well and creepmeter anomalies in the Mojave segment of the San Andreas fault zone, Bull. Seismol. Soc. Am., 80, 1661-1676. Ben-Zion, Y. (1996), Stress, slip and earthquakes in models of complex single-fault systems incorporating brittle and creep deformations, J. Geophys. Res., 101, 5677-5706. Ben-Zion, Y. and Andrews, D. J. (1998), Properties and implications of dynamic rupture along a material interface, Bull. Seism. Soc. Am., 88(4), 1085–1094. Ben-Zion, Y., Z. Peng, D. Okaya, L. Seeber, J. G. Armbruster, N. Ozer, A. J. Michael, S. Baris and Aktar, M. (2003), A shallow fault zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey, Geophys. J. Int. 152, 699-717. Ben-Zion, Y. and Sammis, C. (2003), Characterization of Fault Zones, Pure and Appl. Geophys. 160, 677–715. Ben-Zion, Y. and Shi, Z. (2005), Dynamic rupture on a material interface with spontaneous generation of plastic strain in the bulk, Earth Planet. Sci. Lett., 236, 486-496, DOI: 10.1016/j.epsl.2005.03.025. Ben-Zion, Y. and Lyakhovsky, V. (2006), Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology, Geophys. J. Int. 165, 197-210. 53  Bercovici, D. and Ricard, Y. (2003), Energetics of a two-phase model of lithospheric damage, shear localization and plate-boundary formation, Geophys. J. Int. 152 (3), 581–596 doi:10.1046/j.1365-246X.2003.01854.x Budiansky, B., and O'Connell, R.J. (1976), Elastic moduli of a cracked solid, Int. J. Solids Struct. 12, 81-97. Carter, N.L. and Tsenn, M.C. (1987), Flow properties of continental lithosphere, Tectonophysics, 136, 27-63. Chester, F. M., J. P. Evans, and Biegel, R. L. (1993), Internal structure and weakening mechanisms of the San Andreas Fault, J. Geophys. Res. 98, 771–786. Chester, F. M. (1995), Geologic studies of deeply exhumed faults of the San Andreas System, Southern California: Collaborative research with Saint Louis University and Utah State University: NEHRP annual project summary, award No. 94G2457, v. 37. Chester, F. M. and Chester, J. S. (1998), Ultracataclasite structure and friction processes of the Punchbowl Fault, San Andreas System, California, Tectonophysics 295 199–221. Christensen, N. I., and Mooney, W. D. (1995), Seismic velocity structure and composition of the continental crust: A global view, J. Geophys. Res. 100 (B6), 9761-9788. Cochran, E. S., J. E. Vidale, and Li, Y. G. (2003), Near-fault anisotropy following the Hector Mine earthquake, J. Geophys. Res. 108(B9), 2436, doi:10.1029/2002JB002352. Cowie, P. A., C. Vanneste, and Sornette, D. (1993), Statistical Physics Model for the Spatiotemporal Evolution of Faults, J. Geophys. Res. 98(B12), 21,809–21,821. Cundall, P. A. and Board, M. (1988), A microcomputer program for modeling large-strain plasticity problems, in: Numerical methods in geomechanics, Proc. 6th Int. Conf. Numerical Methods in Geomechanics, Innsbruck, ed. Swoboda, C., Rotterdam, Balkema, pp. 2101-2108. De Paola, N., Holdsworth, R. E., Collettini, C., McCaffrey, K. J. W. and Barchi, M. R. (2007), The structural evolution of dilational step-overs in regional transtensional zones, in Cunningham W. D. and Mann, P. (eds), Tectonics of Strike-Slip Restraining and Releasing Bends. Geological Society, London. Special Publication, 290, pp. 433445. Dieterich, J. H. (1972), Time-dependent friction in rocks, J. Geophys. Res. 77, 3690-3697. Dieterich, J. H. (1978), Time-dependent friction and the mechanics of stick-slip, Pure and Appl. Geophys. 116, 790-805. Dieterich, J. H. (1979), Modeling of rock friction 1. Experimental results and constitutive equations, J. Geophys. Res. 84, 2161-2168. Dolan, J. F., D. D. Bowman, and C. G. Sammis. (2007), Long-range and long-term fault interactions in Southern California, Geology, 35, 855-858. Dor, O., Rockwell, T. K., and Ben -Zion, Y. (2006), Geological observations of damage asymmetry in the structure of the San Jacinto, San Andreas and Punchbowl Faults in Southern California: a possible indicator for preferred rupture propagation direction, Pure appl. Geophys. 163, 301–349, DOI 10.1007/s00024-005-0023-9. 54  Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., Duman, T. Y. (2008), Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction, Geophys. J. Int., doi: 10.1111/j.1365-246X.2008.03709.x. Dunham, E. M. and J. R. Rice (2008), Earthquake slip between dissimilar poroelastic materials, J. Geophys. Res., in press. Eissa, E. A. and Kazi, A. (1988), Relation between static and dynamic Young‟s Moduli of Rocks, Int. J. Rock Mech. Min Sci and Geomech. Abstr., 25 (6), 479-482. Evans, J. P., Shipton, Z. K., Pachell, M. A., Lim, S. J., and Robeson, K. (2000), The Structure and Composition of Exhumed Faults, and their Implication for Seismic Processes, In Proc. of the 3rd Confer. on Tecto. problems of the San Andreas system, Stanford University. Fialko, Y., Sandwell, D., Agnew, D., Simons, M., Shearer, P., and Minster, B. (2002), Deformation on nearby faults induced by the 1999 Hector Mine earthquake, Science, 297, 1858-1862. Fialko, Y. (2004), Probing the mechanical properties of seismically active crust with space geodesy: Study of the co-seismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. Geophys Res. 109, B03307, doi:10.1029/2003JB002756. Finzi, Y., Hearn, E. H., Lyakhovsky, V. and Ben-Zion. Y. (2006), 3D Viscoelastic damage rheology models of strike-slip fault systems and their associated surface deformation, Eos Trans. AGU, 87(52), Fall Meet. Suppl., Abstract T21C-0425. Graymer, R.W., Langenheim, V.E., Simpson, R.W., Jachens, R.C., and Ponce, D.A. (2007), Relatively simple through-going fault planes at large-earthquake depth may be concealed by the surface complexity of strike-slip faults, in Cunningham, W.D., and Mann, Paul, eds., Tectonics of Strike-Slip Restraining and Releasing Bends: Geological Society of London Special Publication, v. 290, p. 189-201. doi: 10.1144/SO290.5 0305-8719/07. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y. and Lockner, D. (2004), A Visco-elastic damage model with applications to stable and unstable fracturing, Geophys. J. Int. 159, 1155–1165. Hamiel, Y., Lyakhovsky, V. and Agnon, A. (2005), Rock dilation, nonlinear deformation, and pore pressure change under shear, Earth Planet. Sci. Lett. 237, 577-589. Hamiel, Y., Katz, O., Lyakhovsky., Reches, Z. and Fialko, Y. (2006), Stable and unstable damage evolution in rocks with implications to fracturing of granite, Geophys. J. Int. 167, 1005-1016. Hamiel, Y., and Fialko, Y. (2007), Structure and mechanical properties of faults in the North Anatolian Fault system from InSAR observations of coseismic deformation due to the 1999 Izmit (Turkey) earthquake, J. Geophys Res. 112,B07412, doi:10.1029/2006JB004777. Harris, R. A., Archuleta, R. J., and Day, S. M. (1991), Fault steps and the dynamic rupture process: 2-D numerical simulations of a spontaneously propagating shear fracture, Geophys. Res. Lett. 18, 893-896. 55  Harris, R. A., and Day, S. M. (1993), Dynamics of Fault Interaction: Parallel Strike-slip Faults, J. Geophys. Res. 98, 4461–4472. Harris, R. A., and Day, S. M. (1999), Dynamic 3D simulations of earthquakes on en echelon faults, Geophys. Res. Lett. 26, 2089-2092. Hickman, S., Sibson, R.H. and Bruhn, R. (1995), Introduction to a special section, mechanical involvement of fluids in faulting. J. Geophys. Res. 100, 12,831– 12,840. Ide, J. M. (1936), Comparison of statically and dynamically determined Young‟s modulus of rocks, Proc. Nat. Acad. Sci., USA, 22, 81-92. Karabulut, H. and Bouchon, M. (2007), Spatial variability and non-linearity of strong ground motion near a fault, Geophys. J. Int. 170, 1, 262-274. Kim, Y. S., Peacock, D. C. P., and Sanderson, D. J. (2004), Fault damage zones, J. of Struct. Geology, v. 26, p. 503-517. King, G. (1986), Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure, Pure Appl. Geophys. 124, 567–585. Kirby, S. H. and Kronenberg, A. K. (1987), Rheology of the lithosphere: selected topics, Rev. Geophys. 25: 1219-1244. Korneev, V.A., Nadeau, R.M. & McEvilly, T.V. (2003). Seismological Studies at Parkfield IX: Fault-one imaging using Guided Wave Attenuation, Bull. Seism. Soc. Am., 93, 14151426. Lewis, M. A., Peng, Z., Ben-Zion, Y., and Vernon, F. (2005), Shallow seismic trapping structure in the San Jacinto fault zone, Geophys. J. Int. 162, 867–881, doi:10.1111/j.1365246X.2005.02684.x, 2005. Li, Y. G., P. Leary, K. Aki, and P. Malin (1990), Seismic trapped modes in the Oroville and San Andreas fault zones, Science, 249, 763-766. Li, Y.-G., Aki, K., Adams, D., Hasemi, A., and Lee, W. H. K. (1994), Seismic guided waves trapped in the fault zone of the Landers, California, earthquake of 1992, J. Geophys. Res. 99(B6), 11,705–11,722. Li, Y.-G., Chen, P., Cochran, E. S., Vidale, J. E., and Burdette, T. (2006), Seismic evidence for rock damage and healing on the San Andreas fault associated with the 2004 M 6.0 Parkfield Earthquake, Bull. Seism. Soc. Am. 96, no. 4B, S349–S363. Liu, Y., Teng T. L., and Ben-Zion, Y. (2004), Systematic analysis of shear wave splitting in the aftershock region of the 1999 Chi-Chi earthquake: Evidence for shallow anisotropic structure and lack of systematic temporal variations, Bull. Seismol. Soc. Am., 94, 23302347. Lyakhovsky, V. and Myasnikov, V. P. (1984), On the behavior of elastic cracked solid, Phys. Solid Earth, 10, 71-75. Lyakhovsky, V. and Myasnikov, V. P. (1985), On the behavior of visco-elastic cracked solid, Phys. Solid Earth, 4, 28-35. 56  Lyakhovsky, V., Ben-Zion, Y. and Agnon, A. (1997a) Distributed damage, faulting, and friction, J. Geophys. Res. 102, 27,635-27,649. Lyakhovsky, V., Reches, Z., Weinberger, R. and Scott, T.E. (1997b), Non-linear elastic behavior of damaged rocks, Geophys. J. Int. 130, 157-166. Lyakhovsky, V. Ben-Zion, Y. and Agnon, A. (2001), Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere, J. Geophys. Res. 106: 4103-4120. Lyakhovsky, V., Ben-Zion, Y. and Agnon, A. (2005), A viscoelastic damage rheology and rateand state-dependent friction, Geophys. J. Int. 161, 179-190. Lyakhovsky, V., and Ben-Zion, Y. (2008), Scaling relations of earthquakes and aseismic deformation in a damage rheology model, Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365-246X.2007.03652.x. Lyakhovsky V., and Ben-Zion, Y. (2009), Evolving fault zone structures in a damage rheology model, in review, Geochem., Geophys., Geosys. Malvern, L.E. (1969), Introduction to the Mechanics of a Continuum Medium, New Jersey, Prentice-Hall, Inc., 713 p. McGuire, J. and Ben-Zion, Y. (2005), High-resolution imaging of the Bear Valley section of the San Andreas Fault at seismogenic depths with fault-zone head waves and relocated seismicity, Geophys. J. Int. 163, 152-164, doi: 10.1111/j.1365-246X.2005.02703.x. Micklethwaite, S. and Cox, S. F. (2004), Fault-segment rupture, aftershock-zone fluid flow and mineralization, Geology, 32, 813–816. Mooney, W. D., and Ginzburg, A. (1986), Seismic measurements of the internal properties of fault zones, Pure Appl. Geophys. 124, 141–157. Nadeau R. Antolik M. Johnson P. Foxall W. and McEvilly T. V. (1994), Seismological studies at Parkfield III: Microearthquake clusters in the study of fault-zone dynamics, Bull. Seism. Soc. Am. 83, 247-263. Oglesby, D. D., Day, S. M., Li, Y.-G., and Vidale, J. E. (2003), The 1999 Hector Mine Earthquake: the dynamics of a branched fault system, Bull. Seismol. Soc. Am., 93, 6, 2459–2476. Oglesby, D. D., 2005. The dynamics of strike-slip step-overs with linking dip-slip faults. Bull. Seismol. Soc. Am. 95, 1604–1622. Olson J., and Pollard D. D., (1989), Inferring paleostresses from natural fracture patterns: A new method, Geology: Vol. 17, No. 4 pp. 345–348. Peng, Z., Ben-Zion, Y., Michael A. J., and Zhu, L. (2003), Quantitative analysis of seismic trapped waves in the rupture zone of the 1992 Landers, California earthquake: Evidence for a shallow trapping structure, Geophys. J. Int. 155, 1021-1041. Peng, Z. and Ben-Zion, Y. (2004), Systematic analysis of crustal anisotropy along the KaradereDuzce branch of the north Anatolian fault, Geophys. J. Int. 159, 253-274, doi:10.1111/j.1365- 46X.2004.02379.x.  57  Peng, Z. and Ben-Zion, Y. (2005), Spatio-temporal variations of crustal anisotropy from similar events in aftershocks of the 1999 M7.4 İzmit and M7.1 Düzce, Turkey, earthquake sequences, Geophys. J. Int., 160(3), 1027-1043, doi: 10.1111/j.1365-246X.2005.02569.x. Peng, Z. and Ben-Zion, Y. (2006), Temporal changes of shallow seismic velocity around the Karadere-Duzce branch of the north Anatolian fault and strong ground motion, Pure Appl. Geophys. 163, 567-600, DOI 10.1007/s00024-005-0034-6. Poliakov, A., Cundall, P., Podladchikov, Y. and Lyakhovsky, V. (1993), An explicit inertial method for the simulation of viscoelastic flow: an evaluation of elastic effects on diapiric flow in two- and three-layers model, In Proc. NATO Advanced Study Institute on Dynamic Modeling and Flow in the Earth and Planets, Runcorn, K.E. and Stone, D. (eds) Dordrecht, Kluwer, pp. 175-195. Poupinet, G., Ellsworth, W. L., and Frechet, J. (1984), Monitoring Velocity Variations in the Crust Using Earthquake Doublets: An Application to the Calaveras Fault, California, J. Geophys. Res. 89(B7), 5719–5731. Powell, R. E., and Weldon, R. J. (1992), Evolution of the San Andreas Fault, Annu. Rev. Earth Planet. Sci. 20, 431-468. Revenaugh, J. (2000), The relation of crustal scattering to seismicity in southern California, J. Geophys. Res. 105(B11), 25,403–25,422. Rice, J.R. and Ben-Zion, Y. (1996), Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci. USA 93: 3811-1818. Rockwell, T. K. and Ben-Zion Y. (2007), High localization of primary slip zones in large earthquakes from paleoseismic trenches: Observations and implications for earthquake physics, J. Geophys. Res. 112, B10304, doi:10.1029/2006JB004764. Rubinstein, J.L. and Beroza, G.C. (2004), Evidence for widespread strong ground motion in the Mw 6.9 Loma Prieta earthquake, Bull. Seism. Soc. Am., 94, 1595–1608. Rudnicki, J. W., and Rice, J. R. (2006), Effective normal stress alteration due to pore pressure changes induced by dynamic slip propagation on a plane between dissimilar materials, J. Geophys. Res. 111, B10308, doi:10.1029/2006JB004396. Saucier, F., Humphreys, E.D., and Weldon, R.J., (1992), Stress near geometrically complex strike-slip faults: application to the San Andreas Fault at Cajon Pass, southern California, J. Geophys. Res., 97, 5081–5094. Schaff, D. P., Bokelmann, G. H. R., Beroza, G. C., Waldhauser, F., and Ellsworth, W. L. (2002), High-resolution image of Calaveras Fault seismicity, J. Geophys. Res. 107(B9), 2186, doi:10.1029/2001JB000633. Schaff, D. P., and Beroza, G. C. (2004), Coseismic and postseismic velocity changes measured by repeating earthquakes, J. Geophys. Res. 109, B10302, doi:10.1029/2004JB003011.  58  Schulz, S. E. and Evans, J. P. (2000), Mesoscopic Structure of the Punchbowl Fault, Southern California and the Geologic and Geophysical Structure of Active Strike-slip Faults, J. of Struct. Geol. 22, 913–930. Segall, P., and Pollard, D. D. (1980), Mechanics of discontinuous faults, J. Geophys. Res. 85, 4337-4350, 1980. Sengor, A. M. C., Tuysz, O., Imren, C., Sakınc, M., Eyidogan, H., Gorur, N., Le Pichon, X. and Rangin, C. (2005), The North Anatolian Fault: A New Look, Annu. Rev. Earth Planet. Sci. 33: 37–112. Sheldon H. A., and Micklethwaite S. (2007), Damage and permeability around faults: Implications for mineralization, Geology: Vol. 35, No. 10 pp. 903–906. Shipton Z. K., and Cowie P. A. (2003), A conceptual model for the origin of fault damage zone structures in high-porosity sandstone, J. Struct. Geol., 25, 3, 333-344 Sibson, R.H. (1985), Stopping of earthquake ruptures at dilational fault jogs, Nature 316, 248– 251. Sibson, R.H. (2003), Thickness of the Seismic Slip Zone, Bull. Seismol. Soc. Am. 93, 3, 11691178. Stirling, M.W., Wesnousky, S.G. and Shimazaki, K. (1996), Fault trace complexity, cumulative slip, and the shape of the magnitude-frequency distribution for strikeslip faults: a global survey, Geophys. J. Int. 124, 833–868. Stierman, D. J. (1984), Geophysical and geological evidence for fracturing, water circulation and chemical alteration in granitic rocks adjacant to major strike-slip faults, J. Geophys. Res. 89, B7 5849-5857. Sylvester, A. G., and Smith, R, (1976), Tectonic transpretions and basement controlled deformation in the San Andreas fault zone, Salton trough, California, AAPG Bull. 60, 2081-2102. Sylvester, A. G. (1988), Strike-slip faults, Geol. Soc. Am. 100, 1666-1703. Tchalenko, J.S. (1970), Similarities between shear zones of different magnitudes, Geological Society of America Bulletin 81, 1625–1640. Templeton, E. L., and Rice, J. R. (2008), Off-fault plasticity and earthquake rupture dynamics, 1. Dry materials or neglect of fluid pressure changes, J. Geophys. Res. Tenthorey, E., Cox, S. F. and Todd, H. F. (2003), Evolution of strength recovery and permeability during fluid-rock reaction in experimental fault zones, Earth Planet. Sci. Lett., 206(1–2), 161–172. Thurber, C., Zhang, H., Waldhauser, F., Hardebeck, J., Michael, A., and Eberhart-Phillips, D. (2006), Three-dimensional compressional wavespeed model, earthquake relocations, and focal mechanisms for the Parkfield, California, Region. Bull. Seis. Soc. of Am., 96, 4B, S38–S49, doi: 10.1785/0120050825 Turcotte, D. L., and Glasscoe, M.T. (2004), A damage model for the continuum rheology of the upper continental crust, Tectonophysics, 383, 71-80. 59  Wesnousky, S. (1994), The Gutenberg-Richter or characteristic earthquake distribution, which is it?, Bull. Seismol. Soc. Am. 84, 1940–1959. Wesnousky, S. G. (2006), Predicting the endpoints of earthquake ruptures, Nature, 444, 358-360. Wilcox, R. E., Harding, T. P., and Seely, D. R. (1973), Basic wrench tectonics, AAPG Bull., 57, 74-96. Wu, C. Z. Peng and Y. Ben-Zion (2009), Non-linearity and Temporal Changes of Fault Zone Site Response Associated with Strong Ground Motion, Geophys. J. Int., 176, 265-278, doi: 10.1111/j.1365-246x.2008.04005.x. Yang, W. and Y. Ben-Zion, (2009), Observational analysis of correlations between aftershock productivities and regional conditions in the context of a damage rheology model, Geophys. J. Int. (2009) 177, 481–490 doi: 10.1111/j.1365-246X.2009.04145.x Yamashita, T. (2007), Postseismic quasi-static fault slip due to pore pressure change on a bimaterial interface, J. Geophys. Res. 112, B05304, doi:10.1029/2006JB004667.  60  3  DAMAGE-ZONE HEALING AND STRIKE SLIP FAULT EVOLUTION: NUMERICAL RESULTS AND GEOPHYSICAL CONSTRAINTS2  Summary Numerical simulations of long-term crustal deformation reveal the important role that damage healing (i.e. fault-zone strengthening) plays in the structural evolution of strike-slip fault systems. Geophysical observations indicate that after an early post-seismic stage of rapid healing, rock damage in fault zones approaches asymptotically values which may vary markedly with location. While variable inter-seismic damage levels are expected for a variety of reasons (e.g., variations in rock type, local stress patterns, heat flow and fluid content), the high uncertainties in estimates of damage levels and fault zone properties at seismogenic depth make it difficult to correlate between fault zone conditions and healing effectiveness. Nevertheless, the observations motivate exploring the sensitivity of fault zone structure and evolution to reasonable variations in healing effectiveness. Here we perform 3D simulations of fault system evolution applying a continuum damage rheology and we find that the spatial extent of damage zones and the long-term geometrical complexity of strike-slip fault systems are strong functions of the healing parameters. Specifically, simulations with highly-effective healing form localized intensively damaged fault cores that interseismically extend only a few kilometers deep, and are bracketed by wide zones of distributed off-fault damage. Ineffective healing yields deeper zones of intense damage that persist throughout the interseismic interval, and narrower zones of distributed off-fault damage. Furthermore, highly-effective healing leads to a rapid evolution of an initially segmented fault system to a simpler through-going fault, while ineffective healing along a segmented fault preserves complexities such as stepovers and fault jogs and results in long-lasting, distributed deformation.  2  A version of this chapter has been submitted for publication. Finzi, Y., Hearn, E.H., Ben-Zion, Y. and Lyakhovsky V. Damage-zone healing and strike slip fault evolution: Numerical results and geophysical constraints.  61  3.1  Introduction In a typical fault zone, the principal slip surface is often surrounded by breccia and  embedded within a damage zone that consists of open and healed fractures, veins, and other secondary features. Nested zones with increasing damage are present: Distributed damage may be found up to several kilometers from the fault core near the Earth‟s surface around large strikeslip faults. Within about 1 km of the active fault core, a zone of more intense damage is present. This zone may produce local seismic anisotropy (e.g., Boness and Zoback, 2004), concentrate coseismic strain (e.g. Fialko et al. 2002), and elevate seismic scattering (e.g., Revenaugh, 2000; Peng and Ben-Zion, 2006). Within this zone, there is a 10-100 m-wide zone of intense damage which can trap seismic waves (e.g., Li et al., 1994; Ben-Zion et al., 2003) and may appear at the surface as a belt of pulverized fault zone rocks (e.g., Dor et al., 2006, 2008). Recent seismic and geodetic observations at several fault zones show variable levels of interseismic damage, suggesting that healing effectiveness varies within fault systems (e.g., Hearn and Fialko, 2009; Cochran et al., 2009). Damage intensity along strike-slip faults is typically higher near stepovers, fault kinks and fault bends, which exhibit strain hardening associated with secondary fractures at different scales (e.g., Kim et al., 2004). Persisting rigidity loss and modification of the stress state at these localities can affect rupture propagation, ground shaking, and fluid flow patterns (Biegel et al., 2008; Bhat et al., 2007, Micklethwaite and Cox, 2006). However, due to difficulties in simulating fault growth in 3D, the long-term structural evolution of fault zone damage has not been thoroughly examined so far. Seismic and structural studies indicate that the complexity (segmentation and geometric disorder) of a fault system reduces with time as the fault accumulates slip (e.g., Wesnousky, 1994; Ben-Zion and Sammis, 2003). Ben-Zion et al. (1999) and Lyakhovsky et al. (2001) suggest that the ratio of the timescale for material healing to the timescale of stress loading governs the evolution of fault zone structure and seismicity patterns, with high ratios (“long memory”) leading to simple localized fault systems and characteristic earthquake statistics and low ratios (“short memory”) producing persisting structural complexities and power law earthquake statistics. Laboratory and field studies indicate that rock healing is logarithmic in time, with significant rapid healing on a timescale of a day or so and slow additional small healing on longer times (e.g., Dieterich and Kilgore, 1996; Johnson and Jia, 2005; Wu et al., 62  2009). Since the loading rates of faults vary typically only by one order of magnitude or less, the ratio of timescales for healing over loading does not provide sufficient information to distinguish between observed variations among different faults.  In the current paper we analyze how material healing effects damage zone structure and fault system evolution. Our simulations show that the damage zone structure along faults, the temporal stability of geometrical complexities, and the evolution of fault configuration and strain distribution patterns are all strongly affected by the effectiveness of healing. Specifically, more effective healing yields wider but shallower damage zones and shorter lived stepovers (and other fault complexities) compared to damage zones formed in conditions less favorable for healing. In the following section we briefly describe the damage rheology model and discuss the seismic and geodetic constraints on the damage healing rate parameters. In section 3.3 we present the effects of healing effectiveness on simulated damage zone structure, fault complexity, and fault system evolution.  3.2  Damage rheology and healing formulation Many current, theoretical fault system studies employ numerical simulations in which  brittle rock is modeled as a rigid elastic-plastic solid, and faults are prescribed surfaces governed by friction. Such frictional models typically assume that deformation at all stages occurs on predefined surfaces, and they do not provide a mechanism for understanding distributed deformation and the evolution of fault zones. Several recent studies that apply frictional surfaces simulate inelastic off-fault yielding as plastic strain (e.g., Andrews, 2005; Ben-Zion and Shi, 2005; Ma, 2008). However, inelastic deformation in the brittle portion of the lithosphere is associated with distributed cracking that modifies its elastic properties (e.g., Lockner et al., 1977; Jaeger and Cook, 1979). Simulating the evolution of fault zones and related deformation fields requires concepts such as damage rheology that accounts for evolving elastic properties, off-fault deformation and spontaneous formation of new fault surfaces.  63  3.2.1 Damage rheology: Theory To integrate fault zone weakening and healing in numerical simulations of large-scale fault zone evolution, we apply concepts of continuum damage mechanics. Based on the crack density formulation of Budiansky & O'Connell (1976), the damage model we use (Lyakhovsky et al., 1997a, 1997b) relates the elastic moduli to a single variable. This scalar damage variable () reflects the extent of strength degradation due to crack formation and opening, with = 0 indicating undamaged rock and = 1 corresponding to a complete strength loss. Our damage rheology model is briefly described below. For further details and selected applications of this model we refer to previous publications (Lyakhovsky et al., 1997a, b, 2005; Hamiel et al., 2004; Ben-Zion and Lyakhovsky, 2006). To account for the evolution of material properties, the damage rheology model introduces a third elastic modulus () and makes the elastic rigidity () a function of the damage state variable, as follows:    = constant;     ·m0–) ≈ 0    m      where  and  are the Lamé parameters of linear Hookean elasticity, and  = I1/√I2 is the strain invariants ratio (I1=kk and I2=ijij are the first and second invariants of the elastic strain tensor  ij). The parameter 0 separates states of material degradation ( > ) and healing ( < ) associated with positive and negative damage evolution, respectively. As damage accumulates (i.e.  increases), the rigidity (shear modulus  decreases and the modulus  increases from 0 (damage free) to m, amplifying the non-linearity of rock elasticity.  In this study we focus on damage healing. At high confining pressure, low shear stress and high temperature, healing of micro-cracks is favored, resulting in a recovery of elastic moduli. Motivated by observations of a logarithmic increase of the friction coefficient with time (e.g., Dieterich 1978), Lyakhovsky et al. (1997a) suggested the following healing function:  64    d  C1  exp  I 2    0  dt  C2   for    0  (3.2)  According to equation (3.2) fault zone healing rate is a function of the material properties (C1, C2 and 0), the state of the fault zone (and I2), and the ratio between confinement and shear (). For a given strain state and level of damage, the rate and the overall effectiveness of healing are primarily determined by the healing rate (material) parameters C1 and C2. Depending on the combination of these parameters, the healing process may be fast or slow, and may yield insignificant to complete healing of the damaged material (Fig. 3.1). Finzi et al. (2009) express the healing parameter space in terms of two descriptive parameters: the healing time scale (h) and the characteristic interseismic damage level (ch). We define the time scale for healing as the time at which the relative change of elastic moduli falls below an arbitrarily chosen rate of 0.1% yr-1 (d/dt=3e-11sec-1). We consider this rate a suitable threshold with slower healing corresponding to a steady state damage level and insignificant healing (i.e. less than 10% healing expected during time intervals of the order of 100 years). The characteristic interseismic damage level ch (termed the “final” damage level f in Finzi et al., 2009), is a first-order predictor of the damage level along a simple fault segment at time t= h. To estimate ch we assume a strain invariant ratio suitable for healing (0 ≈ -1), set the healing rate in equation 3.2 to the chosen threshold, and derive the damage level expected as the healing becomes slower than 0.001yr-1:          ch = C2 ln [3e-11/ C1I2]  (3.3)  It is important to note that the interseismic damage level (both in our numerical simulation and along natural fault zones), is expected to depart from ch values as a result of local variations in stress, in loading history and in fault geometry. Nevertheless, a material with healing parameters corresponding to high ch (e.g., ch > 0.6) is expected to undergo ineffective healing, whereas a material with parameters C1 and C2 corresponding to very low ch (e.g., ch < 0.3) is expected to undergo highly effective healing. The parameter ch also depends on the 65  applied strain and the depth at which the healing occurs (through I2 in Eq. 3.3). In this paper all  ch values are calculated assuming a depth of 3-5 km and a strain invariant ratio suitable for healing. To derive the time scale for healing, we integrate equation 3.2 (assuming constant strain) and get an analytic expression of t which is rewritten to express the healing time required for the damage level to decrease from =1 to =ch:   1   ch   1 exp C 2   h  C  1  I 2 1 exp  C2  C2   (3.4)  The time scale for healing (h) indicates how fast a damaged material heals to a quasisteady state of very slow healing, however it does not reveal the effectiveness of healing or the characteristic damage level (Fig. 3.1). The healing time scales associated with parameters admissible for simulations of natural fault zones range within one order of magnitude, between  h=10 yr and h=102 yr.  As evident from Figure 3.1, to fully determine the healing parameters of a material one would need either (1) the damage level at two different times after failure, or (2) the current damage level, the time elapsed since failure, and the time scale relevant to the specific healing process. Until recently, damage-healing data from fault zones was unavailable and model parameters could only be constrained based on analytical considerations and laboratory experiments. In such an effort, Lyakhovsky et al. (2005) obtained an important relation between the two healing parameters based on empirical rate and state friction results:     2 C1  BC2 exp  0  /  cmp  C2   (3.5)  where B (~1-2 s-1) is a laboratory-determined time scale for the evolution of the friction coefficient with hold time (Dieterich, 1978), cmp is the compaction strain estimated by the ratio between lithostatic stress and the bulk modulus (K), and 0 is the damage level at t=0 ( healing 66  starts immediately after failure with  ≈1). Finzi et al. (2009) used this relation between C1 and C2 to narrow the range of admissible healing parameters for fault evolution models, showing that natural healing processes may be represented by C2 values between 0.015 and 0.035, and C1 values between C1~10-24 s-1 and 10-10 s-1. As detailed in the following subsection, new geophysical data support this conclusion.  Figure 3.1. Damage as a function of time (after failure) for various sets of parameters C1 and C2. The healing processes plotted in red have short healing time scales ( h ~ 16 yr) and therefore exhibit fast initial healing and slow insignificant healing at t >h. The healing processes with long time scales (h ~ 50 yr; blue curves) exhibit slower initial healing but may undergo some possibly significant healing long after t=h. The effectiveness of healing (i.e. the expected typical interseismic damage level) depends primarily on ch and is independent of h. The healing time scales and characteristic damage levels are indicated with diamond symbols on the healing curves.  67  3.2.2 Healing of damaged material: Constraints from fault zone observations Due to a shortage of observations of damage levels in individual fault zones, for a range of different times in the earthquake cycle, previous attempts to constrain healing parameters have been based on laboratory experiments and analytic considerations (e.g., Ben-Zion et al., 1999; Lyakhovsky et al., 2001). Laboratory stick-slip experiments on a wide range of materials show logarithmic fault strength recovery with time (e.g., Dieterich, 1978; Dieterich and Kilgore, 1996). Logarithmic healing is also compatible with various seismic observations of fast active fault-zone healing over time scales of days to months, followed by very little additional healing during the interseismic stage of the earthquake cycle (e.g., Karabulut and Bouchon, 2007; Wu et al., 2009). Based on estimates of fault zone damage levels derived from a limited set of trapped seismic wave and geodetic observations (Peng et al., 2003; Lewis et al., 2005; Fialko et al., 2002; Li et al., 2006; Hamiel and Fialko, 2007) and assuming logarithmic healing, Finzi et al. (2009) constrained the healing time scale (h<40 yrs) and the characteristic interseismic damage level (ch > 0.4). The characteristic interseismic damage was also set a maximum value of  ch = 0.75, based on numerical indications that shallow material (z < 3 km) could not sustain damage levels higher than  = 0.75 for long time periods (Finzi et al., 2009).  In our current simulations, the C1 and C2 values correspond to healing time scales similar to those applied in previous studies (within the range represented in Figures 3.1 and 3.2), but the range of ch values is broader than that previously applied. While values of ch were previously assumed to range between 0.4 and 0.5, we currently select healing parameters that represent ch values between 0.25 and 0.65 based on estimates of fault zone rigidity along the ECSZ (Hearn and Fialko 2009; Cochran et al., 2009; Barbot et al., 2009). This wider range of admissible ch values is a more conservative assumption, given the uncertainties in estimates of interseismic damage levels and the likely variability in healing effectiveness (due to local variation in stress, strain history, fluid flow and rock type). Our appreciation of the natural variability in fault zone healing effectiveness and the wider range of parameters in our current simulations (representing the full range of healing effectiveness) enable us to identify the previously underestimated effect that healing has on damage zone structure and evolution. Our current results (Section 3.3)  68  improve and correct our previous analysis (Chapter 2, Fig. 2.8). The healing parameter space suggested for simulations of natural fault zone evolution is outlined in Figure 3.2.  Figure 3.2. Geophysical, analytical and laboratory-based constraints on healing parameters. Minimum and maximum expected interseismic damage (min= 0.25, max= 0.65) and maximum healing time scale (h=40 yr) superimposed on the parameter space resulting from the analytic considerations (Eq. 3.5). The space between the minimum and maximum damage levels and above the 40 yr healing time scale represents natural fault zone healing processes. While this parameter space is comparable to that previously suggested (Finzi et al., 2009), in the current work we apply a wider range of parameters to better represent the variations in expected interseismic damage levels within this parameter space.  69  3.3  Results  3.3.1 Fault zone structure as a function of healing To investigate the structure of damage zones associated with evolving strike slip fault systems, we use three dimensional realizations of a strike-slip fault within a layered seismogenic crust governed by damage rheology, underlain by viscoelastic layers representing the lower crust and upper mantle. We use diabase and dunite flow laws and assume a fixed geotherm (20°C/km) in our simulations. The modeled region is 100 to 250 km long (along strike), 100 km wide and 50 km deep. More details of our typical model setups and parameters are given by Finzi et al. (2009). The top 3 km of the model represents a sedimentary layer weaker (lower rigidity) than the underlying crystalline rock. As frequent failure of surface elements due to low confining stress would be computationally time consuming, we suppress damage accumulation in this layer. Various test models with damage not suppressed within the surface layer (not presented here) suggest that the surface damage zone associated with a strike-slip fault in such a simulation may be slightly wider, and consist of higher inter-seismic damage levels. A variable force boundary condition (Lyakhovsky and Ben Zion, 2008) is applied to the sides and bottom of the model domain, simulating a constant far field fault parallel velocity with relative rate of 32 mm/yr (corresponding to the San Andreas Fault). The boundary driving forces are set to represent imposed fault zones outside the model domain, and they induce faulting near the centers of the north and south edges of the model (Lyakhovsky and Ben-Zion, 2008; Finzi et al., 2009).  Our simulations output the distribution of strain and damage within the model domain, from which we can calculate surface velocities, rigidity and other related quantities. Figure 3.3 illustrates some features which correspond to observed geological structures such as fault segments, stepovers, and flower structures. Contiguous sets of elements that fail repeatedly, resulting in a higher damage level and lower rigidity than their surroundings, are interpreted as fault segments (Fig. 3.3). Because of their relative weakness, these fault segments are also the centers of high strain rates. Cross-sections through modeled fault segments (Fig. 3.3) display “flower structures” with depth, which comprise of localized damage along the active fault core and a broader zone of distributed damage in the top 3-10 kilometers of the crust (Fig. 3.3). The 70  overall geometry of the “flower-like” damage zones in our simulations is compatible with field observations (e.g., Sylvester, 1988; Rockwell and Ben-Zion, 2007) and with numerical simulations of dynamic ruptures with off-fault yielding in the form of plasticity (e.g., Ma, 2008). Based on these observations, we define two damage sub-zones that are distinct in their evolution patterns, damage level, and spatial distribution. Localized Active Fault (LAF) damage represents the highly localized damage along the active fault cores (Fig. 3.3). LAF damage is coseismically very high along the primary slip zone, but it rapidly heals at depth. Distributed Off Fault (DOF) damage is the sustained, cumulative damage resulting from many earthquakes. The DOF damage develops during the early stages of fault system evolution, and thereafter its spatial extent is stable and the degree of damage within it evolves (Fig. 3.3). The width of the DOF damage is 6-13 km in our simulations, and its average depth extent (away from the active fault core) is up to 7 km (Fig. 3.3). The width of the LAF damage represents the maximum width of the localized damage along the fault core (typically one element wide, indicating plausible localization to slip surfaces or zones narrower than our numerical resolution), and its depth is the maximum depth extent of the damage zone (Fig. 3.3). Descriptive analyses of damage structures along simulated strike-slip fault segments and along fault stepover zones are given in Finzi et al. (2009).  71  Figure 3.3. Distributed off fault (DOF) and localized active fault (LAF) damage along a simulated strike-slip fault. The vertical cross-section (left) shows a simulated flower structure with localized damage at depth (along the active fault) and distributed damage near the surface which is the cumulative effect of past earthquakes. The horizontal sections (right) show persistent distributed damage along fault traces (typically in the top 3-5 km) and at considerable depth within fault stepovers.  While our previous analysis concluded that the flower structure is a robust feature that shows limited sensitivity to variations in the healing parameters C1 and C2 (Finzi et al., 2009, Section 2), we now realize that the effect of healing on damage zone structure is greatly enhanced when co-varying both parameters to represent a wide range of ch values. By thoroughly exploring a broader, geophysically constrained healing parameter space, we show here that the properties of simulated fault zone structures and related deformation field are significantly affected by the effectiveness of healing. The direct result of healing in our simulations is postseismic restrengthening of the faultzone and therefore a reduction of the damage zone extent. The deeper section of the fault zone (the deep LAF damage) is most susceptible to healing due to the relatively high confining stress, and is therefore most sensitive to variations in the healing parameters. Our results indicate that 72  the depth extent of the LAF damage (after the early postseismic interval) ranges from the entire seismogenic crust in fault zones experiencing ineffective healing ( ch > 0.65; Fig. 3.4 bottom panel; Fig. 3.5) to a few kilometers in fault zones experiencing effective healing ( ch < 0.4; Fig. 3.4 top panel; Fig. 3.5). During most of the seismic cycle, intense and contiguous LAF damage is typically limited to the top section of the crust (e.g., Fig. 3.3, Fig. 3.4 top panel; Fig. 3.5). The dimensions of the shallow DOF damage are also sensitive to variations in the healing parameters, with ineffective healing resulting in faster (more complete) strain localization and therefore narrower fault zones and more effective healing yielding wider fault zones.  It should be noted that our simulations probably overestimate the maximum width of fault damage zones as they do not incorporate depth-dependent damage-rate parameters which would yield less damage accumulation in the top few kilometers of the crust (Lyakhovsky et al., 2005). In addition, our model element dimensions (1-3 km wide) may result in overestimated damage zone width, and they do not enable characterization of small scale structures and extreme strain localizations expected along active slip zones (e.g., Chester et al., 1993; Rockwell and Ben-Zion, 2007).  73  Figure 3.4. Cross-sections displaying evolving damage levels along two strike-slip faults representing effective (top, ch=0.4) and ineffective (bottom, ch=0.7) healing conditions. All cross-sections are 25 km deep and 16 km wide.  74  Figure 3.5. Sensitivity of damage zone dimensions to healing effectiveness (ch). The width of the DOF damage decreases in simulations with ineffective healing (higherch), and the depth of the LAF damage (during the interseismic stage) increases in simulations with ineffective healing (higher ch).  3.3.2 Fault zone evolution patterns and temporal stability of fault stepovers Strike-slip fault systems evolve over time. As total offset increases, asperities and other frictional barriers to slip are reduced, segment length increases, en-echelon segments coalesce, and the fault straightens and simplifies, reducing overall resistance to slip. Detailed mapping of active and exhumed fault zones (e.g., Tchalenko, 1970; Chester et al., 1993; Sibson, 2003) indicate that the internal structure of fault zones evolves from distributed deformation and bandlimited fractal structures, through localization to principal slip zones, to mature large-scale faults with tabular damage zones and narrow cores of ultra-cataclasite (Lyakhovsky and Ben-Zion, 2009). The structural evolution of fault systems is also characterized by evolving seismicity 75  patterns and fault slip-rate. The slip rate along any new propagating fault by definition must experience a finite period of acceleration, and as it matures and the slip resistance decreases, the fault will accelerate to some steady state slip rate until the tectonic driving force or the fault configuration changes. Ben-Zion and Sammis (2003) describe fault evolution through localization as an eventual consequence of strain weakening rheology, common in crustal materials. However, due to data limitations various aspects of fault evolution are not well understood, and clear distinctions of evolution stages and definitions of mature and immature faults are lacking.  Our previous simulations of damage zone evolution indicated that newly formed and propagating fault segments undergo a very short stage of complexity increase and DOF damage build-up (Finzi et al., 2009). During this early stage strain is distributed within widening damage zones, and fault segments interact to form damaged linking zones. Based on the simulated evolution of the DOF damage width, we previously suggested that this complexity-increase stage ends when the fault has accumulated an offset of about 0.05-0.1 km and is thereafter followed by a long-term state of complexity decrease. Acknowledging that the above time scale is short, we recently analyzed the depth extent of the LAF damage and the slip-rate evolution as additional indicators of fault evolution. We also examined the sensitivity of these fault evolution time scales to variations in the healing parameters. Our analysis indicates that the LAF damage reaches its maximum (co-seismic) depth extent (i.e. the fault core extends to the bottom of the seismogenic zone during large seismic events) after the fault has accumulated an offset of about 0.1-0.2 km, with little sensitivity to the applied healing parameters. Even though the DOF in our simulations fully forms before the LAF damage reaches its maximum extent, the similarity in time scales of these evolution indicators suggests that damage zone width may be used as an indicator for the establishment of the damage zone at depth. The long-term (tectonic) slip-rate in our simulations is first achieved after the first system-size earthquake on the fault (t=1500 years in Fig. 3.6) followed by a significant decrease in slip rate. After the accumulation of approximately 0.1 km of strike-slip offset (~ 3000 model years), the fault experiences ~2000 years of relatively small slip rate fluctuations around the tectonic slip-rate (Fig. 3.6). This occurs when the fault DZ is almost fully established, and it may indicate that the fault is starting to 76  mature (i.e. the fault is transitioning into a state of complexity decrease). After this transition the slip rate experiences several large fluctuations that slowly subside until the rate stabilizes at the long-term rate or until the system is perturbed by a change in fault configuration. Our results indicate that modeled fault systems maintain a near tectonic slip rate only after deformation has localized along narrow slip-zones and after the faults are fully formed down to the bottom of the seismogenic zone. While the healing parameters do not have a significant effect on the time scale of slip rate evolution in our simulations, the complexity of the fault system does affect the magnitude of slip rate fluctuations. The results shown in figure 3.6 are from a simulation of a segmented (complex) fault system with relatively ineffective healing ( ch=0.5). Simulations with a simple through-going fault yield similar time scales of slip rate evolution (i.e. slip rate stabilizes on the tectonic slip rate when the offset reaches ~0.1 km), but they display much smaller slip fluctuations. Simulations with effective healing are expected to evolve into simple fault systems and therefore should also yield relatively small slip rate fluctuations. Short time scales for damage zone and slip rate maturation are supported by various studies of exhumed fault zones that show slip localization within narrow core zones along small faults with minor (< 1 m) offsets (e.g., Evans et al., 2000). The great structural similarity between these faults and large faults with offsets and lengths greater by several orders of magnitude (e.g., the Punchbowl Fault, Chester et al., 1993) indicates that localization to a narrow core zone occurs at a very early stage of shear deformation.  77  Figure 3.6. Slip-rate and damage zone evolution with time. Slip rate averaged over several earthquake cycles (red), geologic (cumulative) slip rate since the formation of the fault system (green), DOF damage width (yellow) and LAF damage maximum (coseismic) depth extent (blue). The damage zone width and depth are normalized to their value at the end of the simulation. These dimensions and the slip rate initially increase rapidly (in association with the first large earthquake), then the slip rate fluctuates around the tectonic rate eventually approaching the long-term geologic rate (32mm/yr). The healing parameters in this simulation correspond toch=0.5, h=30 yr.  The time scale outlined above describes the structural evolution of a single fault segment, and it does not incorporate the effect of large-scale persistent fault complexities that may impede localization. To study this effect we constructed a series of models simulating deformation patterns associated with a stepover zone between two en-echelon strike-slip fault segments. As with our previous simulations (Finzi et al., 2009), we find that large volumes of damage within 78  the stepover zones exhibit ongoing interseismic deformation and distributed inelastic strain. Our current results show that these deformation patterns suppress earthquake propagation and strain localization within the stepover zone resulting in a significantly slower process of complexity reduction and fault system maturation. Eventually, through-going faults may bridge the stepover zone and the tips of the en-echelon fault segments may re-align and link into a contiguous curved fault, however this process is strongly dependent on the effectiveness of healing processes along the fault zone.  In simulations with parameters that favor effective (near-complete) healing and yield low interseismic damage levels, fault stepovers evolve quickly into a single, contiguous fault segment that rapidly retains the characteristics of a mature fault segment (Fig. 3.7). In simulations with settings that favor ineffective (limited) healing, fault stepovers persist, slowing the simplification of the fault geometry (Fig. 3.7). Figure 3.7 compares the evolution of fault stepovers experiencing effective (top;ch=0.4) and ineffective (bottom; ch=0.7) healing, showing snapshots of the fault structure at 5 km depth during the initial 6000 model years after the formation of the stepover (imposed as initial conditions in our simulations). It is apparent from the top panel in Figure 3.7 that with effective healing the en-echelon faults link rapidly and the stepover is readily gapped changing the original structure to a through-going fault that localizes deformation. This active contiguous fault exhibits damage levels significantly higher than its surrounding residual damage. Setting the healing to be inefficient (bottom panel of Figure 3.7) results in a persistent stepover which displays an array of active secondary faults within the original stepover structure. The original en-echelon faults show no sign of evolution after 5000 years (0.16 km of cumulative offset), suggesting that while damage accumulates within the stepover it does not localize deformation on one through-going fault and does not yield a reduction in fault complexity. Extrapolating this pattern to longer time scales leads us to speculate that with very limited healing, stepovers and other fault complexities may only grow with time and they are not expected to become inactive or to yield localization of deformation on a through-going (mature) fault. This evolution pattern would yield long-lasting basins or push-up swells with ongoing internal deformation, as is widely observed along segmented strike-slip faults. 79  Figure 3.7. Damage maps (at a depth of 5 km) displaying two evolving stepover zones along two segmented  strike-slip  faults  representing  effective  (top,ch=0.4)  and  ineffective  (bottom,ch=0.7) healing conditions. While effective healing (top) leads to rapid smoothing of the fault system and turns the stepover into an inactive structure, ineffective healing keeps the stepover structure active and the fault system remains segmented throughout the duration of the simulation. All maps are 70 km long (along strike) and 50 km wide.  3.4  Discussion and conclusions Using an analysis of our damage rheology healing parameters, fault-zone geophysical  observations, and numerical simulations, we have evaluated how material healing affects fault zone complexity, damage zone structure, and long-term fault system evolution.  80  Healing effectiveness controls both the dimensions of the fault damage zone (Fig. 3.5) and the complexity of the fault system in our simulations (Fig. 3.7). This suggests that observable structural characteristics of a fault system and long-term evolution patterns may be linked with fault zone properties and conditions that control gouge-scale healing processes. Establishing such a link would yield a powerful tool to study fault systems where detailed observations of fault zone properties are not available. Furthermore, such a link could help in evaluating the temporal stability of complex fault configurations and local fault structures such as stepovers. Additional postseismic and interseismic observations of fault zone material properties correlated with maps of fault structures and strain distribution within fault systems are needed to study the spatial and temporal characteristics of such a link. It should also be noted that actual interseismic damage levels are not solely determined by material properties (ch, C1 and C2), and the variability of damage levels within simulated and observed fault systems also reflects local variations in fault zone conditions (e.g., rock type, loading geometry, fluid content, temperature). While we can represent such heterogeneities in our simulations, our numerical model does not account for fluid content and rock-fluid interactions. In fact, our healing formulation generalizes the role material properties and conditions play in the healing process, and it does not describe specific strengthening mechanism such as fracture closure, grain rearrangement, compaction, and pressure solution (e.g., Beeler and Hickman, 2004; Yasuhara et al., 2005; Chester et al., 2007).  The fault stepover simulations here and in Finzi et al. (2009) show extensive damage and elasticity degradation sustained over many earthquake cycles. Simulated tensional stepovers exhibit damage patterns consistent with intense tensile fracturing and dilation, and are expected to sustain enhanced permeability. Such damage patterns are consistent with recent structural evolution models for dilational stepovers (De Paola et al., 2007), and with studies of permeability evolution patterns along segmented faults (Micklethwaite and Cox, 2004; Sheldon and Micklethwaite, 2007). We suggest that healing effectiveness has an important effect on the structural properties and temporal stability of fault stepovers, and therefore may have important implications for earthquake propagation and seismic hazard studies. We further suggest that 81  long-lived stepover zones represent unfavorable conditions for healing that should also correlate with significant material weakening down to the lower parts of the seismogenic zone.  Ben-Zion et al. (1999) and Lyakhovsky et al. (2001) suggested that the ratio between healing and loading timescales controls fault zone complexity, where relatively fast healing and slow loading results in a disordered fault system. However, our results indicate that for fault zone rocks with healing timescales that vary by one order of magnitude, and even less variability in loading time scales, healing effectiveness plays a more dominant role in determining the properties of the evolving damage zones. In our simulations, effective healing (associated with low inter-seismic damage levels along fault segments) produces wide, shallow damage zones. However, our simulations indicate that faults with effective healing rapidly supersede preexisting stepovers, whereas ineffective healing (high residual damage levels) produces fault systems with long-lived stepovers. In the effective healing case, a stepover is permanently weak to only a shallow depth. Below a few kilometers, confining stresses are large and the fault zone heals rapidly to the (low) characteristic damage level. The stepover has less of an influence on the propagation of large ruptures, which break through it, eventually allowing a through-going fault to form. Residual damage persists in this case only in the top few kilometers of the shallow, inactive stepover zone. Our results, based on fully 3D simulations and a wide range of healing parameters, provide important refinements to the earlier general findings of Lyakhovsky et al. (2001) and Ben-Zion et al. (1999), which were based on simpler 2D simulations. The current study also improves on the analysis of Finzi et al. (2009), by demonstrating that co-varying both healing parameters (C1 and C2 of equations 3.2-3.5) to obtain a large range of values of the characteristic residual damage level ch, profoundly affects the modeled damage zone structure and fault system evolution.  While possible links between observable fault structures and long-term deformation patterns require additional support from fault zone observations and further numerical analysis, the following features of fault zone structure and evolution appear to be robust in our simulations and are well supported by field observations:  82    Fault zone healing is very rapid with most (>50%) of the healing occurring within hours to days after an earthquake. In addition, fault zones exhibit a wide range of quasi-steady state interseismic damage levels, with  between (approximately) 0.25 and 0.65 (i.e., a rigidity reduction of 25%-65%, and shear wave velocity reduction of 15%-40%).    Healing effectiveness is an important factor in the evolution of fault damage zones. Conditions favoring complete healing yield wider DOF damage in the top few kilometers of the crust and shallower fault-core LAF damage, whereas conditions favoring ineffective (limited) healing yield narrower DOF damage and deeper LAF damage zones.    Healing effectiveness is an important factor in the temporal stability of fault stepovers. Fault complexities such as stepover zones are expected to be shorter-lived if fault zone conditions favor more complete healing.    Active fault stepovers typically consist of damaged material down to the lower seismogenic zone, and should be detectable seismically and geodetically. Obtaining more information regarding material properties within stepover zones would aid future analysis of fault zone healing and fault system evolution.  83  Bibliography Andrews, D.J., 2005. Rupture dynamics with energy loss outside the slip zone. J. Geophys. Res., 110, B01307, doi:10.1029/2004JB003191. Barbot, S., Fialko, Y., Sandwell, D., 2009. Three-dimensional models of elasto-static deformation in heterogeneous media, with application to the Eastern California Shear Zone. Geophys. J. Int., 179, 500-520, DOI: 10.1111/j.1365-246X.2009.04194.x Beeler, N.M., Hickman, S.H., 2004. Stress-induced, time-dependent fracture closure at hydrothermal conditions, . J. Geophys. Res. 109, B02211, doi: 10.1029/2002JB001782. Ben-Zion, Y., Dahmen, K., Lyakhovsky, V., Ertas D., Agnon, A., 1999. Self-driven mode switching of earthquake activity on a fault system. Earth Planet. Sci. Lett., 172/1-2, 1121. Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J.G., Ozer, N., Michael, A.J., Baris, S., Aktar, M., 2003. A shallow fault zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey. Geophys. J. Int. 152, 699-717. Ben-Zion, Y., Sammis, C., 2003. Characterization of fault zones. Pure and Appl. Geophys. 160, 677–715. Ben-Zion, Y., Shi, Z., 2005. Dynamic rupture on a material interface with spontaneous generation of plastic strain in the bulk. Earth Planet. Sci. Lett., 236, 486-496, DOI: 10.1016/j.epsl.2005.03.025. Ben-Zion, Y., Lyakhovsky, V., 2006. Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology. Geophys. J. Int. 165, 197-210. Bhat, H.S., Dmowska, R., King, G.C.P., Klinger, Y., Rice, J.R., 2007. Off-fault damage patterns due to supershear ruptures with application to the 2001 Mw 8.1 Kokoxili (Kunlun) Tibet earthquake. J. Geophys. Res., 112, B06301, doi:10.1029/2006JB004425. Boness, N.L., Zoback, M.D., 2004. Stress-induced seismic velocity anisotropy and physical properties in the SAFOD Pilot Hole in Parkfield, CA. Geophys. Res. Lett., 31, L15S17, doi:10.1029/2004GL019020. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81-97. Chester, F.M., Evans, J.P., Biegel, R.L., 1993. Internal structure and weakening mechanisms of the San Andreas Fault. J. Geophys. Res. 98, 771–786. Chester, F.M., Chester, J.S, Kronenberg, A.K., Hajash, A., 2007. Subcritical creep compaction of quartz sand at diagenetic conditions: Effects of water and grain size. J. Geophys. Res. 112, B06203, doi: 10.1029/2006JB004317. Cochran, E.S., Li, Y.G., Shearer, P.M., Barbot, S., Fialko, Y., Vidale, J.E., 2009. Seismic and geodetic evidence for extensive, long-lived fault damage zones. Geology, v. 37, no. 4, p. 315-318, DOI: 10.1130/G25306A.1.  84  De Paola, N., Holdsworth, R.E., Collettini, C., McCaffrey, K.J.W., Barchi, M.R., 2007. The structural evolution of dilational step-overs in regional transtensional zones, in: Cunningham W.D., Mann, P. (Eds), Tectonics of Strike-Slip Restraining and Releasing Bends. Geological Society, London. Special Publication., 290, pp. 433-445. Dieterich, J.H., 1978. Time-dependent friction and the mechanics of stick-slip. Pure and Appl. Geophys. 116, 790-805. Dieterich, J.H., Kilgore, B.D., 1996. Imaging surface contacts: power law contact distributions and contact stresses in quarts, calcite, glass, and acrylic plastic. Tectonophysics, 256, 219-239. Dor O., Ben-Zion, Y., Rockwell, T.K., Brune, J.N., 2006. Pulverized rocks in the Mojave section of the San Andreas Fault Zone. Earth Planet. Sci. Lett. 245, 642-654 doi:10.1016/j.epsl.2006.03.034. Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., Duman, T.Y., 2008. Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction. Geophys. J. Int., 173, 483–504, doi: 10.1111/j.1365246X.2008.03709.x. Evans, J.P., Shipton, Z.K., Pachell, M.A., Lim, S.J., Robeson, K., 2000. The structure and composition of exhumed faults, and their implication for seismic processes, In Proc. of the 3rd Confer. on Tectonic Problems of the San Andreas System, Stanford University. Fialko, Y., Sandwell, D., Agnew, D., Simons, M., Shearer, P., Minster, B., 2002. Deformation on nearby faults induced by the 1999 Hector Mine earthquake. Science, 297, 1858-1862. Finzi Y., Hearn, E.H., Lyakhovsky, V., Ben-Zion, Y., 2009. Structural properties and deformation patterns of evolving strike-slip faults: Numerical simulations incorporating damage rheology. Pure and Appl. Geophys., 166, 1537-1573, doi: 10.2007/s00024-0090522-1. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y., Lockner, D., 2004. A Visco-elastic damage model with applications to stable and unstable fracturing. Geophys. J. Int. 159, 1155– 1165. Hamiel, Y., Fialko, Y., 2007. Structure and mechanical properties of faults in the North Anatolian Fault system from InSAR observations of coseismic deformation due to the 1999 Izmit (Turkey) earthquake. J. Geophys Res. 112,B07412, doi:10.1029/2006JB004777. Hearn, E.H., Fialko, Y., 2009, Can compliant fault zones be used to measure absolute stresses in the upper crust?, J. Geophys. Res., 114, B04403,DOI: 10.1029/2008JB005901. Jaeger, J.C., Cook, N.G.W., 1979. Fundamentals of Rock Mechanics, 3rd ed. Chapman and Hall, New York, 593 p. Johnson, P.A., Jia, X., 2005. Nonlinear dynamics, granular media and dynamic earthquake triggering. Nature, 473, 871874. Karabulut, H., Bouchon, M., 2007. Spatial variability and non-linearity of strong ground motion near a fault. Geoph. J. Int. 170, 1, 262-274. 85  Kim, Y.S., Peacock, D.C.P., Sanderson, D.J., 2004. Fault damage zones. J. of Struct. Geology, v. 26, p. 503-517. Lewis, M.A., Peng, Z., Ben-Zion, Y., Vernon, F., 2005. Shallow seismic trapping structure in the San Jacinto fault zone. Geophys. J. Int. 162, 867–881, doi:10.1111/j.1365246X.2005.02684.x, 2005. Li, Y.-G., Aki, K., Adams, D., Hasemi, A., Lee, W.H.K., 1994. Seismic guided waves trapped in the fault zone of the Landers, California earthquake of 1992. J. Geophys. Res. 99(B6), 11,705–11,722. Li, Y.-G., Chen, P., Cochran, E.S., Vidale, J.E., Burdette, T., 2006. Seismic evidence for rock damage and healing on the San Andreas Fault associated with the 2004 M 6.0 Parkfield Earthquake. Bull. Seism. Soc. Am. 96, no. 4B, S349–S363. Lockner, D., Walsh, J., Byerlee, J., 1977. Changes in seismic velocity and attenuation during deformation of granite. J. Geophys. Res., 82, 5374-5378. Lyakhovsky, V., Ben-Zion, Y., Agnon, A., 1997a. Distributed damage, faulting, and friction. J. Geophys. Res. 102, 27,635-27,649. Lyakhovsky, V., Reches, Z., Weinberger, R., Scott, T.E., 1997b. Non-linear elastic behavior of damaged rocks. Geophys. J. Int. 130, 157-166. Lyakhovsky, V. Ben-Zion, Y., Agnon, A., 2001. Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere. J. Geophys. Res. 106: 4103-4120. Lyakhovsky, V., Ben-Zion, Y., Agnon, A., 2005. A viscoelastic damage rheology and rate- and state-dependent friction, Geophys. J. Int. 161, 179-190. Lyakhovsky, V., Ben-Zion, Y., 2008. Scaling relations of earthquakes and aseismic deformation in a damage rheology model. Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365246X.2007.03652.x. Lyakhovsky V., Ben-Zion, Y., 2009. Evolving fault zone structures in a damage rheology model. in review, Geochem., Geophys., Geosyst. Ma, S., 2008. A physical model for widespread near-surface and fault zone damage induced by earthquakes. Geochem. Geophys. Geosyst., 9, Q11009, doi:10.1029/2008GC002231. Micklethwaite, S., Cox, S.F., 2004. Fault-segment rupture, aftershock-zone fluid flow and mineralization. Geology, 32, 813–816. Micklethwaite, S., Cox, S.F., 2006. Progressive fault triggering and fluid flow in aftershock domains: Examples from mineralized Archaean fault systems. Earth and Planet. Sci. Lett., Vol. 250, Issues 1-2, 318-330. Peng, Z., Ben-Zion, Y., Michael A.J., Zhu, L., 2003. Quantitative analysis of seismic trapped waves in the rupture zone of the 1992 Landers, California earthquake: Evidence for a shallow trapping structure. Geophys. J. Int. 155, 1021-1041. Peng, Z., Ben-Zion, Y., 2006. Temporal changes of shallow seismic velocity around the Karadere-Duzce branch of the north Anatolian fault and strong ground motion. Pure and Appl. Geophys. 163, 567-600, DOI 10.1007/s00024-005-0034-6. 86  Revenaugh, J., 2000. The relation of crustal scattering to seismicity in Southern California. J. Geophys. Res., 105, 25,403–25,422. Rockwell, T.K., Ben-Zion Y., 2007. High localization of primary slip zones in large earthquakes from paleoseismic trenches: Observations and implications for earthquake physics. J. Geophys. Res. 112, B10304, doi:10.1029/2006JB004764. Sheldon H. A., Micklethwaite S., 2007. Damage and permeability around faults: Implications for mineralization. Geology: Vol. 35, No. 10 pp. 903–906. Sibson, R.H., 2003. Thickness of the seismic slip zone. Bull. Seismol. Soc. Am. 93, 3, 11691178. Sylvester, A.G., 1988. Strike-slip faults. Geol. Soc. Am. 100, 1703-1777. Tchalenko, J.S., 1970. Similarities between shear zones of different magnitudes. Geol. Soc. of Am. Bull. 81, 1625–1640. Wesnousky, S., 1994. The Gutenberg-Richter or characteristic earthquake distribution, which is it? Bull. Seismol. Soc. Am. 84, 1940–1959. Wu, C., Peng, Z., Ben-Zion, Y., 2009. Non-linearity and temporal changes of fault zone site response associated with strong ground motion. Geophys. J. Int., 176, 265–278, doi: 10.1111/j.1365-246X.2008.04005.x. Yasuhara, H., Marone, C., Elsworth, D., 2005. Fault zone restrengthening and frictional healing: the role of pressure solution. . J. Geophys. Res. 110, B06310, doi: 10.1029/2004JB003327.  87  4  EFFECT OF LATERAL RHEOLOGICAL CONTRASTS ON DAMAGE ZONE  STRUCTURE  AND  SURFACE  DEFORMATION  AROUND  STRIKE-SLIP FAULTS3  Summary Numerical simulations incorporating damage evolution in the crust and various rheologies at depth are used to explore the effects of lateral contrasts in lithosphere rheology on the growth and evolution of strike-slip fault systems. Using long-term fault evolution models, we address recent hypotheses regarding asymmetric interseismic velocity profiles, the effects of crustal elasticity contrasts on fault propagation and damage distribution. We assess whether rheological contrasts cause asymmetric shear and damage zone structure, whether evolving strike-slip faults migrate toward and propagate along rheological interfaces, and how these effects may influence interseismic surface velocities. We find that reasonable contrasts in lithosphere elasticity have only a minor effect on the symmetry of damage or the propagation of faults. Contrasts in effective plate thickness or viscosity structure attract propagating faults, which tend to locate about 1 to 5 km from the interface, in the thinner plate (or the side with lower viscosities). Brittle damage is more extensive on the thin-plate (or low-viscosity) side of the fault. In the lower crust and mantle, viscosity contrasts induce the formation of broad but high strain-rate shear zones, which may be offset from the material interface by up to 10 km. In addition, our results demonstrate how models with reasonable contrasts in viscoelasticity explain marked asymmetry in interseismic GPS velocities (for example, across the North Anatolian Fault in the Marmara Sea), and that extreme contrasts in elasticity are not required. Finally, we show that models of asymmetric interseismic deformation in which the fault trace is assumed to coincide with the material interface may lead to an erroneous interpretation of GPS surface velocity data.  3  A version of this chapter will be submitted for publication. Finzi, Y. and Hearn, E.H. Effect of lateral rheological contrasts on damage zone structure and surface deformation around strike-slip faults.  88  4.1 Introduction Major faults often juxtapose materials that have different physical properties. Some strike-slip faults, such as the San Andreas, Alpine and Philippine faults, juxtapose oceanic lithosphere and continental lithosphere, and are associated with material contrasts over very large (lithosphere) scales (Le Pichon et al., 2005). For large offsets, such large-scale material interfaces could be due to the action of the fault itself, moving different materials large distances along strike. However, such interfaces may occur with small total offsets if propagating ruptures migrate toward bi-material interfaces (Brietzke and Ben-Zion, 2006; Armijo et al., 1998; Sengor et al., 2005; Lyakhovsky et al., 1994). A smaller-scale bi-material interface may result from large damage zones generated by repeated earthquakes, because faults tend to localize on one side of such damage zones (e.g., Ben-Zion and Andrews, 1998; Ben-Zion and Sammis, 2003). Geodetic measurements of deformation around strike-slip faults appear consistent with such asymmetry and deviates from the predictions of uniform or layered elastic dislocation models (Meade et al., 2002; Savage et al., 2004; Le Pichon et al., 2005; Fialko, 2006; Wdowinski et al, 2007). Some studies have suggested that the observed, asymmetric deformation patterns may result (or partly result) from the existence of intensive damage along large fault zones (Le Pichon et al., 2003; Fialko, 2004). Jolivet et al. (2008) and Le Pichon et al. (2003) suggest that along certain sections of the Altyn Tagh Fault in Tibet, and the North Anatolia Fault Zone in Turkey, the ratio of rigidity of the juxtaposed rocks is between 5 and 10 (based on the assumption that the asymmetry is due to a contrast in lithosphere rigidity). As this ratio is far greater than expected rigidity contrast across continental strike-slip faults (Ben-Zion and Andrews, 1998; Hauksson, 2000), other conditions probably contribute to the asymmetric pattern of deformation. Such conditions may include contrasts in viscosity across the faults, variations in the effective plate thickness, effects of an asymmetric gouge and damage zone, variations in fault dip, and enhancements of the rigidity ratio by dynamic effects (e.g., Li and Rice, 1987; Lisowski et al., 1991; Fialko et al., 2001, 2005; Schmalzle et al., 2006). To investigate how material interfaces with variable rigidity, viscosity and effective elastic plate thickness affect surface deformation, damage distribution and fault migration we perform a series of numerical simulations applying damage rheology. In the following sections we test three central hypotheses regarding the effects that bi-material interfaces have on the 89  structure and deformation patterns of strike-slip faults. Our simulations test whether fault-zone structures along various types of material interfaces are symmetric or asymmetric (section 4.3.1), and whether the evolving faults tend to migrate toward and propagate along the interfaces (section 4.3.2). Simulations are also used to characterize the surface velocity profiles across strike-slip faults located along material interfaces (section 4.3.3).  4.2 Methods and model set-up 4.2.1 Numerical model and damage rheology We apply damage rheology to simulate evolving strike-slip faults in a 3D layered model domain. Our numerical model (Lyakhovsky et al., 1997a, b, 2001) is based on the emerging interdisciplinary view that the brittle portion of the Earth's lithosphere is “damaged” in the sense that it contains an internal distribution of microcracks, joints and faults. During continuing deformation (e.g. shearing), the material around such structures is subjected to large stress concentrations which necessarily lead to increased cracking (damage accumulation) and a reduction of the elastic strength. Subsequently, when shear stresses drop (e.g., after a seismic failure of a fault), the internal damage undergoes healing which increases the elastic strength of material. This model provides a self-consistent physical framework for studying the coupled evolution of fault structures and earthquakes in an upper crust that exhibits inelastic processes such as creep and off fault damage, and that is affected by viscoelastic processes in the underlying lower crust and mantle. In this model, an evolving damage variable () quantifies the extent of internal damage and modifies the elastic properties of crustal material, simulating the creation and healing of faults. This scalar variable () reflects the extent of strength degradation due to crack formation and opening, with = 0 indicating undamaged rock and = 1 corresponding to a complete strength loss. The damage rheology and our numerical model are briefly described below. For further details and selected applications of this model we refer to previous publications (Lyakhovsky et al., 1997a, b, 2005; Hamiel et al., 2004; Ben-Zion and Lyakhovsky, 2006; Finzi et al., 2009).  90  To account for the evolution of material properties, the damage rheology model introduces a third elastic modulus () and makes the elastic rigidity () a function of the damage state variable, as follows:    = constant;     ·m0–) ≈ 0    m    where  and  are the Lamé parameters of linear Hookean elasticity, and  = I1/√I2 is the strain invariants ratio (I1=kk and I2=ijij are the first and second invariants of the elastic strain tensor  ij). The parameter 0 separates states of material degradation ( > ) and healing ( < ) associated with positive and negative damage evolution, respectively. As the damage variable  increases, the shear modulus  decreases, Poisson‟s ratio  increases, and the modulus  increases from 0 (damage free) to m, amplifying the non-linearity of rock elasticity.  A typical model setup consists of a layered seismogenic crust governed by damage rheology, underlain by a viscoelastic upper mantle governed by a power-law rheology. The material properties in our simulations reflect the common lithologies of the continental lithosphere (e.g., quartz-diorite (upper crust), diabase (lower crust) and olivine (upper mantle)), and simplified Newtonian rheologies representative of the mantle (Table 4.1). The boundary conditions apply a variable force on the sides and bottom of the model domain to simulate a lithospheric-scale plate boundary outside the model domain with a variable degree of strain localization where the plate boundary intersects the model domain, and far-field plate velocities comparable to the San Andreas Fault (Lyakhovsky et al., 2001). The driving forces are proportional to the mismatch (slip-deficit) between the far field plate motion and displacement of the boundary nodes, multiplied by a boundary stiffness parameter (Finzi et al., 2009). Unlike the constant velocity or constant stress conditions, our variable force boundary condition accounts for stress accumulation and strain-rate decrease during the interseismic periods, as well as abrupt stress reduction at the model boundary during seismic events. 91  4.2.2 Model set-up To investigate the lithospheric conditions in which asymmetric deformation patterns form, we simulate fault evolution in model domains with five types of material interfaces. In these simulations we investigate the effects of contrasts in rheology, which we impose by varying both rheologies and crustal thickness. In addition to idealized elastic and Newtonian viscoelastic materials with a range of viscosities, we model several experimentally determined power-law rheologies for the crust and mantle. The material parameters for these flow laws are given in Table 4.1, and profiles of effective viscosities produced by these flow laws are shown in Figure 4.1 (for high and low strain rates). Simulated effective viscosity values in our models are spatially variable (e.g., lower in areas of high strain rates), and do not vary simply with depth and temperature, as in Figure 4.1. Results from simulations with lateral rheological contrasts are compared with deformation patterns in laterally homogeneous simulations. The latter simulations involve a layered model with a 3 km thick sediment layer (material S in Table 4.1), diabase crust down to 20 or 40 km (material C in Table 4.1) and a Newtonian mantle (material M1 in Table 4.1) extending down to the bottom of the model domain (z=50 Km). In these simulations the boundary conditions induce shear localization (and fault nucleation) near the center (x=40 km) of the front/back faces of the domain boundary (at y= 0 and y=150 km).  Five model categories incorporate the following types of lateral contrasts:  A. Material interface in the crust. In these simulations the entire crust (modeled to 20 km depth) is represented by a diabase flow law, but a rigidity contrast is imposed (that is, materials C1 and C2 are juxtaposed, see Table 4.1). The mantle is modeled as laterally homogeneous (using either material M4 or M2 in Table 4.1). In other simulation categories (AB and D), we juxtapose quartz-diorite and diabase rheologies in the crust, resulting in a contrast in both the rigidity and the effective viscosity. The resulting rheological contrast is described below.  92  B. Material interface in the mantle. For simple, idealized models, the mantle is treated as Newtonian, and mantle viscosities of 1019 Pa S and 1023 Pa S are juxtaposed (i.e., materials M1 and M2 in Table 4.1). In these models the crust (material C in Table 4.1) is 35 km thick. A second suite of models in this category juxtaposes wet and dry olivine rheologies in the mantle (M3 and M4, Table 4.1). In these models the crust is 20 or 35 km thick. Figure 4.1 shows that these power law rheologies should produce a significant contrast in effective viscosity across the interface (1021 Pa S and 2 x 1022 Pa S at the Moho). We indicate the crust thickness and the rheologies for specific models whose results are shown in the figures, and address how parameter variations influenced results in the discussion section.  AB. Material contrasts in both the crust and the mantle. A crustal rigidity contrast (between materials C1 and C2) and Newtonian mantle rheologies (M1, M2) were assumed in one suite of models of this category. In another suite of models, lateral contrasts in power-law crust (C1, C3) and mantle (M3, M4) rheologies were modeled. At a low strain rate (10-15 /s) the diabase rheology in our simulations (C-C2) displays a high effective viscosity (above 1023 Pa s) to a depth of 22 km, tapering to 1020 Pa s at 30 km depth (Fig. 4.1, note that the effective viscosities are lower for higher strain rates). The quartz-diorite rheology (C3) produces effective viscosities that are about 20 times smaller at a given depth (Fig. 4.1). Juxtaposing rheologies C1 and C3 produces a 3-4 km step in the effective elastic plate thickness, where the “plate” is defined as crust with an effective viscosity greater than 1023 Pa s or a differential stress exceeding 100 MPa. (This assumes a strain rate of 10-15 /s. Since strain rate varies spatially in our models, the lateral rheological transition will be more complex). Rheologies we use in specific AB models are indicated where results from these models are shown and discussed.  C. An interface in crust thickness, juxtaposing a thin crust (Moho at 15 km) against a thick crust (Moho at 35 km). The crust is modeled with material C. The mantle in these simulations was chosen to be a low-viscosity, Newtonian rheology (M2 in Table 4.1). The effective elastic thickness is 15 km on the thin-crust side and about 22 km on the side 93  with the thicker crust. The thickness and density of the sedimentary layer were calculated with a condition of lithostatic equilibrium across the interface.  D. A weakened crustal fault-zone and low-viscosity mantle shear zone embedded in the lithosphere. The width of this zone was set at 2 km in the crust and 5 km in the mantle, so that it is represented by 1 to 3 model elements. The crust in these simulations is 35 km thick. Rheologies used in these simulations were of types C3 and M2 for the low-viscosity zone, and C and M1 for the crust and mantle outside this zone. The shear zone effective viscosity in the mantle is 104 times lower than that of the surrounding mantle, and in the crust it is lower than that of the surrounding crust by about a factor of 20. In the seismogenic crust (top 15-20 km), the C and C3 flow laws both yield very large viscosity values, so brittle failure dominates, and the rigidity contrast becomes more important than the viscosity contrast.  In all simulations we set the boundary conditions to simulate relative plate motions of ~30 mm/yr. In models we run to characterize the damage zone structure and the surface velocity field around faults at material interfaces (Section 4.3.1), we align the interface with the fault trace predicted by models that were run without lateral contrasts. In many of these simulations, the boundary conditions concentrate strain near the centers of the north, south, and bottom model boundaries. This generally results in a single, long fault segment, rather than a segmented fault with numerous steps and bends. The material interface in such simulations is situated along the line connecting the fault nucleation sites on the model boundaries. In other simulations we apply boundary conditions that induce a linear velocity gradient on the north and south domain boundaries, allowing fault nucleation and localization to occur anywhere within the domain. As with the models with localized strain at the model boundaries, we align the material interface along the longest fault segment that formed in the homogeneous simulation. To investigate the conditions under which a fault zone might propagate to and then along a rheological contrast (Section 4.3.2), we vary the orientation of the various interfaces and their proximity to locations of imposed, high strain rates (i.e., seeded faults) at the model boundaries. In these simulations we typically initialize the simulation with the material interface misaligned 94  with the expected trace of the fault (based on simulations with homogeneous lithosphere), offset by 1-10 km from the expected fault nucleation sites at the model boundaries, or we apply boundary conditions that do not localize strain (allowing fault nucleation to occur anywhere within the domain).  ID  description  E (GPa)    A (MPa-nS-1)  n  Q (KJmol-1)  S  sediment  40-50  0.3  1.00E-13  1  0  C  crust  80  0.3  0.126  3.05  276  C1  „strong‟ crust  88-92  0.3  0.126  3.05  276  C2  „weak‟ crust 1  60-72  0.3  0.126  3.05  276  C3  „weak‟ crust 2  60-80  0.3  1.25  2.4  212  M1  Newtonian „strong‟ mantle  150-180  0.3  1.00E-17  1  0  M2  Newtonian „weak‟ mantle  100-150  0.3  1.00E-13  1  0  M3  Power-law „strong‟ mantle  150-180  0.3  7.00E+04  3  520  M4  Power-law „weak‟ mantle  100-150  0.3  1.90E+03  3  420  Table 4.1 Material properties used in our numerical simulations. To test the effect of a range of rigidity ratios we varied the Young‟s Modulus of crust and mantle rheologies as shown in the table. Crustal material C is based on a diabase flow law (Carter & Tsenn, 1987), crustal material C3 is based on a quartz-diorite flow law (Carter & Tsenn, 1987), mantle material M3 is based on an upper mantle (olivine) flow law (Kirby and Kronenburg, 1987), and mantle material M4 is based on a wet olivine flow law (Karato et al., 1986).  95  Figure 4.1 Effective viscosity profiles for wet and dry crust (a) and mantle (b), for flow laws used in our models (Table 4.1). A geothermal gradient of 20 degrees per km is assumed, and profiles for two different strain rates are shown for each rheology. Our models yield spatially heterogeneous strain rates, and at a given depth, effective viscosities will be lower near the fault zone.  96  4.3 Results  We simulate fault (and associated damage zone) evolution for the five categories of laterally heterogeneous models described above. In most cases, the main fault evolves in the softer material, adjacent to the material interface. Damage and elevated shear strain rates at depth are primarily confined to the softer side of the interface, resulting in asymmetric surface deformation.  4.3.1 Damage zone asymmetry and fault location relative to the material interface Because long fault segments follow the material interface in these models, the damage zone structure is fairly consistent along strike for distances that are long relative to the width of the damage zone. Still, the distribution of shallow damage along natural fault zones such as the Eastern California Shear Zone does vary to some degree along strike (Cochran et al 2009, Barbot et al., 2009) and our earlier models show such variations as well (Finzi et al., 2009). Our model elements may fail in small earthquakes at any time, leading to local heterogeneities in damage levels at a particular snapshot in time. To minimize these effects, we average several profiles to characterize surface velocity profiles and lateral offsets between the fault core and the material interface (Table 4.2). To quantify the asymmetry of the damage zone we compare the extent of distributed shallow damage on (1) the weak and strong sides of the material interface, and (2) either side of the fault core. The damage zone widths on the weak and strong sides (Dw and Ds, respectively) are used to define Si, a ratio of damage symmetry relative to the material interface (Si=Dw/(Dw+Ds)). Si is 1 if all the damage is on the weak side of the interface, and 0.5 if the damage distribution is symmetric (relative to the interface). As the fault core does not typically coincide with the material interface, the ratio defined above does not indicate whether the damage is symmetric relative to the fault core. For this we also evaluate the damage zone symmetry relative to the fault core, Sfc = Dl/(Dl+Dr), where symmetry ratios of Sfc>0.5 indicate that more damage evolved away from the interface (on the left side of the faults in Figure 4.2), and symmetry ratios of Sfc<0.5 indicate that more damage evolved on the faults side nearer to the 97  interface (right side of the fault in Figure 4.2). In Table 4.2 we indicate the range of symmetry ratios and offset between the fault core and the material interface for each category of model.  Figure 4.2 Determining damage zone asymmetry measures Si (a) and Sfc (b). The shaded region indicates the stronger material. O is fault offset relative to the material interface.  98  symmetry interface type  relative to interface (Si)  A  symmetry relative to fault offset (km) core (Sfc)  0.55-0.6  0.5-0.55  B  0.6-0.75  0.3-0.5  (mantle interface)  (0.55-0.85)  (0.15-0.6)  0.85  0.35-0.5  (0.65-1)  (0.3-0.55)  0.85  0.55-0.6  (0.7-1)  (0.45-0.75)  0.5-0.55  0.5-0.55  (crustal interface)  AB  C (contrast in crustal thickness)  0-2  Crustal rigidity contrasts of ~20% have limited effect on damage zone (DZ) structure (Fig. 4.3a).  0-5  The fault zone and distributed damage evolve on the low viscosity side of the interface (Fig. 4.3b, 4.3c).  2-10  The fault zone and distributed damage evolve on the weak side of the interface (Fig. 4.3d).  2-5  The fault zone and widely distributed damage evolve on side with the thin crust (Fig. 4.3e).  0  Values are similar to those of a fault in a homogeneous model (Fig. 4.3f).  D (fault and shear zone)  comments  Table 4.2 Asymmetry of fault damage zones along material interfaces. Where large variations in the symmetry ratio are observed along strike, the range of the measured ratio is given in parentheses. Measured offsets indicate the horizontal distance between the main fault zone (at seismogenic depth) and the material interface.  Figure 4.3 shows representative cross-sections showing damage distribution along simulated faults in all model categories. Bi-material interfaces near or along active faults can induce asymmetric distribution of damage relative to the fault core and the material interface (Table 4.2). In simulations with no viscosity contrasts but with a crustal rigidity contrast (category A models), only a subtle damage zone asymmetry was observed, even where significant rigidity contrasts were applied (22.5 GPa and 37.5 GPa for the juxtaposed crustal 99  materials). The damage zone was slightly more asymmetric (Sfc reached 0.55-0.6) in models with extreme rigidity contrasts (15-19 GPa and 45-60 GPa). Our results were insensitive to variations in crustal rigidity within the more realistic parameter range (i.e., with rigidities stated in Table 4.1). For the model shown in Figure 4.3a, the modeled rigidities were 27 GPa and 35 GPa. Category B simulations (for both Newtonian and power-law rheologies) yielded asymmetric damage zones with Si values of 0.55-0.85, and with along-strike variations in the location of the fault core relative to the shallow damage and the material interface (S fc values of 0.15-0.6; see Fig. 4.3b and 4.3c). We find that in simulations where the fault core is significantly offset from the material interface (simulations of type B and AB), the damage zone is asymmetric relative to both the fault core and the interface as the shallow damage typically evolves between the fault core and the material interface (e.g., Fig. 4.3c). An important exception was found in some category AB simulations with power-law crustal rheologies (C1, C3, M3, and M4). In these simulations, several sections of the fault were offset from the interface by 6-10 km, and displayed near-symmetric damage distributions around the fault core (Si=0.8-1, Sfc=0.45-0.55; e.g., Fig. 4.3d). Simulations of category AB without a viscosity contrast in the lower crust (C1, C2, M1, M2) formed damage zones centered between the fault core and the interface and displayed smaller offsets (2-6 km), lower Si values (0.65-0.8) and lower Sfc values (0.3-0.5). For the model shown in Figure 4.3d (type AB), the modeled rigidities were 27 GPa and 33 GPa. Category C simulations typically formed a fault core 2-4 km from the interface (on the thin crust side) with shallow damage being widely distributed and asymmetric relative to the fault and interface (Si=0.8-0.9, Sfc=0.5-0.6; e.g., Fig. 4.3e). Simulations with an embedded localized fault and shear zone (category D models) produce symmetric damage zones which were not offset relative to the fault core.  100  Figure 4.3 Damage zone structures along various types of material interfaces. Cross-sections of damage level along interfaces that constitute: (a) crustal rigidity contrast (type A in Table 4.2), (b, c) mantle viscosity and rigidity contrasts (type B), (d) rigidity and viscosity contrast in crust and mantle (type AB), (e) contrast in lithospheric structure (type C), (f) embedded fault and shear zone with low rigidity and viscosity (type D). All cross-sections display an area 18 km wide and 26 km deep, and dashed lines indicate the location (x coordinate) of the interfaces. The wide damage distributions at depth (panels d and f) are not realistic fault features, they are partly due to extrapolation (in plotting) and a coarse grid used in these simulations (2.5-4 km wide). 101  4.3.2 Fault propagation toward and along bi-material interfaces Our simulations of category A with a crustal rigidity interface (20-40% contrast) indicate that faults that nucleate along such interface will tend to propagate along it as a simple continuous fault. Comparing such a simulation (Fig. 4.4a) with a homogeneous but otherwise similar simulation (Fig. 4.4b), shows that propagation along such interface may suppress the formation of stepovers that would offset the fault from the interface (in the absence of a material interface). In the homogeneous simulation (Fig. 4.4a) the faults that nucleate on the north and south boundaries of the domain form a stepover linking zone rather than linking directly as observed in the simulation with material interface (Fig. 4.4b). Other simulations (not in Figure 4.4) show that faults will: (a) diverge from such an interface where it is not aligned with the fault trace determined by regional stress (even if the fault nucleated on the interface and the interface deviates by only 15° from the its eventual strike), (b) not change the direction in which they propagate where they intersect such an interface, (c) not propagate to such an interface where it is 5-7 km from their nucleation site (even if the interface is aligned with the faults propagation track). A simulation consisting of an extreme rigidity contrast (with juxtaposed rigidities of 15 and 60 GPa) displayed fault evolution patterns consistent with the above observations (i.e. did not form a fault constrained to evolve along the interface). The variety of results described above lead us to conclude that a material interface consisting of a large (but realistic) crustal rigidity contrast does not have a significant effect on the geometrical configuration of evolving strike slip fault systems, and would not attract or modify the propagation path of individual faults. In simulations with a material interface in the mantle (categories B and AB) we do not observe fault propagation to the material interface, except for minor readjustment of the fault where it is nucleated near the interface and on the high viscosity high rigidity side of it (Fig. 4.4d). However we do observe greater effect on the long term evolution of fault systems. The locus of deformation along segmented faults with distributed damage and faults that propagated sub-parallel to the interface seem to migrate (over time scales of several thousands of years) toward the interface, leaving other parts of the damage zone inactive (Fig. 4.4c). Such evolution patterns would also explain the wide distribution of shallow damage and the large range of offset between the fault core and the interface observed in these simulations (Fig. 4.3b-d, Table 4.2). Once a fault zone reaches such an interface (type B or type AB) which is aligned with its general 102  propagation path, the fault and damage zone will evolve along the interface with strain concentrated adjacent to the mantle interface (on the low viscosity side) and shallow damage distributed above the interface (type B) or on the side with the weaker crust (type AB). In simulations with a contrast in crustal thickness (category C) evolving faults propagate to and along the material interface. This is observed in simulations where faults were nucleated on both sides of the interface at an offset of 5-7 km from the interface (Fig. 4.3e). In simulations where strain was not localized on the model boundaries, the faults nucleated near the lithospheric interface and propagated adjacent to the interface (on the side with the thin crust). Similar fault propagation patterns are observed in simulations with embedded fault and shear zones (type D). Figure 4.4f shows a fault propagating along such interface that is not aligned with the fault trace observed in homogeneous simulations (not in Fig. 4.4). In this simulation the interface guided fault forms a right stepping, compressional stepover where the homogeneous simulation produces only left-stepping tensional stepovers. Figure 4.4f also shows how the faults propagated 7-15 km from their nucleation sites to the material interface before they started propagating along this interface. To get more information about the condition in which a type (D) interface controls fault propagation, we set up several variations of these simulations. First we found that a high viscosity, high rigidity channel has a similar effect on fault propagation, inducing fault propagation along (though not within) the embedded channel of strong material. We also found that if the embedded fault and shear zone are significantly unaligned with the expected fault strike (~20° from the fault trace in the homogeneous simulation) and if they are embedded ~10 km from the fault nucleation site in its compressional quadrant, the fault propagation path would not be significantly affected by the material heterogeneity. Insignificant fault-interface interactions were also observed in a simulation without an embedded crustal fault zone and with a slightly more subtle viscosity contrast (10 19 and 1020 PaS, for the shear zone and mantle, respectively). However, other simulations with a mantle shear zone not overlain by a crustal fault zone, but where the viscosity contrast was more significant (1018 and 1023 PaS), the crustal fault formed directly above the mantle shear zone even though it was not aligned with the expected fault propagation path (forming a fault path identical to that in Fig. 4.4f).  103  We conclude that various types of material interface have a significant effect on the geometrical configuration of evolving strike slip fault systems and on the propagation path of individual faults. Particularly, simulated faults tend to propagate to and along contrasts in effective plate thickness, and where lithospheric-scale heterogeneities such as shear zones are embedded within the model crust and mantle. Mantle shear zones significantly affect the crustal deformation above and may induce fault propagation to and along their path even where they do not align with regional stress. However, our results also indicate that significant rigidity contrasts (up to ~20%) in the crust may not have a significant effect on fault propagation and fault-system evolution.  Figure 4.4 Damage maps of simulated faults along various types of material interfaces. (a, b) Homogeneous and crustal rigidity contrast model (type A) showing some limited effect of the material interface on fault geometry. (c) The currently active fault (higher damage level) runs along the material interface while the initial fault (to its left) is less active. (d) In simulations with material contrast in crust and mantle, faults that were nucleated in the stronger side propagate into the weaker side and align with the material interface. (e, f) Faults that nucleate away from material interfaces of type C and D propagate to and then along the interface. In all panels the dashed white line indicates the position of the interface (left side weak), and the width of the area shown is 50 km across and 70 km (top three panels) or 150 km (bottom panels) along strike. All panels show damage at 6-10 km depth.  104  105  4.3.3 Asymmetric deformation patterns and surface velocities In this subsection we examine strike-slip deformation patterns along material interfaces by calculating the strain rate and fault parallel velocity associated with faults that evolve along such interfaces. We describe our observations of strain localization along the various types of material interface, and quantify asymmetric distribution of deformation where it is observed. The interseismic strain rate and velocity cross-sections in Figure 4.5 display significant localization of deformation along (or adjacent to) all types of simulated material interfaces. Regardless of the type of simulation (including homogeneous and all variations of heterogeneous models), the highest strain rates are observed in the seismogenic crust (typically at depths of 3-20 km). These strain localizations represent faults in our numerical simulations. Comparing Figures 4.5a and 4.5b we see that localization is enhanced along the material interfaces in the crust (Fig. 4.5a) and mantle (Fig. 4.5b). Particularly in the upper mantle (20-30 km depth) and in the lower crust (14-20 km depth), deformation is significantly more localized in simulations with a material interface. An additional observation from simulations with interface types A, B and AB (Fig. 4.5ac) is that while crustal deformation localizes on the material interface (or immediately to the weak side of it), mantle deformation may localize up to 10 kilometers from the interface. In Figure 4.5c, where mantle strain localization is offset by ~8 km and the weak lower crust displays distributed flow (z=22-30 km), crustal deformation is relatively distributed and asymmetric, with significant surface strain rates (3*10-14 sec-1) extending much farther into the weak side (up to 9 km from the interface). In simulations with a contrast in crustal thickness (interface type C), highly localized deformation (high strain rates in Fig. 4.5d) occurs in the thin crust near the interface with the thicker crust. The thick crust displays very little strain confined to the surface sediments near the interface. Above and below the thin crust, strain-rate and interseismic velocity are highly asymmetric relative to the material interface (and the fault core). At the surface, shear is highly asymmetric with some low gradient deformation (strain rates > 3*10-14 sec-1) extending as far as 30 km from the interface on the weak side. At 15-30 km depth, where mantle and crust are juxtaposed, strain is localized within a ~10 km shear zone centered about 5 km from the interface. In simulations with an embedded fault and shear zone (type D, Fig. 4.5e), deformation 106  is localized in the crust and upper mantle, it slightly diverges near the surface and it is widely distributed at the deeper part of the simulated mantle (z=40-50 km). The deformation around this type of material interface is symmetric at all depths. To quantify surface velocity patterns we calculate the velocity symmetry ratio (defined below), we measure the offset between the location of maximum surface strain rate and the material interface, and we measure the distribution of surface velocity. As absolute surface velocity values directly depend on the actual boundary velocities (which may slightly differ from the far-field velocity applied in our boundary conditions, see Finzi et al., 2009), we refer to the velocities as percentages of the relative boundary velocity for each of the simulations. In this way, the area around a fault rupture which is bracketed by velocity contours of 12 mm/yr (in a simulation with boundary velocity of 14 mm/yr) is referred to as the shear zone in which 80-85% of the deformation occurs (see Fig. 4.6). In a similar way the shear zone in which ~55% of deformation occurs (assuming 14 mm/yr on the boundaries) is bracketed by the 8 mm/yr contours (Fig. 4.6). These two measures of the surface velocity distribution (shear zones encompassing ~55% and 85% of deformation) are used to determine the surface velocity symmetry ratio. The velocity symmetry ratio is calculated relative to the location of maximum (near-surface) strain rate, which best represents the location of the surface trace of a fault (Fig. 4.6). As with the damage symmetry ratio, a velocity ratio of 0.5 corresponds to symmetric distribution. The velocity patterns observed in simulations of the various material interface types are summarized in Table 4.3.  107  Figure 4.5 Strain rate and velocity cross-sections of simulated faults along material interfaces. (a) Type A model showing high strain rates along the material interface and symmetric deformation. (b) Type B model showing increased shear localization in the mantle (compare to a), and slight asymmetry particularly in the far field (upper mantle and near surface). (c) Type AB model showing asymmetric deformation patterns. (d) Type C model showing shear localization within the thin crust (near the interface) and wide, asymmetric distribution of deformation in the upper mantle and near the surface of the thin crust side. (e) Type D model showing localized symmetric deformation within the embedded fault and shear zone. In all panels the dashed white line outlines the structure simulated, and the width of the area shown is 50 km across and 50 km deep.  108  109  110  Figure 4.6 Determining velocity profile asymmetry R85 (or R55). Black curve represents a velocity profile across the fault, in plan view. The shaded region of width C indicates where 85% of the relative velocity of the model boundaries (Vo) is accommodated. This region is defined so half of the relative motion is accommodated on each side of the fault (dotted line). A (in the softer material) is greater than B (on the stiffer side), and the asymmetry ratio R85 is A/C. A similar calculation for a region accommodating 55% of the relative plate motion yields R55.  Symmetry ratio  Shear zone width (km)  (R55)  (R85)  (C55)  (C85)  Offset from interface (km)  A  0.55-0.6  0.5-0.55  3-4  6-8  0  B  0.5  0.6-0.7  5  8-10  0-5  A+B  0.55-0.65  0.6-0.7  6-7  8-11  0-7  C  0.5-0.6  0.65-0.75  5-8  10-13  2-8  D  0.5  0.5-0.55  4-6  9-11  0  Homogeneous  0.5  0.5  4-5  6-9  -  Interface type  Table 4.3 Fault-parallel surface deformation patterns in simulations with various bi-material interfaces. Velocity symmetry ratio, width of shear zone (encompassing ~55% and 85% of deformation), and offsets between the location of maximum (near-surface) strain rate and the material interface as measured in our simulations. In each simulation we average symmetry ratios measured at different locations along strike, and we present the range of average values for each type of simulation. 111  The surface velocity measurements reveal the different effects that crustal and mantle interfaces have on the deformation patterns in our simulations. Where faults follow a crustal material interface (types A and AB) the near field deformation is asymmetric and is typically centered along the interface (i.e., R55 is 0.55-0.6, and the offset is ~0 km). Along interfaces of type A, where the interface is only in the crust, the far field velocity pattern is more symmetric with R85=0.5-0.55. In simulations with a mantle interface (types AB, B) the far-field surface velocity field is significantly asymmetric (R85 ratios of 0.6-0.7). This is observed even in simulations where the near-field velocity is symmetric (in simulations of type B, without a crustal interface). Faults in simulations with a mantle interface (types B and AB), and with a thick crust (Moho at 30-35 km depth), often form parallel faults or elongated stepovers with >30 km of overlapping faults (these situations are not described in Table 4.3). In these simulations, the surface deformation is widely distributed (80% shear and damage levels above 0.4 within 1618 km wide shear zone), and is asymmetric (velocity and damage symmetry ratios of 0.6-0.8). Simulations of type C, which consist of a significant lateral contrast in the elastic plate thickness, display asymmetric velocity patterns (ratios of R55=0.5-0.6 and R85=0.65-0.75) with large offsets and a widely distributed velocity (Table 4.3, Fig. 4.5d). The side with the thicker crustal layer and greater effective elastic plate thickness in these simulations is not deforming significantly. Brittle deformation in these simulations is highly localized within the thin crust (alongside the interface with the thick crust). Below and above the thin crust, the low viscosity mantle and the weak (strain hardening) surface layer, both display widely distributed deformation (Fig. 4.5d). In simulations with embedded fault and shear zones (type D), the velocity distribution is typically symmetric with a relatively narrow shear zone and fault core within the embedded heterogeneity. While the width of shear zone in these simulations is narrower than that in the homogeneous simulations, still the minimum width in both model types is grid dependent and both may represent near-complete localization of deformation onto a fault plane.  112  4.4 Discussion and conclusions Rheological contrasts in the lithosphere have a significant effect on the formation and evolution of strike-slip faults. They tend to localize strain and damage on the weak side of the interface, and they may induce a somewhat asymmetric damage distribution (Figures 4.3 and 4.4).  4.4.1 Crustal interfaces, mantle interfaces and crust-mantle interactions While a significant rheology contrast in the crust may affect the near-field deformation patterns (with damage and strain more readily occurring on the weak side of the interface), a lower crust or upper mantle viscosity contrast may induce asymmetric deformation over a larger spatial scale (Fig. 4.5c). Mantle material interfaces in models of type B and AB induce the formation of shear zones which have a significant effect on the distribution of crustal deformation and surface velocity patterns. The mantle shear zones (defined here as the zone accommodating about 80% of the relative motion across the model) are typically 10-18 km wide and are centered in the lower viscosity material, about 5-10 km from the interface (Table 4.2). The mantle shear zone typically induces strain localization in the overlying crust, with some distributed strain in the crust above the low viscosity mantle (Figures 4.5b and 4.5c). Where the lower crust is weak, the mantle shear zone may induce a decoupling zone in the lower crust, characterized by a relatively low effective viscosity and widely distributed deformation (Fig. 4.7). The lower-crust decoupling zone, which is about 15 to 30 km wide in our models, induces a broad distribution of upper crust deformation manifested in the formation of a complex fault-system with parallel faults and/or a broad crustal shear zone (Fig. 4.7, also Fig. 4.3). This finding is consistent with results from previous viscoelastic models which did not incorporate damage and fault evolution in the upper crust (Roy, 1998). In simulations with a crustal thickness contrast (type C), a very strong viscosity contrast where the relatively stiff lower crust is juxtaposed against the low-viscosity (M2) mantle induces strain localization in the mantle under the thin crust (offset by approximately 5-10 km from the interface). Crustal strain is highly localized in the thin crust near the interface with the thick 113  crust. In these models we note significant distributed deformation in the surface sediment layer (material S in Table 4.1), above the thin plate (Fig. 4.5d). This may be due to the material interface between sediments on this side of the model and diabase (crustal material C in Table 4.1) on the side with a thick plate. The strength contrast between the crust and sediments and the transition from crustal strain weakening (enhanced by damage processes) to strain hardening in the sediments (devoid of damage in our simulations) may explain the distinct strain localization in seismogenic depth and widening of strain distribution near the surface.  114  Figure 4.7 Velocity (a) and effective viscosity (b) across a fault zone of type B. The contrast in effective viscosity induces the formation of a mantle shear zone (80% of shear concentrated in a zone less than 20 km wide). Above the mantle interface in the lower crust (24-32 km depth), the effective viscosity is relatively low (of the order of 2*1020 Pa s) and deformation is diffuse. The crust exhibits reduced effective viscosity along two parallel faults. The weak (low viscosity) lower crust appears to be locally de-coupled from the mantle above the interface. Both panels are 80 km across fault and 50 km in the vertical direction (major ticks shown every 10 km). 115  Our models of type D show that either strong or weak, finite-width heterogeneities localize strain, and fault zones form adjacent to such features. Our findings are consistent with earlier research based on 2D (plane stress and plane strain) models of heterogeneous, viscous lithosphere under compression (Tommasi et al., 1995). Though these models do not incorporate fault formation in brittle upper crust, they do indicate strain localization associated with both strong and weak heterogeneities. Our models and the results from Tommasi et al. (1995) and references therein underline the importance of inherited structure to lithosphere deformation patterns.  4.4.2 Interpretation of surface velocity data along material interfaces As the simulated surface deformation patterns may be asymmetric and offset from the material interface, they necessarily deviate from predictions of uniform dislocation models with a vertical strike-slip fault. Interpretation of such surface velocity patterns requires careful attention to lateral contrasts in the lithospheric structure. While various studies acknowledge the importance of incorporating material interfaces when modeling geodetic data (e.g., Schmalzle et al., 2006), other studies attribute the asymmetry solely to a crustal elasticity contrast (e.g., LePichon et al., 2003 and 2005; Jolivet et al., 2008). To address possible pitfalls associated with using elastic dislocation models to interpret GPS velocity profiles, we extract interseismic surface velocity profiles from our fault evolution simulations and fit these profiles with asymmetric, elastic dislocation models. Interseismic velocity profiles may be interpreted using elastic dislocation models with a range of crustal rigidity contrasts. The magnitude (and even the sense) of the contrast depends on the assumed location of the dislocation (and interface; Le Pichon et al., 2003). Figure 4.8a illustrates this with a comparison of symmetric and asymmetric elastic dislocation models of velocities across the NAFZ in the Marmara Sea. Both a uniform model (Meade et al., 2002) and a highly asymmetric model (Le Pichon et al., 2003) can fit the available GPS velocity data, but in the former case, the NAFZ must be 15 km to the south of what appears to be its true position (Le Pichon et al., 2003). Figure 4.8a shows a situation with sparse GPS velocity data in the near field (due in part to the presence of the Marmara Sea). Given the fault location, the elasticity contrast  116  is varied to fit the GPS site velocities relative to the fault (that is, the long-spatial scale asymmetry of the velocity profile). The tradeoff between elasticity asymmetry and fault position is also apparent in Figure 4.8b, which shows the interseismic velocity profile from model C and a few examples of dislocation models that can fit these velocities. In Figure 4.8b, the sense of strength asymmetry in the elastic dislocation model is actually opposite that prescribed in our fault evolution models (i.e., with the strong side to the left). If the elastic dislocation were forced to be at the position of the material contrast we assumed in model C, the elasticity contrast would be large, with the stronger material on the right. Model C and our other models produce highly localized, interseismic surface strain. Locking depths for elastic dislocation models tuned to fit interseismic velocities for laterally homogeneous test models are of the order of 4 km. This may be due in part to continuously occurring, small earthquakes in these simulations (that is, the fault is never fully locked). This is a feature of the fault evolution modeling method because the fault geometry and earthquake history are evolved, not prescribed. If the exact position of a crustal strike-slip fault is not known, but a GPS velocity profile is available, this could lead to a misinterpretation of the inferred fault position or the magnitude of the inferred elasticity contrast. In the most typical scenario, the interseismic velocity profile would be interpreted using a uniform elastic dislocation model, and the estimated fault position would be off by 4 km (Fig. 4.8b). Note that even if one does not attempt to fit the shape of the velocity profile in the near field this error will result: there is simply more strain to one side of the evolved fault than the other. Figure 4.8c shows how our modeled velocity profiles compare with the results of dislocation models when we require the elastic dislocation and the evolved fault to be in the same location. Figure 4.8c shows that realistic contrasts in viscous rheology due to juxtaposition of different rock types can result in strongly asymmetric velocity profiles. Models AB, B and C may yield interseismic velocity profiles which may only be fit to elastic dislocation models (with the dislocation located at the fault zone) if extreme contrasts in elasticity are assumed.  117  Figure 4.8 Modeled and GPS velocity profiles for the NAFZ. The top portion of each figure shows the fault location (prescribed in dislocation models, or evolved in our models). (a) Two elastic dislocation models for the North Anatolian Fault in the Marmara Sea. Slip rates are 23-24 mm/yr and locking depths are indicated. (b) Four examples of elastic dislocation models that produce asymmetric velocity profiles similar to our model C. The dislocation models require locking depths of about 4 km, probably because of ongoing small earthquakes over the interseismic interval in our models. (c) Interseismic velocity profiles from models AB, B, and C, and NAFZ GPS velocities from panel a. Our modeled velocities have been scaled to match a relative plate rate of 23 mm/yr. x is 0 at the fault trace. Model B evolved two major, parallel faults, so x = 0 at an average of their positions, and velocity profiles from epochs when each fault was dominant are shown. The fit to GPS data is improved by modeling a lower NAFZ slip rate (which is closer to the geologic rate).  118  119  It might be argued that some of the difference in inferred rigidity may be justified as resulting from the frequency dependence of the elastic moduli of rocks. Geodetic observations sample static (vanishing frequency) elastic properties and seismic observations sample dynamic (high-frequency) elastic properties. Published estimates of elastic parameters for various rock types are often based on seismic observations and experiments involving high-frequency energy transmission through small rock samples (e.g., Paterson and Wong, 2005). However, various compilations show that the static rigidity is consistently lower than the dynamic rigidity (e.g., Eissa and Kazi, 1988; Budiansky and O’Connell, 1976, Ciccotti and Mulargia, 2004; Ide, 1936), particularly for rocks near the Earth’s surface, which contain more microcracks. If material on one side of the fault is more damaged than that on the other side, large contrasts in elasticity may result (Ben Zion et al., 2003). However, as distributed fault damage is typically confined to the top few kilometers of the crust (Peng and Ben-Zion, 2004; Lewis et al., 2005; Dor et al., 2006, 2008; Finzi et al., 2009), extreme rigidity contrasts to great depths are not a viable explanation for asymmetric interseismic surface velocity data, and a contrast in mantle and/or lower-crust viscosity should be considered as a preferred explanation. Our results suggest that in locations where several parallel strike-slip faults are present, lateral variations in viscosity structure (or alternatively, elastic plate thickness) may result in incorrect inferences of slip rate for each fault. A similar point was made by Chery (2008) based on the observation that steady-state strain rate should correlate with plate thickness variations rather than individual fault slip rates. The Chery (2008) model relies on steady-state deformation and essentially no shear tractions on the base of the upper crust (extremely low effective viscosity for the mid-to lower crust). We reach a similar conclusion without relying on this assumption, which is inconsistent with most observations of postseismic deformation (e.g., summary in Burgmann and Dresen, 2008).  4.4.3 Propagation of faults to material interfaces The observed tendency of simulated faults to migrate toward material interfaces consisting of large scale contrasts in viscosity (and effective plate thickness) is consistent with various fault-system evolution studies. The evolving North Anatolian Fault (NAF) is suggested to have been attracted toward the margins of the Sea of Marmara, the north Aegean and the Saros 120  troughs which constitute a strong material heterogeneity (bi-material interface), where the basement is vertically offset several kilometers by the fault and the crust in the trough is highly sheared and faulted (Le-Pichon et al., 2005; Le Pichon et al., 2003). In addition studies by Armijo et al. (1996) and Sengor et al. (2005) suggest that the NAF propagation path has been affected by the material interface along the Black Sea (where continental and oceanic lithospheres are juxtaposed by ancient plate boundary structures), and by the pre-existing Aegean rift. Another plate boundary that is suggested to have evolved along a pre-existing material interface is the Dead Sea transform which propagated along a continental margin with significant lateral variations in effective plate thickness (Lyakhovsky et al., 1994). The above studies are consistent with our numerical results, showing that contrasts in effective plate thickness and preexisting weak zones may significantly influence the propagation path of evolving faults. Our result, suggesting that a crustal interface consisting of a rigidity contrast does not significantly attract propagating faults may appear to contradict previous numerical models of dynamic rupture propagation patterns at the vicinity of such material interfaces (Brietzke and Ben Zion, 2006). In these simulations Brietzke and Ben Zion show that ruptures are attracted to rigidity contrasts and pre-existing weak zones. The differences between static and dynamic strength could explain why rigidity contrasts attracted propagating ruptures in their work but had only a minor effect on fault propagation path in our (static) simulations.  4.4.4 Conclusions We have shown that asymmetric surface velocities around strike-slip faults may be due to contrasts in effective viscosity. Where elastic dislocation models require non-realistic rigidity ratios (>2) to explain observed asymmetric surface velocity profiles, this asymmetry is likely due to a lateral contrast in lower crust and/or upper mantle viscosity. Our results also show the great sensitivity of geodetic models to the location of the main active fault, and manifest the great care needed in interpretation of asymmetric velocity observations. Propagating faults are attracted to rheological contrasts, especially to channels of rheologically distinct material. This suggests that old lithosphere sutures may be more important to the propagation path of a new fault than a contrast in rheology between the two blocks of material that have been sutured together. Our models suggest that faults do not locate at the 121  material contrast, but that they form a few kilometers from the interface, in the softer material. This suggests that dislocation and earthquake-cycle models of fault zones in which the fault is co-located with the material contrast may yield inaccurate velocities in the near field. We also find that lateral contrasts in viscosity structure cause an asymmetric distribution of shallow damage around fault segments. This suggests caution in interpreting asymmetric damage distribution in terms of rupture propagation direction, though more work is required to determine whether our finding holds for the highly localized, extreme damage noted at smaller spatial scales which are interpreted as indicating directivity.  122  Bibliography Armijo, R., B. Meyer, G. C. P. King, A. Rigo, and D. Papanastassiou, (1996). Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean. Geophys. J. Int. 126, 11–53. Barbot, S., Fialko, Y., Sandwell, D., (2009), Three-dimensional models of elasto-static deformation in heterogeneous media, with application to the Eastern California Shear Zone. Submitted to Geophys. J. Int., 179, 500-520, DOI:  10.1111/j.1365-  246X.2009.04194.x Ben-Zion,Y.and D. J. Andrews, (1998), Properties and implications of dynamic rupture along a material interface, Bull. seism. Soc. Am., 88(4), 1085–1094. Ben-Zion, Y. and Sammis, C. (2003), Characterization of Fault Zones, Pure and Appl. Geophys. 160, 677–715. Ben-Zion, Y., Z. Peng, D. Okaya, L. Seeber, J. G. Armbruster, N. Ozer, A. J. Michael, S. Baris and Aktar, M. (2003), A shallow fault zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey, Geophys. J. Int. 152, 699-717. Ben-Zion, Y. and Lyakhovsky, V. (2006), Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology, Geophys. J. Int. 165, 197-210. Burgmann, R. and G. Dresen (2008), Rheology of the Lower Crust and Upper Mantle: Evidence from Rock Mechanics, Geodesy, and Field Observations, Ann. Rev. Earth Plan. Sci., 36, 531-567 Brietzke, G. B. and Y. Ben-Zion, 2006. Examining tendencies of in-plane rupture to migrate to material interfaces, Geophys. J. Int., 167, 807–819. Budiansky, B. and R. O’Connell (1976), Elastic Modulii of a cracked solid, Int. J. Solids Structures, 12, 81-97. Carter, N.L. and M.C. Tsenn (1987), Flow properties of continental lithosphere, Tectonophysics, 136, 27-63.  123  Chery, J. (2008), Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate), Earth and Planetary Science Letters, 269, 352-365. Ciccotti, F. and M. Mulargia (2004), Differences between static and dynamic elastic moduli of a typical seismogenic rock, Geophys. J. Int., 157, 474-477. Cochran, E.S., Li, Y.G., Shearer, P.M., Barbot, S., Fialko, Y., Vidale, J.E., (2009), Seismic and geodetic evidence for extensive, long-lived fault damage zones. Geology, v. 37, no. 4, p. 315-318, DOI: 10.1130/G25306A.1. Dor O., Ben-Zion, Y., Rockwell, T.K., Brune, J.N., (2006), Pulverized rocks in the Mojave section of the San Andreas Fault Zone. Earth Planet. Sci. Lett. 245, 642-654, DOI:10.1016/j.epsl.2006.03.034. Dor, O., Yildirim, C., Rockwell, T.K., Ben-Zion, Y., Emre, O., Sisk, M., Duman, T.Y., (2008), Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction. Geophys. J. Int.,  173, 483–504, doi: 10.1111/j.1365-  246X.2008.03709.x. Eissa, E. A. and A. Kazi (1988), Relation between static and dynamic Young’s Moduli of Rocks, Int. J. Rock Mech. Min Sci and Geomech. Abstr., 25 (6), 479-482. Fialko, Y. (2004), Probing the mechanical properties of seismically active crust with space geodesy: Study of the co-seismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake. J. Geophys. Res. 109, B03307, doi:10.1029/2003JB002756. Fialko, Y., Simons, M. and D. Agnew, (2001), The complete (3-D) surface displacement field in the epicentral area of the 1999 Mw 7.1 Hector Mine earthquake, southern California, from space geodetic observations. Geophys. Res. Lett. 28, 3063-3066. Fialko, Y., Sandwell, D., Simons, M. & Rosen, P. (2005), Three-dimensional deformation caused by the Bam, Iran, earthquake and the origin of shallow slip deficit. Nature 435, 295–-299.  124  Fialko, Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: Study of the co-seismic deformation due to the 1992 Mw7.3 Landers (southern California) earthquake, J. Geophys Res., 109, B03307, doi:10.1029/2003JB002756. Fialko, Y., (2006), Interseismic strain accumulation and the earthquake potential on the southern San Andreas fault system, Nature, 411, DOI: 10.1038/nature04797,968-971. Finzi, Y., Hearn, E. H., Y., Ben-Zion. And Lyakhovsky, V. (2009), Structural properties and deformation patterns of evolving strike-slip faults: Numerical simulations incorporating damage rheology, Pure Appl. Geophys., 166, 1537-1573, DOI: 10.1007/s00024-0090522-1. Karato, S.-I., M. S. Paterson, and J. D. FitzGerald (1986), Rheology of synthetic olivine aggregates: Influence of grain size and water, J. Geophys. Research, 91, 8151-8176. Kirby, S. H. and Kronenberg, A. K., (1987), Rheology of the lithosphere: selected topics. Rev. Geophys., 25: 1219-1244. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y. and Lockner, D. (2004), A Visco-elastic damage model with applications to stable and unstable fracturing, Geophys. J. Int. 159, 1155–1165. Hauksson, E. (2000), Crustal structure and seismicity distribution adjacent to the Pacific and North America plate boundary in southern California. J. Geophys. Res. 105, 1387513903. Ide, J. M. (1936), Comparison of statically and dynamically determined Young’s modulus of rocks, Proc. Nat. Acad. Sci., USA, 22, 81-92. Jolivet R., R. Cattin, N. Chamot-Rooke, C. Lasserre, and G. Peltzer, (2008), Thin-plate modeling of interseismic deformation and asymmetry across the Altyn Tagh fault zone, Geophys. Res. Lett., 35, L02309, DOI:10.1029/2007GL031511. Kirby, S. H. and Kronenberg, A. K. (1987), Rheology of the lithosphere: selected topics, Rev. Geophys. 25: 1219-1244.  125  Le Pichon, X., N. Chamot-Rooke, C. Rangin, and A. M. C. Sengor, (2003), The North Anatolian fault  in  the  Sea  of  Marmara,  J.  Geophys.  Res.,  108(B4),  2179,  doi:10.1029/2002JB001862. Le Pichon, X., C. Kreemer, and N. Chamot-Rooke (2005), Asymmetry in elastic properties and the evolution of large continental strike-slip faults, J. Geophys. Res., 110, B03405, doi:10.1029/2004JB003343. Lewis, M.A., Peng, Z., Ben-Zion, Y., Vernon, F., (2005), Shallow seismic trapping structure in the San Jacinto fault zone. Geophys. J. Int. 162, 867–881, doi:10.1111/j.1365246X.2005.02684.x, 2005. Li, V. C. & Rice, J. (1987), Crustal deformation in great California earthquake cycles. J. Geophys. Res. 92, 11533-11551. Lisowski, M., Savage, J. and W.H. Prescott (1991), The velocity field along the San Andreas fault in central and southern California, J. Geophys. Res. 96, 8369–-8389. Lyakhovsky, V., Z. Ben-Avraham, M. Achmon, (1994), The origin of the Dead Sea rift. Tectonophysics, 240: 29-43. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon, (1997a), Distributed damage, faulting, and friction, J. Geophys. Res., 102: 27,635-27,649. Lyakhovsky, V., Z. Reches, R. Weinberger, and T. E. Scott, (1997b), Non-linear elastic behavior of damaged rocks, Geophys. J. Int., 130: 157-166. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon, (2001), Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere, J. Geophys. Res., 106: 4103-4120. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon, (2005), A visco-elastic damage rheology and rateand-state-dependent friction. Geophys. J. Int.. 161, 179-190, doi:10.1111/j.1365246X.2005.02583.x. Meade B., B. H.Hager, S. C. McClusky, R. E. Reilinger, S. Ergintav, S. Lenk, A. Barka, and H. Ozener, ( 2002), Estimates of seismic potential in the Marmara Sea region from block models of secular deformation constrained by global positioning system measurements, Bull. Seismol. Soc. Am., 92, 208-215. 126  Molnar P. and C. H. Jones (2004), A test of laboratory based rheological parameters of olivine from an analysis of late Cenozoic convective removal of mantle lithosphere beneath the Sierra Nevada, California, USA. Geophys. J. Int., 156, 555-564. Paterson, M.S. and T.-F. Wong (2005), Experimental Rock Deformation - the brittle field, Second Edition, Springer Verlag, Berlin, 347 pp. Peng, Z. and Ben-Zion, Y. (2004), Systematic analysis of crustal anisotropy along the KaradereDuzce branch of the north Anatolian fault, Geophys. J. Int. 159, 253-274, doi:10.1111/j.1365- 46X.2004.02379.x. Ranalli, G. (1995), Rheology of the earth, 2nd edition, Chapman & Hall, London, 413 pp. Roy M. (1998), Evolution of fault systems at a strike-slip plate boundary: A viscoelastic model, Geophys. Res. Lett., 25, 15, 2881-2884. Savage J. C., W. Gan, W. H. Prescott, and J. L. Svarc, (2004), Strain accumulation across the Coast Ranges at the latitude of San Francisco, 1994-2000, J. Geophys. Res., 109, B03413, DOI: 10.1029/2003JB002612. Schmalzle, G., T. Dixon, R. Malservisi, and R. Govers, (2006), Strain accumulation across the Carrizo segment of the San Andreas Fault, California: Impact of laterally varying crustal properties, J. Geophys. Res., 111, B05403, doi:10.1029/2005JB003843. Sengor, A. M. C., Tuysz, O., Imren, C., Sakınc, M., Eyidogan, H., Gorur, N., Le Pichon, X. and Rangin, C. (2005), The North Anatolian Fault: A New Look, Annu. Rev. Earth Planet. Sci. 33: 37–112. Tommasi, A., A. Vauchez, and B. Daudré (1995), Initiation and propagation of shear zones in a heterogeneous continental lithosphere, J. Geophys. Res., 100(B11), 22,083–22,101. Wdowinski, S., B. Smith-Konter, Y. Bock, and D. Sandwell, 2007. Diffuse interseismic deformation across the Pacific–North America plate boundary. Geology, 35,4; 311–314; doi: 10.1130/G22938A.1  127  5  CONCLUDING CHAPTER  As a summary of my PhD research, the following sections present the key contributions of my thesis (Section 5.1), updates to the research described in my thesis chapters and recent advances in my field of study (Section 5.2), and a discussion of the potential limitations of my research and future research opportunities to further develop ideas presented in the thesis (Section 5.3).  5.1 Main contributions The topic of this thesis is the characterization of strike-slip fault structure and associated deformation patterns in the context of evolving fault systems. To numerically simulate faultsystem evolution I apply a damage rheology model that accounts for important aspects of fault evolution that are typically overlooked in numerical studies. In particular, damage rheology enables simulation of fault formation, propagation and long-term evolution of fault systems without pre-determination (numerical prescription) of fault planes, propagation paths and seismicity patterns. The damage rheology also accounts for off-fault deformation and fault zone material degradation and healing, and therefore provides a valuable tool for studying strain distribution and evolving deformation patterns within fault systems. To establish the applicability of the damage rheology model and numerical code to studies of fault evolution, several aspects of the code were benchmarked using analytic solutions and widely accepted numerical codes that do not apply damage rheology (Chapter 2, Appendix B, and Appendix C). An important research effort was dedicated to relating the numerical results to field observations and geophysical data in order to constrain model parameters based on natural fault processes (Chapters 2, 3). In particular, seismic and geodetic constraints on faultzone healing parameters improved the analysis of damage zone structure along fault segments and fault stepovers (Chapter 2), and further analysis of the healing parameter space revealed the important effects healing has on the structure of fault damage-zones, fault system evolution and on associated deformation patterns (Chapter 3).  128  Numerical simulations of evolving strike-slip faults were used to characterize damage and strain distribution across fault segments and within fault stepover zones (Chapter 2), fault system evolution patterns in a laterally homogeneous lithosphere (Chapters 2 and 3) and in the presence of lateral rheological contrasts (Chapter 4). They were also used to investigate how such lateral rheological contrasts affect surface velocity patterns, and to demonstrate the importance of accurate analysis of lithospheric structure and fault location in interpretation of geodetic data around strike-slip faults (Chapter 4). The main contributions of my thesis are described in the following sub-sections.  5.1.1 Establishing damage rheology applicability to studies of natural fault systems To establish the credibility of my numerical models and to encourage widespread incorporation of damage rheology in crustal deformation modeling I construct several suites of simplified models to test key numerical and physical (rheological) aspects of the numerical framework used in my work. Simulations without damage evolution are used to benchmark viscous deformation calculations, to establish the advantages of applying a variable force boundary condition, and to characterize a range of admissible stiffness parameters for these boundary conditions (Appendix B, Finzi et al., 2007). An analysis of the healing parameter space constrained by geophysical observations of shear modulus reduction in damaged fault zones is used to calibrate model parameters C1 and C2 (Chapters 2, 3). The constrained range of healing parameters (C1 = 10-24 to 10-10 s-1 and C2 = 0.015 to 0.035) is significantly reduced relative to previous studies (e.g., Lyakhovsky et al., 2005). As the healing parameters are important factors in the simulation of damage-zone mechanical properties and fault system evolution patterns their calibration greatly improves the applicability of damage rheology models to studies of natural fault systems. Finally, I demonstrate the applicability of my models by documenting features in the model results that are characteristic of natural strike-slip fault systems. These include flower structures, branching faults, formation of fault stepovers, and fault evolution patterns such as segment linkage and strain localization.  129  5.1.2 Damage-zone structure and deformation patterns along segmented strike-slip faults Fault segments in our models are surrounded by damage zones, which develop early in the evolution of the fault system (before about 50 meters of total offset has accumulated). Each damage zone resembles a flower structure in which damage is highly localized at depth but widely distributed in the top 3-10 kilometers of the crust. The formation of this flower structure is enhanced by rapid and near-complete postseismic healing of damage at depths exceeding 5-10 km. Simulations also show that off-fault strain distribution correlates with the distribution of permanent shallow damage along fault segments. The fault stepover zones in our models exhibit extensive damage and elasticity degradation sustained during many earthquake cycles. Such damage patterns, associated with intense tensile fracturing and dilation, are consistent with recent structural evolution models for dilational stepovers (De Paola et al., 2007), and with mineral exploration studies that relate hydrothermal ore deposits to long-lasting permeability increase within fault stepovers (Micklethwaite and Cox, 2004; Sheldon and Micklethwaite, 2007). An important implication of these extensive damage zones concerns the reduction of interseismic stress accumulation in the weakened stepover region, which could aid earthquake rupture arrest. We also note that the extensive rock damage near stepovers could amplify seismic ground motion, and should be detectable with focused seismic and geodetic studies. We also find that lateral rheological contrasts in the lithosphere (particularly in the viscosity structure) have a significant effect on the formation and evolution of damage zones. Strain and damage localize along the weak side of the interface, and an asymmetric damage distribution relative to the position of the fault may result. a significant rheology contrast in the crust affects near-field deformation patterns, with more damage and strain on the weak side of the interface. A viscosity contrast deeper in the crust or in the upper mantle may induce asymmetric crustal deformation over a broad spatial scale, but little asymmetry of the shallow damage zone (in the near field). In the upper mantle and lower crust, such viscosity contrasts may induce the formation of a shear zone extending as much as 10 kilometers from the interface (into the low viscosity side).  130  5.1.3 Revealing the important role healing efficiency plays in fault system evolution Recent geophysical fault zone observations suggest that interseismic (and postseismic) fault zone damage levels vary significantly within fault systems (e.g., Cochran et al., 2009; Barbot et al., 2009). Since the characteristic interseismic damage levels (and hence healing effectiveness) vary significantly within natural fault settings, models with combinations of damage model parameters which yield a wide range of interseismic damage levels are geologically admissible and should be investigated. We note that in past studies including ours (Finzi et al., 2006, 2007), model healing parameters were selected to allow a range of healing time scales constrained by estimations of early postseismic healing rates (Lyakhovsky et al., 2005; Lyakhovsky and Ben-Zion, 2008), but they did not represent a wide range of interseismic damage levels. This resulted in a partial representation of the healing parameter space and incomplete analysis of how healing affects fault zone evolution. Fault evolution simulations applying the full range of healing parameters (as determined in Chapters 2 and further constrained in Chapter 3) show that healing effectiveness controls both the dimensions of the fault damage zone and the complexity of the fault system (Chapter 3, Finzi et al., 2008). Conditions favoring complete healing yield wider distributed off-fault (DOF) damage in the top few kilometers of the crust and shallower fault-core damage, whereas conditions favoring ineffective (limited) healing yield narrower DOF damage and deeper faultcore damage zones. In addition, healing effectiveness is an important factor in the temporal stability of fault complexities (e.g., stepovers), with shorter-lived complexities expected where conditions favor more complete healing. Healing effectiveness also has an important effect on the structural properties of fault stepovers, with deeper more intensive damage expected where conditions favor partial (less effective) healing. Based on this, we further suggest that large longlived stepover zones represent unfavorable conditions for healing that should also correlate with significant material weakening down to the lower parts of the seismogenic zone. The above findings suggest that observable structural characteristics of a fault system and long-term evolution patterns may serve as first order indicators of gouge-scale healing processes. Establishing such a link would yield a powerful tool to study fault systems where detailed observations of fault zone properties are not available. Finally, as an important factor in the  131  evolution of fault complexity and fault structure (particularly at stepovers), healing effectiveness may have important implications for earthquake propagation and seismic hazard studies. . 5.1.4 Asymmetric surface deformation due to lateral contrasts in rheology, and interpretation of geodetic data We have shown that asymmetric surface velocities around strike-slip faults may be due to contrasts in effective viscosity. Elastic dislocation models require unrealistic rigidity ratios (>2) to explain observed asymmetric surface velocity profiles, but this asymmetry is more likely due to a lateral contrast in lower crust and/or upper mantle viscosity. Our results also show the great sensitivity of geodetic models to the location of the main active fault, and demonstrate that great care is needed in the interpretation of asymmetric velocity observations. If the exact position of a crustal strike-slip fault is not known, but a GPS velocity profile is available, lateral contrasts in viscosity stricture could lead to a misinterpretation of the inferred fault position or the magnitude of the inferred elasticity contrast. In regions with multiple fault strands, this could lead to an incorrect assignment of fault slip rates, which is of obvious importance to seismic hazard assessments.  5.2 Updates on working hypotheses and on current knowledge in the field As stated above, the healing parameter study described in detail in chapter 3 is an improvement of the analysis initially presented in chapter 2. While the results and conclusions given in chapter 2 remain valid to date, the improved healing parameter analysis in chapter 3 indicates that the effect of healing on damage zone structure and fault system evolution was greatly underestimated in chapter 2. The research described in chapter 3 was inspired by new observations of postseismic and interseismic fault zone strength within the Eastern California Shear Zone (e.g., Cochran et al., 2009; Barbot et al., 2009; Hearn and Fialko, 2009), which made us realize that healing effectiveness may vary significantly, and these variations should be represented in our model parameters within the existing parameter space. Characterizing the healing parameter space in terms of the characteristic interseismic damage level (ch) and the healing time scale enabled us to reproduce the observed variation in interseismic damage levels and to better understand how healing affects damage zone structure and fault system evolution. 132  Various studies in the past couple of years have focused on the topics presented in this thesis, and support our findings or complement them by providing new insights on aspects related to my thesis. The following is a review of several such studies. Using a damage rheology model and numerical simulations similar to those used in my thesis, Lyakhovsky and Ben-Zion (2009) compare fault evolution on a scale of a single crack and on the scale of lithospheric deformation. They associate the formation of initially segmented (enechelon) strike-slip faults with the asymmetric generation of new damage at the tip of existing fault segments (correlative to the formation of wing cracks at the tip of a small-scale crack under mixed-mode loading). Lyakhovsky and Ben-Zion (2009) then show how initially segmented fault systems evolve with continuing deformation to a through-going localized structure, in accordance with my results on fault system evolution (Chapters 2 and 3, Finzi et al., 2007). However, their analysis of how healing affects the evolution of a segmented fault is done in terms of the ratio between healing and loading time scales. It does not consider the entire admissible range of variations in healing effectiveness because healing parameters are co-varied in a way that represents a narrow range of healing effectiveness (and interseismic damage) levels. This leads to their interpretation that slower healing (relative to the loading rate) yields more efficient fault zone localization and reduction of geometrical complexity. This is counter intuitive as slow healing would result in long-lasting fault segments (and stepovers), and while the initial en-echelon fault configuration fails to heal and become inactive, the localization of deformation on a through-going fault and the smoothing of the fault system would be delayed or prevented. Chapter 3 of my thesis presents a thorough analysis of the healing parameters and a more complete numerical parameter space study. In this chapter I conclude that within the parameter space representative of natural settings, effective healing results in shorter-lived fault complexities (more efficient fault zone localization and reduction of geometrical complexity) and ineffective healing results in longer lasting fault complexity (slower and inefficient reduction of complexity), regardless of the healing time scale. Finally, in a simulation consisting of lithospheric-scale heterogeneity in the viscosity structure, Lyakhovsky and Ben-Zion (2009) show the formation of shear zones along the material interface in the lower crust and mantle. The material interface at the top of the crust (where viscosity is high and therefore viscosity contrasts have a lesser effect on fault evolution) 133  induces long-lived distributed damage along the interface exhibiting asymmetric patterns with most of the damage occurring within the weaker sedimentary layer (sediment rigidity is about half of the crustal rock rigidity in their model). These results, which are comparable to those presented in chapter 4, were achieved using a significantly different model of lithospheric structure and material interface and different boundary conditions. The distribution of damage and strain within the sediment layer in the Lyakhovsky and Ben-Zion (2009) model is much broader than that observed in our simulations with extreme crustal rigidity contrasts. This difference may be attributed to the different lithospheric structures simulated, as their interface (between sediment and crustal-rock) is a shallow feature (3-8 km deep) and wider damage distributions are expected at that depth (compared to the highly localized damage expected at mid seismogenic depth).  Recent observational studies address the structural evolution of fault systems and try to quantify various aspects of this process. Ferrill et al. (2008) use photogeologic maps, InSAR data and various time markers in recent alluvial deposits around faults ruptured in the 1992 Landers, California earthquake to infer that damage zone formation occurs during the very early stage of fault evolution (i.e. during the accumulation of several centimeters of fault offset) followed by strain localization and reduction of the geometric complexity of the active fault. This time frame for damage build-up is an order of magnitude shorter than that suggested in chapter 2 of my thesis (i.e. the mapped damage zones reached their maximal extent (width) much faster than faults in my simulations). The discrepancy in time scales may be attributed to the order of magnitude difference between observed and simulated damage zone dimensions. As the simulated faults and damage zones are much larger it seems reasonable that the time frame for damage zone formation would be longer. De Joussineau and Aydin (2009) analyze geometrical features of strike-slip faults in the Valley of Fire State Park, Nevada and conclude that they display a negative power law relation between the number of fault steps and the cumulative slip, and a positive power law relation between cumulative slip and both the mean segment length and the mean step dimensions. In agreement with previous studies (see review by Ben-Zion and Sammis, 2003) and with my results (Chapter 2), their analysis of fault trace complexity indicates overall coalescence and 134  smoothing with progressive fault growth. The formulation of quantitative relations between complexity indicators and total fault offset is an important contribution that may be used to determine stages in fault system evolution and to evaluate the maturity of fault systems. Such quantitative relations may also be used to calibrate numerical models of fault evolution (such as the one used in my research) with observations from natural fault systems.  5.3 Limitations of reported research and outline of potential future research 5.3.1 Loose constraints on material parameters and sparse observation of in situ damage evolution As described in chapter 3, a significant weakness of the damage rheology model is the numerous model parameters, which are sometimes difficult to relate to observable phenomena. In recent years, much effort has gone into relating these parameters to results of brittle deformation experiments and thermodynamic theory (e.g., Lyakhovsky et al., 2001, 2005; Hamiel et al., 2004, 2006; Lyakhovsky and Ben-Zion, 2008). Damage model parameters are fairly well constrained, however many of these constraints are based on laboratory experiments (e.g., Hamiel et al., 2009) and on simplified numerical tests that reproduce empirical or theoretical scaling relations (e.g., Lyakhovsky et al., 2005; Ben-Zion and Lyakhovsky, 2006). The applicability of these parameters to studies of lithospheric deformation processes may not be taken for granted and additional work is needed to calibrate damage model parameters based on fault zone observations. Unfortunately, quantitative descriptions of fault zone material (mechanical) properties are sparse and typically limited to the top kilometer of the crust and to inactive exhumed faults. To address this problem we synthesized a large volume of recent seismic and geodetic measurements from various strike-slip faults to better constrain the healing rate parameters, C1 and C2 (Chapters 2 and 3). In our analysis of the healing rate we also discuss the shortage of continuous or repeating measurements of fault zone properties that could potentially detect damage evolution (healing) occurring after the very early postseismic stage. A potential research project that could greatly improve our knowledge of fault damage evolution patterns in natural settings, concerns the damage zones at fault stepovers. Our numerical simulations predict that fault stepover zones exhibit extensive damage (consistent with intense fracturing and dilation) and elasticity degradation sustained during many earthquake 135  cycles. The stepover zones are therefore expected to exhibit local modifications of surface deformation patterns that should be detectable with detailed seismic and geodetic imaging (e.g., Hamiel and Fialko 2007, Hearn and Fialko 2009, Cochran et al., 2009). The internal structure and mechanical properties within the stepover zones could be determined by analysis of geodetic (InSAR) data combined with field mapping and geophysical fault zone surveys (e.g., tomography or analysis of trapped waves). Surface deformation patterns from simulations tailored to represent the structure and properties of fault stepovers could then be compared with those measured in the field in order to calibrate model parameters (e.g., healing rate parameters and the damage related viscosity Cv). The calibrated models could then be used to study various aspects of fault system deformation and evolution that concern fault stepover. Specifically, longterm simulations could shed light on the conditions in which such stepover zones are expected to arrest propagating earthquakes of various magnitudes.  5.3.2 Computational limitations As the damage rheology model we employ describes the evolving properties of crack distribution rather than those of individual cracks (Lyakhovsky and Myasnikov, 1984 and 1985), it is applicable to volumes with a sufficiently large number of relatively small cracks. Furthermore, our numerical models are inherently discrete with an abrupt transition from initial (static) to final (dynamic) stress levels during brittle failure (e.g., Rice and Ben-Zion, 1996). For these reasons, several aspects of our results are expected to be grid-dependent. For example, minimum stress drop, earthquake magnitude and damage zone dimensions depend on grid characteristics. An analysis of our entire set of damage zone simulations indicates that while the spatial extent of damage is somewhat grid -dependent, the average level of damage within voluminous damage zones (e.g., stepover zones) is probably not grid-size dependent. The analysis described in chapter 2 also indicates that the spatial extent of damage zones in our models with element dimensions between 0.6 and 4 km are not significantly sensitive to element size. Yet, localized damage such as that typical of the fault core at depth is portrayed in our simulations as a deformation band one element wide (i.e. our model cannot accurately represent strain localization on fault planes and shear zones narrower that our element widths). In our  136  simulations, discontinuous damage distributed on a well defined plane with maximum width of one model element is interpreted as an indication of highly localized strain along a fault plane. Another weakness of our numerical models concerns the computational load resulting from the wide range of spatial and temporal scales in our simulations. A lithosphere-scale simulation with model elements smaller that 0.5 km and a basic (minimum) time step of a fraction of a second run very slowly and become numerically unstable when large total strains accumulate. These simulations typically terminate (crash) after 10-20 days of CPU time, during which the simulated damage zones do not reach a stable width (the simulations run for 15001800 model years and reach a total slip of approximately 0.05-0.06 km). Because of these problems, our models cannot conclusively predict the details of surface damage patterns, the width of the fault core, or the geometry of small faults and fractures within flower structures and stepover regions. Significant computational improvements such as optimizing the code for parallel computing would greatly enhance the applicability of the numerical model to studies of faultsystem evolution. Such improvements would enable more detailed simulations of fault structures and damage patterns at the surface, and would also allow the incorporation of additional deformation mechanisms and coupling of damage evolution with fluid flow. One potential research project that would become plausible is the quantitative analysis of evolving fault zone complexity and the relation between cumulative slip and geometrical fault-trace characteristics (to be compared with field observations of De Joussineau and Aydin, 2009). While the current framework is able to simulate the evolution of a strike-slip fault consisting of 3-5 segments in a time frame of about 10,000 years, an improved code would enable simulations of larger model domains with more detailed structures (more segments and clearer distinction from grid induced kinks in fault segments) and longer simulated time frame which could result in a statistic description of the evolution process. Such numerical simulations calibrated by field observations would enable a study of how lithospheric conditions (e.g., strain rate and heat flow) and rheological properties (e.g., healing effectiveness) affect the structural evolution and deformation patterns of strike-slip fault systems. Another important research project that would become plausible after the numerical code is optimized concerns the evolution of seismicity patterns within and around fault stepovers. Our 137  current research indicates that fault stepovers experience ongoing seismicity during the entire earthquake cycle (Chapter 2). Detailed simulations of segmented strike-slip faults could characterize the spatial and temporal evolution of seismicity and test whether seismicity patterns in stepover zones change during various stages of the earthquake cycle on the main fault segments. In addition, these simulations could be used to study the conditions in which fault stepovers arrest rupture propagation, though better representation of dynamic rupture processes in the code would be necessary for a thorough study of the interaction between earthquake ruptures and fault stepovers. Yet another important research direction that would become plausible after numerical optimization is the simulation of coupled fluid-flow and deformation processes along natural faults. The interaction between fracture and fluid flow plays an important role in mechanical and hydrological crustal processes such as earthquakes, faulting and mineralization. Evolving permeability and pore pressure are also expected to play an important role in damage zone structural evolution and fault-system evolution processes. Such simulations would apply the analytic model of Hamiel et al., 2005 which incorporates damage evolution and classical poroelastic formulation to evaluate porosity and permeability evolution during deformation.  .  138  Bibliography Barbot, S., Fialko, Y., Sandwell, D., 2009, Three-dimensional models of elasto-static deformation in heterogeneous media, with application to the Eastern California Shear Zone, Geophys. J. Int., 179, 500-520, DOI: 10.1111/j.1365-246X.2009.04194.x Ben-Zion, Y., Sammis, C., 2003. Characterization of fault zones. Pure and Appl. Geophys. 160, 677–715. Ben-Zion, Y., Lyakhovsky, V., 2006. Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology. Geophys. J. Int. 165, 197-210. Cochran, E.S., Li, Y.G., Shearer, P.M., Barbot, S., Fialko, Y., Vidale, J.E., 2009. Seismic and geodetic evidence for extensive, long-lived fault damage zones. Geology, v. 37, no. 4, p. 315-318, DOI: 10.1130/G25306A.1. de Joussineau G., and A. Aydin, 2009. Segmentation along strike-slip faults revisited, Pure and Appl. Geophys., 166, 1575-1594, doi: 10.2007/s00024-009-0511-4. De Paola, N., Holdsworth, R.E., Collettini, C., McCaffrey, K.J.W., Barchi, M.R., 2007. The structural evolution of dilational step-overs in regional transtensional zones, in: Cunningham W.D., Mann, P. (Eds), Tectonics of Strike-Slip Restraining and Releasing Bends. Geological Society, London. Special Publication., 290, pp. 433-445. Ferril, D.A., K.J. Smart, and M. Necsoiu, 2008, Displacement-length scaling for single-event fault ruptures: insights from Newberry Springs Fault Zone and implications for fault zone structure, in: Wibberley C.A.J., W. Kurz, J. Imber, R.E. Holdsworth, and C. Collettini (eds.) The Internal Structure of Fault Zones: Implications for Mechanical and Fluid-Flow Properties, 299, 113-122. DOI: 10.1144/SP299.7. Finzi Y., Hearn E., Lyakhovsky V. and Ben-Zion Y., 2008, Fault zone healing and structural evolution of fault systems - observations and numerical simulations. Eos Trans. AGU, 89(53), Fall Meet. Suppl., Abstract T43E-08. Finzi Y., Hearn E., Lyakhovsky V. and Ben-Zion Y., 2007, Damage zone structure and deformation patterns along segmented strike slip faults. Southern California Earthquake Center Annual Meeting Procedures and Abstracts Vol. XVII, Palm Springs, California. Finzi Y., Hearn E., Lyakhovsky V. and Ben-Zion Y., 2006, 3D viscoelastic damage rheology models of strike-slip fault systems and their associated surface deformation. Eos Trans. AGU, 87(52), Fall Meet. Suppl., Abstract T21C-0425. Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y., Lockner, D., 2004. A Visco-elastic damage model with applications to stable and unstable fracturing. Geophys. J. Int. 159, 1155– 1165. Hamiel, Y., Lyakhovsky, V. and Agnon, A. (2005), Rock dilation, nonlinear deformation, and pore pressure change under shear, Earth Planet. Sci. Lett. 237, 577-589.  139  Hamiel, Y., Katz, O., Lyakhovsky., Reches, Z. and Fialko, Y. (2006), Stable and unstable damage evolution in rocks with implications to fracturing of granite, Geophys. J. Int. 167, 1005-1016. Hamiel, Y., Fialko, Y., 2007. Structure and mechanical properties of faults in the North Anatolian Fault system from InSAR observations of coseismic deformation due to the 1999 Izmit (Turkey) earthquake. J. Geophys Res. 112,B07412, doi:10.1029/2006JB004777. Hamiel, Y., V. Lyakhovsky, S. Stanchits, G. Dresen, and Y. Ben-Zion, 2009, Brittle deformation and damage-induced seismic wave anisotropy in rocks, Geophys. J. Int., 178,901-909, DOI: 10.1111/j.1365-246X.2009.04200.x. Hearn, E.H., Fialko, Y., 2009, Can compliant fault zones be used to measure absolute stresses in the upper crust?, J. Geophys. Res., 114, B04403,DOI: 10.1029/2008JB005901. Lyakhovsky, V. and Myasnikov, V. P., 1984, On the behavior of elastic cracked solid, Phys. Solid Earth, 10, 71-75. Lyakhovsky, V. and Myasnikov, V. P., 1985, On the behavior of visco-elastic cracked solid, Phys. Solid Earth, 4, 28-35. Lyakhovsky, V. Ben-Zion, Y., Agnon, A., 2001. Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere. J. Geophys. Res. 106: 4103-4120. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon, 2005, A visco-elastic damage rheology and rateand-state-dependent friction. Geophys. J. Int.. 161, 179-190, doi:10.1111/j.1365246X.2005.02583.x. Lyakhovsky, V., Ben-Zion, Y., 2008. Scaling relations of earthquakes and aseismic deformation in a damage rheology model. Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365246X.2007.03652.x. Lyakhovsky V., Ben-Zion, Y., 2009. Evolving fault zone structures in a damage rheology model. in review, Geochem., Geophys., Geosyst. Micklethwaite, S., Cox, S.F., 2004. Fault-segment rupture, aftershock-zone fluid flow and mineralization. Geology, 32, 813–816. Rice, J.R. and Ben-Zion, Y., 1996, Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci. USA 93: 3811-1818. Sheldon H. A., Micklethwaite S., 2007. Damage and permeability around faults: Implications for mineralization. Geology: Vol. 35, No. 10 pp. 903–906.  140  APPENDIX A  NUMERICAL PROCEDURES  The numerical code we use for our modeling utilizes the Fast Lagrangian Analysis of Continua (FLAC) method (e.g., Cundall and Board, 1988; Poliakov et al., 1993). This formulation is explicit in time and it continuously updates the shape functions of tetrahedral elements allowing large deformations to be simulated. The general procedure involves solving the equations of motion to determine nodal velocities and to calculate the elastic component of the element strains. The total strain tensor is calculated by summing the elastic strain component, a viscous strain component and a component representing the damage-related inelastic strain ijv. Hamiel et al. (2004) introduced the gradual inelastic strain, ijv, whose accumulation rate is proportional to the rate of damage accumulation:  d C v dt  ij  dt  0   d iji  d 0 dt d 0 dt  (A1)  Cv is a material constant and ij = ij – kkij /3 is the deviatoric stress tensor. The compliance, or inverse of viscosity (Cv∙d/dt), relates the deviatoric stress to the rate of irreversible strain accumulation. Ben-Zion and Lyakhovsky (2006) connected the rate of irreversible strain accumulation with partitioning between seismic and aseismic deformation in the seismogenic zone, and showed that the fraction of elastic strain released seismically, referred to as the seismic coupling coefficient , can be estimated as:    1 1 R  (A2)  The damage rheology constitutive and kinetic relations provide element stresses and local values of damage, which are used to update material properties (moduli  and ) and to calculate the nodal force balances applied in the next time step. This procedure continues through many time steps until damage in any element reaches the critical level =cr associated with brittle instability (Lyakhovsky and Ben-Zion, 2008). In our numerical procedure we assume that with the onset of instability the damage variable increases instantaneously to =1, leading to dynamic 141  rupture. The dynamic stress drop associated with brittle instability produces a reduction of the elastic strain in a rapid process, and abrupt increase of plastic deformation (corresponding to slip in simpler models with planar faults) is calculated following the associated plasticity law of Drucker (1949). The dynamic failure is arrested and post-failure material healing starts when the strain invariant ratio drops below a value d (d ≤ 0) and the system regains its stability. This is equivalent to setting the stress drop in frictional frameworks according to the difference between static and dynamic friction. The real dynamic process, wherein waves are generated by the stress drop, is not simulated here. However, a quasi-dynamic procedure is applied to simulate a rupture front propagation. This procedure is based on the fact that elastic waves propagating in a solid with increasing damage are amplified when the damage is close to the critical level defined by the convexity conditions (Ben-Zion & Lyakhovsky, 2006). While the initiation of brittle instability leading to macroscopic failure is at a critical level c of the damage state variable, the critical level of damage for propagation of the instability is lower than the initiation level by a dynamic weakening factor (Ben-Zion & Lyakhovsky, 2006). This is accomplished by recalculating the stress field after each stress drop in every element involved in the rupture process and reducing the critical value of the damage parameter, c, to d: 𝛼𝑑 = 𝛼𝑐 − 𝜏𝑟  𝑑𝛼 𝑑𝑡  (A3)  where r is a material parameter indicating the characteristic timescale for seismic wave damping in a damage free material. The right term in equation (A3) may be considered as a dynamic weakening of the material during the occurrence of brittle instability.  142  Bibliography Ben-Zion, Y. and V. Lyakhovsky, 2006. Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology, Geophys. J. Int., 165, 197- 210. Cundall, P. A. and M. Board, 1988. A microcomputer program for modeling large-strain plasticity problems, in: Numerical methods in geomechanics, Proc. 6th Int. Conf. Numerical Methods in Geomechanics, Innsbruck, ed. Swoboda, C., Rotterdam, Balkema, pp. 2101-2108. Drucker, D. C., 1949. Some implications of work-hardening and ideal plasticity. Quart. Appl. Math. 7, 411-418. Hamiel, Y., Y. Liu, V. Lyakhovsky, Y. Ben-Zion, and D. Lockner, 2004. A visco-elastic damage model with applications to stable and unstable fracturing. J. Geophys. Int. 159: 1-11. Lyakhovsky, V., and Y. Ben-Zion, 2008. Scaling relations of earthquakes and aseismic deformation in a damage rheology model, Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365-246X.2007.03652.x. Poliakov, A., P. Cundall, Y. Podladchikov, and V. Lyakhovsky, 1993. An explicit inertial method for the simulation of viscoelastic flow: an evaluation of elastic effects on diapiric flow in two- and three-layers model, In Proc. NATO Advanced Study Institute on Dynamic Modeling and Flow in the Earth and Planets, Runcorn, K.E. & Stone, D. (eds) Dordrecht, Kluwer, pp. 175-195.  143  APPENDIX B  BOUNDARY CONDITIONS  Numerical models of evolving fault systems typically incorporate either constant velocity or constant stress boundary conditions. Lyakhovsky and Ben-Zion (2008) show that these boundary conditions may generate very different velocity distributions during the interseismic period and conclude that both conditions are inappropriate for simulations of long-term fault system evolution. To perform such simulations, Lyakhovsky and Ben-Zion (2008) introduced a modified boundary condition in which forces are proportional to a stiffness of virtual springs multiplied by the mismatch (slip-deficit) between the far field plate motion and displacement of the boundary nodes (Fig. B.1). Unlike the constant-velocity or constant-stress boundary conditions, this boundary condition accounts for stress accumulation and strain-rate decrease during the inter-seismic periods, as well as abrupt stress reduction at the model boundary during seismic events. The variable force condition realistically adjusts the forces applied on the model domain boundaries according to the evolution of elastic properties, stress state, and seismic events within the simulated domain (Lyakhovsky and Ben-Zion, 2008). We constructed a series of test models to demonstrate the performance of this boundary condition and to study the spring stiffness parameter. Using these models we compared postseismic deformation in large domain models with fixed fault parallel boundaries, to results obtained from narrow-domain models with the variable force boundary condition. In the large models the domain width was set to be twenty times larger than the prescribed rupture depth (so that little deformation is expected to occur near the fault parallel boundaries), and the faultparallel boundaries were fixed. In the narrow-domain models, the fault perpendicular width was 140 km (i.e., fault parallel boundaries were 70 km from the fault). In both the narrow- and widedomain models the rheology was set to be viscoelastic without damage evolution, and with realistic material properties (=30 GPa, =0.25, elastic plate thickness 20 km, and Newtonian substrate viscosity of 1019 Pa s). The analysis of stiffness values and their effect on boundary forces and velocities shows that using a very large spring stiffness (>5 ∙106 MPa) results in a constant velocity on the boundary nodes. Free boundary conditions are approximated with the use of very low (<5∙102 MPa) stiffness values. Using intermediate stiffness values, the force applied on boundary nodes varies as a response to deformation processes within the model 144  domain, except for rare instances in which a constant force is applied on boundary nodes that experience a very large mismatch between the plate motion and the boundary displacement. To best represent constant velocities far from the fault zone (i.e., at the boundaries of the very large models) we suggest that the spring stiffness in narrow-domain models should be set to 5∙103 to 104 MPa. This is similar to the value used in the models presented in chapter 2 (104 MPa).  Figure B.1 Schematic illustration of the variable force boundary condition. The model domain is coupled to a constant, fault-parallel velocity in the far-field, with the degree of coupling depending on a stiffness parameter (springs). Red lines indicate shear localization imposed by the boundary conditions as they were applied in our models (localized shear is also applied on the bottom boundary of our models, not illustrated here). Tectonic velocities and domain dimensions are for illustration and may be varied.  145  Bibliography Lyakhovsky, V., and Y. Ben-Zion, 2008. Scaling relations of earthquakes and aseismic deformation in a damage rheology model, Geophys. J. Int. 172, 651-662, doi: 10.1111/j.1365246X.2007.03652.x.  146  APPENDIX C  CODE VERIFICATION: VISCOELASTIC COMPONENT  The rheology of the continental crust and upper mantle is typically described by empirical constitutive relations that express the strain rate as a function of stress and temperature for the relevant lithologies (e.g., Goetze, 1978; Brace and Kohlstedt, 1980; Kirby, 1983). The surface sediments, crust and upper mantle layers of our numerical simulations are all modeled as viscoelastic materials, obeying the following constitutive equation:  ε = A σn exp −  Q+PV ∗ RT  (C1)  where A is an empirical constant,  is stress, Q is activation energy, V* is activation volume, P is pressure, T is temperature, and R is the gas constant (R=8.3144 J/mol oK). While steady-state plastic flow of rocks can take place by a variety of mechanisms, direct observations of ductile shear zones and indirect analyses of seismic anisotropy in mantle rocks demonstrate that dislocation creep dominates at stresses above 10 MPa. For this type of flow mechanism the stress exponent n is typically between n=3 and n=5 (e.g., Karato, 1989; Weertman, 1978). The material parameters in our models are based either on common lithologies that are assumed to control continental ductile deformation (power-law rheology) or on common simplifications of lithospheric materials (Newtonian rheology; see Table 4.1 in text). To verify that the numerical procedure properly models viscoelastic deformation, we simulated postseismic surface deformation following a kinematically imposed earthquake, using a model in which no damage evolution was allowed. We compared our results with those from a viscoelastic finite-element model (GAEA, Saucier et al., 1992). This extends previous comparisons which yielded interseismic velocities comparable to analytic solutions (Lyakhovsky et al., 2001). The two models were set up with identical domain size (300 × 250 × 100 km), fault depth (20 km), coseismic stress drop (20 MPa), lithospheric structure, and material properties (=30 GPa, =0.25, elastic plate thickness 20 km, and Newtonian substrate viscosity of 1019 Pa s). In both models the fault-parallel side and bottom boundaries were fixed, and the top boundary was stress-free. To simulate stress-free fault-perpendicular end boundaries our damage-disabled model incorporated the variable force boundary condition with a low stiffness, which is approximately equivalent to the stress-free boundary condition used in the FEM. Figure C.1 147  shows surface velocities calculated with both models, for several time epochs after the modeled earthquake (1, 2, 5, and 18 years). The minor differences between the models are due in part to applying the variable force boundary condition with spring stiffness slightly too high (103 MPa) and therefore not simulating completely stress-free end boundaries. This slightly suppressed the velocities away from the fault.  Figure C.1. Surface velocities from our damage-disabled model (points) compared to a finite element solution of viscoelastic post-seismic deformation (solid lines). The various lines and points represent velocities after 1 year (red), 2 years (black), 5 years (purple) and 18 years (green) of postseismic deformation. Maxwell time m=110 years.  148  Bibliography Brace, W. F. and D. L. Kohlstedt, 1980. Limits on lithospheric stress imposed by laboratory experiments. J. Geophys. Res., 85: 6,248-6,252. Goetze, C., 1978. The mechanisms of creep in olivine. Philos. Trans. R. Soc. London A, 288: 99119. Karato, S. I., 1989. Seismic anisotropy: mechanisms and tectonic implications. In: S. I. Karato and M. Toriumi (eds), Rheology of Solids and of the Earth. Oxford Univ. Press, New York, 393-422 p. Kirby, S. H., 1983. Rheology of the lithosphere. Rev. Geophys., 21: 1458-1487. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon, 2001. Earthquake cycle, faults, and seismicity patterns in rheologically layered lithosphere, J. Geophys. Res., 106: 4103-4120. Saucier, F., Humphreys, E.D., and R.J. Weldon, 1992. Stress near geometrically complex strikeslip faults: application to the San Andreas Fault at Cajon Pass, southern California, J. Geophys. Res., 97, 5081–5094. Weertman, J., 1978. Creep laws for the mantle of the Earth. Philos. Trans. R. Soc. London A, 288: 9-26.  149  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052947/manifest

Comment

Related Items