Viscoelastic finite-element models of the earthquake cycle along plate-boundary faults by Ali Vaghri 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) August, 2011 © Ali Vaghri 2011 Abstract In this thesis I describe three related projects on the mechanics of plate bounding faults and fault systems. All of the projects involve viscoelastic earthquake cycle modeling, using the finite element method, and investigating how surface deformation is affected by viscoelastic relaxation or viscous shear zone creep. In the first project I develop test models to study how contrasts in material properties across a fault can affect surface deformation. I model contrasts in elastic plate thickness and in viscosity profiles across a strike-slip fault. I find that an offset between the surface trace of a fault and its position at depth (where it creeps), due to a non-vertical dip, may make interseismic deformation look asymmetric. Rheological contrasts do not yield dramatic asymmetry for models with realistic viscosity profiles and the sense of asymmetry in viscoelastic models actually reverses during the interseismic interval. In the second project, I model a set of parallel strike-slip faults, which allow relative motion of the North American and Pacific Plates in northern California. Geologic information allows me to take the timing of large earthquakes on these faults into account in the modeling, and I compare modeled surface deformation with GPS data from the region. I find that to explain the GPS velocity field, I must incorporate the dip of the faults into my model. The required dip for the Green Valley Fault in particular is consistent with new, double-difference hypocenter data from the region. The final project involves modeling Cascadia Subduction Zone earthquake cycles. This subduction zone has been modeled in the past, but always using a kinematic device (―backslip‖) for modeling slip on the subduction interface. I developed a detailed profile model of the Cascadia Subduction Zone, based on data from the Lithoprobe project. Since the Subduction Zone has a curved profile and its slip is driven by relative motion of the model boundaries, the models produced unrealistic uplift along the North America plate boundary. I tested several approaches for preserving the geometry of the subduction ii interface. I suggest that a new implementation of split nodes in the finite-element method would be required to make this correction work. iii Preface Chapter 2 of this dissertation has been accepted for publication in the Bulletin of the Seismological Society of America (―Can lateral viscosity contrasts explain asymmetric interseismic deformation around strike-slip faults?‖, Manuscript D-10-00347, by A. Vaghri and E. Hearn). I was responsible for developing the meshing codes and models, and for testing the central hypothesis of this paper. I also prepared all figures and tables and wrote the manuscript. My PhD advisor, Elizabeth Hearn, edited the final manuscript and supplied the initial idea for this research project. Chapter 3 of this dissertation will be submitted this summer (2011) for publication to a long-format geophysical journal. It will be entitled ―Earthquake cycle models of the PacificNorth America Plate Boundary at Point Reyes, California‖, and the authors will be A. Vaghri and E. Hearn. For this project, I was again responsible for developing the meshing codes and models, and for testing the central hypothesis of this paper. I also prepared all figures and tables and wrote the manuscript. Elizabeth Hearn edited the final manuscript, and Dr. Wayne R. Thatcher of the United States Geological Survey suggested this research project. Chapter 4 presents a detailed exploration of my effort to model an earthquake-cycle problem in which the fault orientation was not parallel to the model side boundaries (the Cascadia subduction zone), without resorting to the back-slip technique. I did not find a satisfactory solution to this technical problem. It is my hope that by presenting my work at past meetings (e.g., the annual SCEC Community Finite Element Modeling Workshops in Golden Colorado), and by outlining it in this dissertation, I may save other modelers time they might have spent on similar approaches to this problem. iv Table of Contents Abstract .................................................................................................................................................. ii Preface ................................................................................................................................................... iv Table of Contents ................................................................................................................................... v List of Tables ....................................................................................................................................... viii List of Figures ....................................................................................................................................... ix Acknowledgements ............................................................................................................................... xi Chapter 1 Introduction ........................................................................................................................ 1 1.1 Plate boundaries...................................................................................................................... 1 1.2 Mechanical models of faults and deformation ....................................................................... 2 1.2.1 Geometric complexity and material heterogeneity ......................................................... 3 1.2.2 Time-dependence of deformation: the earthquake cycle ................................................ 4 1.3 Global positioning system (GPS) ........................................................................................... 5 1.4 Thesis outline ......................................................................................................................... 6 Chapter 2 Can lateral viscosity contrasts explain asymmetric interseismic deformation around strike-slip faults? .................................................................................................................................... 8 2.1 Introduction ............................................................................................................................ 8 2.2 Models and methods ............................................................................................................... 9 2.2.1 Asymmetric meshes...................................................................................................... 11 2.2.2 Spin-up and benchmarking ........................................................................................... 14 2.3 Results .................................................................................................................................. 15 2.3.1 Lateral contrasts in plate thickness: Effects on velocity profiles.................................. 15 2.3.2 Lateral contrasts in viscosity structure ......................................................................... 17 2.3.3 Inferred elasticity contrasts from forward-modeled profiles ........................................ 19 2.3.4 Locally thinned plate models........................................................................................ 23 2.4 Discussion and conclusions .................................................................................................. 26 2.4.1 Comparison with previous studies................................................................................ 26 2.4.2 Ramifications: Likely causes of deformation asymmetry ............................................ 29 2.5 Data and resources ................................................................................................................ 30 Chapter 3 Earthquake cycle models of the Pacific-North America plate boundary at Point Reyes, California .......................................................................................................................................... 31 3.1 Introduction .......................................................................................................................... 31 v 3.2 Data ...................................................................................................................................... 31 3.2.1 Slip rates, creep, interseismic intervals, and timing of earthquakes ............................. 32 3.2.2 Material properties and fault geometry ........................................................................ 36 3.2.3 GPS data ....................................................................................................................... 38 3.2.4 Seismicity data.............................................................................................................. 40 3.3 Method.................................................................................................................................. 41 3.4 Results .................................................................................................................................. 45 3.4.1 Locking depth ............................................................................................................... 45 3.4.2 Viscosity of shear zone ................................................................................................. 46 3.4.3 Viscosity of lower crust and upper mantle ................................................................... 47 3.4.4 West Napa Fault ........................................................................................................... 49 3.5 Discussion ............................................................................................................................ 51 3.6 Conclusions .......................................................................................................................... 56 Chapter 4 Earthquake cycle models of the Cascadia subduction zone ............................................. 64 4.1 Introduction .......................................................................................................................... 64 4.2 Subduction zones .................................................................................................................. 64 4.2.1 The Cascadia subduction zone ..................................................................................... 66 4.2.2 Deformation models of subduction zone ...................................................................... 67 4.3 Data ...................................................................................................................................... 70 4.3.1 Seismic wave velocity and density data ....................................................................... 71 4.3.2 Calibration data: GPS velocities and long-term uplift rate ........................................... 73 4.4 Finite element modeling of the Cascadia subduction zone .................................................. 74 4.4.1 Elastic structure ............................................................................................................ 74 4.4.2 Geometry of the Cascadia subduction interface ........................................................... 76 4.4.3 Meshing the model ....................................................................................................... 77 4.4.4 2D Cascadia elastic model: Mesh testing and benchmarking ...................................... 78 4.4.5 Incorporating viscoelastic relaxation and stress-driven slip ......................................... 80 4.4.6 The problem of tectonic loading and displacement of the subduction interface .......... 84 4.4.7 Corrections to the FE model equations to preserve the geometry of the subduction interface ...................................................................................................................................... 88 4.5 4.5.1 Recommended correction method and future research ........................................................ 95 Physics of aseismic deformation along the subduction zone........................................ 96 vi 4.5.2 Evolution of stresses throughout the earthquake cycle, and seismic hazard ................ 97 4.5.3 Inferred locations of episodic tremor and slip events ................................................... 98 Chapter 5 Conclusion ........................................................................................................................ 99 5.1 Contributions ........................................................................................................................ 99 5.1.1 Asymmetric deformation around strike-slip faults is not due to material contrasts ..... 99 5.1.2 Incorrect assumptions about fault geometry may bias slip rate estimates from geodetic data .................................................................................................................................... 100 5.1.3 Subduction zone earthquake cycle models without backslip: still an unresolved problem .................................................................................................................................... 101 5.2 Limitations.......................................................................................................................... 102 5.3 Future work ........................................................................................................................ 103 Bibliography ....................................................................................................................................... 105 Appendix A Finite element modeling ............................................................................................. 119 Appendix B Mesh software for automatic mesh generation ........................................................... 127 Appendix C Lagrange multiplier .................................................................................................... 133 Appendix D Input file for mesh generating code (Mesh) ............................................................... 135 vii List of Tables Table 2.1. Model parameters ...................................................................................................................... 13 Table 2.2. Inferred locking depth and elasticity contrast (asymmetric models) ......................................... 19 Table 2.3. Inferred locking depth and slip rate (symmetric models) .......................................................... 25 Table 3.1. Slip rate and earthquake timing values. ..................................................................................... 35 Table 3.2. Fault parameters ......................................................................................................................... 56 Table 3.3. Velocities of continues GPS sites ............................................................................................. 58 Table 3.4. Campaign GPS site velocities .................................................................................................... 63 viii List of Figures Figure 2.1. Interseismic GPS velocity profiles across .................................................................................. 8 Figure 2.2. Model diagrams and parameters ................................................................................................. 9 Figure 2.3. Model mesh (profile view) ....................................................................................................... 10 Figure 2.4. Benchmarking results for a symmetric viscoelastic test model ................................................ 12 Figure 2.5. Effct of the plate thickness contrast .......................................................................................... 15 Figure 2.6. Effect of the substrate viscosity on models with a plate-thickness contrast ............................. 16 Figure 2.7. Effect of transition width w on models with a plate-thickness contrast ................................... 17 Figure 2.8. Velocity profiles for models with a viscosity contrast ............................................................. 18 Figure 2.9. Velocity profiles for nonlinear viscoelastic models ................................................................. 22 Figure 2.10. Velocity and strain rate profiles for thinned-plate models ..................................................... 24 Figure 3.1 Study area, showing active faults .............................................................................................. 32 Figure 3.2. Rock types and material property of the blocks ....................................................................... 38 Figure 3.3. GPS velocities north of the San Francisco Bay ........................................................................ 39 Figure 3.4. Profile of plate-boundary parallel GPS site velocities .............................................................. 40 Figure 3.5. Earthquake epicenters within 45 km of the model profile ........................................................ 41 Figure 3.6. Model schematic and mesh in the vicinity of the model .......................................................... 43 Figure 3.7. WRSS of the model as function of lucking depth .................................................................... 45 Figure 3.8. WRSS of models within a 10 km locking depth....................................................................... 47 Figure 3.9. Contour plot of WRSS for different viscosities of the lower crust and mantle ........................ 48 Figure 3.10. Model velocity profiles, showing sensitivity of WRSS to location of the GVF..................... 49 Figure 3.11. Sensitivity of WRSS to relative slip rate for the NF and the GVF. ....................................... 50 Figure 3.12. A comparison of elastic dislocation and viscoelastic finite element models .......................... 52 Figure 3.13. Epicenters of earthquakes within 45 km of the modeled profile ............................................ 54 .Figure 3.14. Earthquake hypocenters projected onto the profile plane. ..................................................... 55 Figure 4.1. Schematic illustration of a subduction zone ............................................................................. 65 Figure 4.2. Geometry of the Cascadia Subduction Zone ........................................................................... 66 Figure 4.3. Profile showing lithosphhere densities for the Cascadia subduction zone .............................. 71 Figure 4.4. Location of the project profile across the Canadian Southwest Cordillera .............................. 72 Figure 4.5. Locations of GPS sites on Vancouver Island and southwest British Columbia ...................... 73 Figure 4.6. Bulk modulus values for the modeled profile.......................................................................... 75 Figure 4.7. Modeled coseismic surface subsidence .................................................................................... 76 Figure 4.8. Part of the unstructured Cascadia profile mesh ........................................................................ 78 Figure 4.9. A comparison of numerical and analytical solutions ................................................................ 80 ix Figure 4.10. Schematic trajectory of split node displacement over a seismic cycle ................................... 86 Figure 4.11. Elastic displacement, model-net (permanent) displacement and total displacement .............. 87 Figure 4.12. Schematic of the model when it is separated to two parts. ..................................................... 89 Figure 4.13. Vertical displacement of the model with displacement correction. ........................................ 90 Figure 4.14 Vertical uplift of the model with force correction term on interface nodes............................. 92 Figure 4.15. Schematic of elements which are touching each other on the fault interface ......................... 93 Figure 4.16 Vertical displacement of model with force correction term on the all nodes of the elements which touch the interface. ........................................................................................................................... 94 Figure 4.17 vertical displacement for the model with force correction term computed by modeling slab with boundary condition on the interface nodes. ........................................................................................ 95 x Acknowledgements First and foremost, I'd like to thank my supervisor, Dr. Elizabeth Hearn, for her generous support throughout my time at the graduate school. Dr. Hearn's unique way of supervision always gave me enough freedom to shape my research, while providing me with insightful feedback and encouragement to bring this work to completion. Many thanks go to the members of my doctoral supervisory committee, Dr. Garry Clarke and Dr. Ronald Clowes. Their advice was always helpful, particularly during Dr. Hearn's sabbatical, when I needed it most. It's been a real pleasure to have been a member of the geophysics group and the department of Earth and Ocean Sciences at UBC. I have had a marvelous time with amazing people in the department. I am thankful to my friends and office mates: Dr. Jounada Ouety, Dr. Yaron Finzi, and Joel Podgorski. I would particularly like to acknowledge Dr. Phil Hammer, Dr. Christian Schoof, Dr. Timothy Creyts and Dr. Valentina Radic. I am overwhelmed by the everlasting support and encouragement of my parents, Salman Vaghri and Arefeh Aghajani and my in-laws Dr Aziz Yahyavi and Minoo Ameli . I am grateful to my caring wife Ladan and our two beautiful daughters Epak and Hannah, who have always filled our home with love and laughter. I am appreciative of my sister Dr. Ziba Vaghri, whose help was pivotal during the process of my admission to UBC. My gratitude also goes to my other two sisters, Simin and Nastaran Vaghri, for their support. xi To Ladan and my lovely daughters Epak and Hannah. xii Chapter 1 Introduction Tectonic plates on the earth‘s surface continuously move relative to each other, and the type of deformation at their boundaries depends on the mechanical properties of the plates and their sense of relative motion. Although the relative motion of plates is fairly steady at timescales of thousands of years, it may be strongly time dependent at plate boundaries. This is because of friction in the upper crust; plates are locked together in the upper crust. As relative motion continues, elastic stresses build until the locked interval releases stress in an earthquake. Several million earthquakes occur in the world every year and cause huge loss of life and property. In 2010, earthquakes killed 296,800 people, affected nearly 208 million others, and cost US$110 billion (UN report, 24 January 2011). Understanding the mechanism of a natural disaster is the very first, important step in the process of saving lives and resources. To study earthquakes, we need to study the geometry and behavior of active faults, as well as the mechanical properties of the earth‘s outermost layers, particularly in plate boundary regions. 1.1 Plate boundaries According to plate tectonic theory, there are two types of convergent plate boundaries. In subduction zones, either an oceanic or a continental plate overrides a subducting oceanic plate. At collisional boundaries, on the other hand, lithospheric plates move toward each other without subducting and the resulting compression causes thickening of the crust. Subduction zones can experience both shallow and deep earthquakes, and are the sources of very strong megathrust earthquakes, such as the 1960 magnitude 9.5 Chile earthquake, the 2004 magnitude 9.3 Sumatra earthquake, and the recent magnitude 9.0 earthquake in northern Japan. Megathrust offshore earthquakes may generate tsunamis which can cause far more casualties and damage than the earthquake itself, as seen recently in Japan and Sumatra. The active crustal deformation zone surrounding the Pacific Ocean, known as the Ring of Fire, is the locus of the majority of the earth‘s subduction zones, for example, the Cascadia subduction zone at the boundary of the Juan de Fuca and North America plates, the Chilean subduction zone at the boundary of the Nazca 1 and South America plates, and the Japan Trench at the boundary of Pacific and North America plates in northern Japan. Transform boundaries occur where two plates laterally slide past each other and create a shearing force as they move. Most transform faults are hidden in the deep oceans where they form a series of short zigzags accommodating relative motion of the seafloor and displacing the spreading ridge. A smaller number of transform faults, which are more well-known, cut continental lithosphere; an example of this is the San Andreas Fault of western North America. The San Andreas fault connects a divergent boundary in the Gulf of California with the Cascadia subduction zone. Another example of a transform boundary on land is the Alpine Fault of New Zealand, which is a strike slip fault between the Australian and the Pacific plates. Continental transform faults can produce moderate to large, devastating earthquakes. Recent examples include the magnitude 7.1 and 6.3 Darfield and Christchurch New Zealand earthquakes, and the magnitude 7 Haiti earthquake in 2010. Divergent plate boundaries, where two plates move away from each other, are typically the source of small to moderate earthquakes. Examples include the East Africa Rift and mid-ocean ridges. In this dissertation, I focus on transform and subduction plate boundaries. 1.2 Mechanical models of faults and deformation At the simplest level, we think of plate-boundary faults as slipping in earthquakes between the surface and a depth of about 20 km, and creeping aseismically at a constant rate at greater depths. Simple analytical solutions for deformation around a dislocation in an elastic halfspace (for example, Freund and Barnett, 1976; Savage and Burford, 1973) may be used to calculate displacements (hence velocities) of points on the surface, in the case of a single, long fault. However, the Earth is more complicated than this. Geometrically complex networks of faults in non-uniform lithosphere accommodate the relative motion of plates, and between earthquakes, the relative motion is not steady-state. 2 1.2.1 Geometric complexity and material heterogeneity Some tectonic plate boundaries are simple, because relative motion of the plates is accommodated mainly by a single fault. Examples include the Cascadia Subduction Zone west of British Columbia, the central section of the San Andreas Fault in California, the North Anatolian Fault in Turkey, and the Dead Sea Fault in the Middle East. Ocean-ocean plate boundaries are also fairly simple, with mid-ocean ridge segments separated by transform faults. Many continental plate boundaries are complex, including the Pacific-North America Plate boundary (i.e. western North America), where relative plate motion occurs along hundreds of active faults within about 300 km of the plate boundary. Other classic examples of diffuse plate boundaries are the India-Eurasia collisional plate boundary, as well as parts of the Andes, New Zealand, and Japan (e.g., Gordon, 1998). In diffuse plate boundaries, it is hard to untangle the contributions of distinct faults to the GPS velocity field, and to infer slip rates on individual faults. Debate continues on whether fault zones extend into the mantle as discrete surfaces, or whether the lithosphere deforms as a continuum below the upper crust, in part because of this ambiguity (Thatcher, 1995). For a vertical shear dislocation creeping below a ―locking depth‖ of D in an elastic Earth, the width of the zone accommodating most of the strain is about 2D (Savage and Burford, 1973). If D is 20 km, then for faults spaced 20 km apart, untangling individual slip rates becomes a challenge. Elastic block models are traditionally used to infer slip rates and locking depths in geometrically complex areas (e.g., Bennett et al., 1996), but non-unique solutions for slip rate and D are common. Hence we need other strategies to understand what is really occurring. One is using both geologic and geodetic information to solve for slip rates. Another is to base locking depths on the maximum depth of seismicity, rather than solving for this. Spatial variations in elasticity of the lithosphere may also alter surface deformation associated with faulting (e.g., Le Pichon et al., 2005). Elastic block models do not usually account for this heterogeneity. 3 1.2.2 Time-dependence of deformation: the earthquake cycle Another issue with the traditional use of elastic block models to infer slip rates and locking depths on faults is that deformation should vary with time between large earthquakes. We expect that at high temperatures and pressures (e.g. in the middle crust and below), geologic materials transition from brittle-elastic to viscoelastic behaviour. If a viscoelastic material is present in the lithosphere, its strain rate depends on stress, so perturbations to stress by earthquakes may affect GPS strain rates measured at the Earth‘s surface. Creeping frictional faults will also respond to a stress perturbation with a change in sliding speed. The result is accelerated deformation just after an earthquake (postseismic deformation) and slower deformation late in the interseismic interval than one would predict using a steady-state, elastic dislocation model to represent the fault. This can lead modelers to infer slip rates incorrectly among active faults in a tectonically complex region. An understanding of the mechanics of faults is crucial for knowing why geologic and GPSbased fault slip rates sometimes disagree with each other. Slip rate for the Owens Valley fault zone in eastern California based on geodetic data and elastic half-space models, 5-7 mm/yr, are faster than longer term geologic rates, 2-3 mm/yr (Dixon et al., 2003). Slip rates in the San Andreas fualt in southern California based on block model constrained by GPS data varies from 35.9 0.5 mm/yr in the Carrizo Plain to 5.1 1.5 mm/yr through the San Benadino segment (Brendan and Bradford, 2005). Reilinger et al., (2006) present a GPS-derived velocity field for the 16 faults of the zone of interaction of the Arabian, Africa and Eurasian plates and show their difference with geological slip rate. On the base of geodetic data slip rate of North Anatolian and East Anatolian faults are respectively 30 2 and 15 3 mm/yr where based on geological estimation these values respectively are 5-15 and 4-7 mm/yr (Reilinger et al., 1997). Thatcher (2009) shows huge discrepancy between GPS velocity rate and geologic slip rate of faults in Tibet. Analytical solutions exist for time-dependent interseismic surface deformation around an infinitely long fault in an elastic plate above a Maxwell viscoelastic halfspace (e.g., Savage and Prescott, 1978). Analytical solutions have also been developed for more complex linear rheologies and for irregular sequences of earthquakes (e.g., Hetland and Hager, 2005). For lithosphere with spatially variable properties, or nonlinear (stress-dependent) viscosity, 4 interseismic deformation must be modeled numerically. Modeling time-dependent interseismic deformation for a complex network of faults remains a great challenge for our field. Typically this is done with the boundary element method, which restricts the amount of heterogeneity which can be modeled, and which does not deal with nonlinear rheology (e.g. Johnson and Segall, 2005). To account for material heterogeneity and nonlinear rheology, as well as complex fault geometries, the finite-element method is required. 1.3 Global positioning system (GPS) The Global Positioning System (GPS) is a space-based global navigation system made up of a network of 32 satellites which provides location and time information. It is operated by the U.S. Department of Defense for military purposes, but it is limitedly accessible by anyone with a GPS receiver. The GPS project was developed in the 1970‘s to overcome the limitations of previous navigation systems. The current GPS consists of three major components. These are the space component, a control component, and a user component. The space component consists of the orbiting GPS satellites, the control component synchronizes the satellites‘ atomic clocks and adjusts the parameters of each satellite's orbital model, and the user component is military and civil GPS receivers. (http://www.gps.gov/systems/gps/) GPS satellites continually transmit messages which include the time of the message, precise orbit information, satellite health and general information. A GPS receiver calculates its position by ranging (i.e., from signal travel times and satellite positions). It must be locked on to the signal of at least three satellites to calculate a 2D position (latitude and longitude) and track movement. (http://www.colorado.edu/geography/gcraft/notes/gps/gps_f.html) Improvements to receiver technology and data analysis procedures have improved the accuracy of GPS in recent years. One of the most significant developments is the creation of the International GPS Service for Geodynamics. This service provides data to accurately post process GPS data. Now the accuracy of horizontal coordinates of sites appears to be <10 mm on 5 a global scale. On a regional scale, precisions of 1-3 mm for horizontal components and <10 mm for the vertical component are available for stable GPS sites (Herring, 1999). Nowadays GPS is a very important geodetic tool which is used widely in geoscience. One application is a continually operating global network of GPS sites which is used to study global earth dynamics. Local GPS networks (which may comprise both permanent and survey-mode sites) provide data on current crustal deformation on a regional scale. One target of such networks, and the topic of this dissertation, is fault-related deformation. 1.4 Thesis outline In Chapter 2, I develop a suite of finite-element models to investigate asymmetric interseismic deformation around strike-slip faults. This work is inspired by the fact that major faults can juxtapose different lithosphere (for example, continental and oceanic). It is also inspired by recent research suggesting that asymmetric deformation across strike-slip faults is due to dramatic contrasts in elastic properties. In this project I study the effect of contrasts in plate thickness and viscosity on asymmetric deformation. Chapter 3 focuses on the Point Reyes area located north of the San Francisco Bay in California. I model a set of parallel strike-slip faults which allow relative motion of the North American and Pacific Plates in this region. Because several parallel strike-slip faults are active, large earthquakes on these faults occur over different interseismic intervals and are staggered in time. I vary fault parameters such as locking depth and the viscosity of shear zones that extend these faults into the mantle, as well as how the Pacific-North America slip rate is partitioned among the faults. I also test several values of effective viscosity for the lower lithosphere. In addition to the geodetic data, I use seismic data to more accurately define the fault geometry. Chapter 4 deals with modeling of another form of plate boundary on the west coast of North America. This chapter is a study of the Cascadia subduction zone on Vancouver Island where Juan de Fuca oceanic plate slips beneath the North America continental plate. Since the forces driving flexure of the subducting plate are not implicitly modeled, the models produce unrealistic long-term uplift of the continental plate. In this chapter, I describe my efforts to add correction 6 terms to the finite-element force balance equations in order to correctly model earthquake-cycle deformation with stress-driven interseismic slip, while not explicitly including every force driving flexure of the subducting plate. 7 Chapter 2 Can lateral viscosity contrasts explain asymmetric interseismic deformation around strike-slip faults? 2.1 Introduction Asymmetric surface velocities have been documented around the San Andreas Fault (Schmalzle et al., 2006, Fialko, 2006, Lisowski et al., 1991, Prescott and Yu, 1986) and several other major, essentially vertical strike-slip faults (e.g. the North Anatolian and Sumatra faults, Le Pichon et al., 2005; and the Altyn Tagh Fault, Jolivet et al., 2008). (Figure 1 shows profiles of GPS fault-parallel velocities along transects crossing the San Andreas and North Anatolian Faults.) This asymmetry has often been interpreted in terms of an elasticity contrast across the fault, due to the juxtaposition of different materials (e.g. Lisowski et al., 1991; Le Pichon et al., 2005; Jolivet et al., 2008; Schmalzle et al., 2006). Figure 2.1. Interseismic GPS velocity profiles across (a) the San Andreas Fault at the Carrizo Plain and (b) the North Anatolian Fault in the Eastern Marmara Sea. Profiles show fault-parallel velocities for sites located in a 100-km-wide swath crossing each fault. Arctangent velocity profiles for uniform and asymmetric (quarter-space) elastic models are shown (Savage and Burford, 1973; Rybicki and Kasahara, 1977). For the San Andreas Fault and North Anatolian Fault, we assume a locking depth of 18 km and slip rates of 35 and 24 mm/yr, respectively. Thin, heavy, and dotted lines represent elasticity ratios (left:right) of 0.5:1, 1:1, and 2:1, respectively. 8 Contrasts in elastic plate thickness (e.g., Huang and Johnson, 2011) or substrate viscosity (Malservisi et al., 2001, Lundgren et al., 2009) across the fault have also been suggested as explanations, as well as a horizontal offset in the position of the creeping fault at depth relative to the surface trace (Jolivet et al., 2008; Fialko, 2006; Meade et al., 2002). Here, we use viscoelastic earthquake-cycle models to assess systematically how lateral contrasts in plate thickness and viscosity structure influence surface deformation. For a variety of parameter ranges we estimate the extent to which the resulting asymmetric deformation would be interpreted in terms of contrasts in elasticity across the fault. We also address the extent to which a locally thinned elastic plate affects surface deformation around a strike-slip fault (and hence, the inferred slip rate and locking depth). 2.2 Models and methods We used the code Geophysical Finite Element Simulation Tool (GeoFEST), version 4.7, for finite-element modeling in this study (Parker et al., 2003). The models comprise an elastic layer overlying a viscoelastic half space (which may be rheologically stratified), and are designed with lateral contrasts in elastic plate thickness and viscosities. Figure 2.2 provides a schematic of the three suites of models we address. The FE meshes extend 500 km horizontally from the fault and are 500 km deep to minimize the effect of imposed boundary conditions on the solutions. They extend 100 km in the along-strike direction. Figure 2.2. Model diagrams and parameters. Diagrams show (a) plate thickness contrast models (b) effective viscosity contrast models (depth-dependent and power-law viscosities are also modeled) (c) thinned-plate models. 9 Hexahedron (brick) elements are used, and elements range in dimension from 1 km near the fault to tens of km in the far field (Figure 2.3). In all of the models, an infinite strike slip fault fully penetrates the elastic plate and the side boundaries move at a relative rate of 40 mm/year parallel to the fault. Earthquakes are modeled at 200-year intervals (TC = 200 years). The bottom model boundary is free to move horizontally and vertically, and the top boundary is free (i.e., tractionless). For all models, elasticity is uniform, with Poisson‘s ratio =0.25 and shear modulus G=40 GPa. To benchmark some of the GeoFEST results, we developed a viscoelastic earthquake cycle model identical to the laterally homogeneous models presented here, using the finiteelement code GAEA (Saucier and Humphreys, 1993). Figure 2.3. Model mesh (profile view). (a) and (b): model with plate thickness transition width (w) of 70 km, D1 = 10 km, and D2 = 25 km. (c) and (d) model with a locally thinned plate. 10 2.2.1 Asymmetric meshes Three categories of earthquake-cycle models were developed for this study (Figures 2.2 and 2.3). For one set of models, the plate thickness varies across the fault (D1 on the left, at x < 0 and D2 on the right, at x > 0), with a linear taper from D1 to D2 over a distance of w. For each scenario, we tested uniform viscosity values of 3.5x10, 3.5x10, and 3.5x10Pa s. Given the shear modulus (40 GPa) these correspond to Maxwell times (TM = 2 η /G) of 5, 50, and 500 years. For a 200-year earthquake cycle, the Savage parameter (TC/TM) values are 40, 4, and 0.4, respectively. Values of D1, D2, w, and η are given in Table 1, together with the model names. For another set of models, localized thinning of the plate around the fault (inspired by Chery, 2008) was tested. For these models, we did not vary the geometry of the plate-thickness ―divot‖, which is 40 km wide, 20 km deep, and centered at the fault (see Figure 2.2c). Viscosity values 18 19 20 of 3.5x10 , 3.5x10 , and 3.5x10 Pa s were modeled in the substrate beneath the plate. We also modeled the effects of viscosity contrasts for layered linear Maxwell and power-law viscoelastic rheologies, using a suite of models in which a contrast in viscosity (or power-law parameter AT) was modeled at x = 0. For power-law flow, the stress-strain rate relation (e.g. Bürgmann and Dresen, 2008) is of the form: A dn d m f Hr O e Q PV RT 2 (2.1) A is a material constant, d is differential stress, n is the stress exponent, fH2O is water fugacity, r is the fugacity exponent, d is the grain size, m is the grain size exponent, Q is the activation energy, V is the activation volume, R is the molar gas constant, and T is absolute temperature. We define the power-law coefficient AT as: AT A 1d m f H2rO e Q PV RT (2.2) - This gives 11 dn AT (2.3) This expression may be rewritten in terms of effective viscosity. e AT d1n (2.4) Models with AT values ranging from 1036 to 1037 Pan s are presented here (n = 3). Models with higher AT did not yield cycle-invariant results within 100 earthquake cycles and models with lower AT were numerically unstable for our model meshes and time-stepping parameters (minimum step length of one day and time step multiplier set by GeoFEST). In the power-law models, the plate thickness D was set at 16 km. Values of viscosity for the Maxwell viscoelastic models are shown on Table 2.1. For the power-law models, Table 2.1 provides values of AT and stress exponent n. Figure 2.4. Benchmarking results for a symmetric viscoelastic test model. Model parameters are given in text. The Maxwell time (TM) is 50 years and the interseismic interval (Tc) is 200 years. Velocity profiles are shown for seven time snapshots at twenty-year intervals (i.e., t/Tc = 0.1 to 0.7). 12 Table 2.1. Model parameters Layered Models Model D1 (km) D2 (km) w (km) η (Pa s) CC001 CC002 CC003 10 10 10 10 10 10 n/a n/a n/a 3.5e18 Pa s 3.5e19 Pa s 3.5e20 Pa s Model D1 (km) D2 (km) w (km) η (Pa s) CC111 CC121 CC131 CC112 CC122 CC132 75CC113 CC123 CC133 10 10 10 10 10 10 10 10 10 16 16 16 16 16 16 16 16 16 30 50 70 30 50 70 30 50 70 3.5e18 Pa s 3.5e18 Pa s 3.5e18 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e20 Pa s 3.5e20 Pa s 3.5e20 Pa s CC211 CC221 CC231 CC212 CC222 CC232 75CC213 CC223 CC233 10 10 10 10 10 10 10 10 10 25 25 25 25 25 25 25 25 25 30 50 70 30 50 70 30 50 70 3.5e18 Pa s 3.5e18 Pa s 3.5e18 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e20 Pa s 3.5e20 Pa s 3.5e20 Pa s CC311 CC321 CC331 CC312 CC322 CC332 75CC313 CC323 CC333 10 10 10 10 10 10 10 10 10 40 40 40 40 40 40 40 40 40 30 50 70 30 50 70 30 50 70 3.5e18 Pa s 3.5e18 Pa s 3.5e18 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e19 Pa s 3.5e20 Pa s 3.5e20 Pa s 3.5e20 Pa s Model D (km) η1 (Pa s) η2 (Pa s) interval (km) 75AV13 15 3.5e18 3.5e20 >10 75AV23 15 3.5e19 3.5e20 >10 75AV3L13 15 3.5e18 3.5e20 3.5e17 3.5e19 3.5e21 3.5e18 15 to 40 40 to 60 >60 75AV3L23 15 3.5e18 3.5e20 3.5e17 3.5e20 3.5e22 3.5e19 15 to 40 40 to 60 >60 Models with plate thickness contrast Models with Viscosity Contrast 13 Table 2.1. Model parameters, continued. Models with Locally Thinned Plate Model D (km) η1 (Pa s) 20C19 40C20 75C21 see text and Fig. 2 3.5e18 3.5e19 3.5e20 Power Law Models Model D (km) n (km) log AT(1) n (Pa s) log AT(2) n (Pa s) 100NAV66 15 3 36 36 100NAV67 15 3 36 37 _______________________________________________________________________________ Note: AT is defined in Equation 2. A model with AT = 37 on both sides was also run, but it did not achieve a cycleinvariant state. In general 40 earthquake cycles were run (20 for η = 3.5e18 Pa s). Leading numbers in some filenames indicate the number of earthquake cycles that were run. 2.2.2 Spin-up and benchmarking Figure 2.4 shows a comparison of the analytical solution to two FE solutions for interseismic surface velocities around a fault. The results are for a laterally homogeneous model with a 19 substrate viscosity of 3.5x10 Pa s. Both of the finite-element models have difficulty reproducing the earliest near-field interseismic velocities predicted by the Savage and Prescott (1978) solution. GAEA performs better for small values of t/Tc and GeoFEST performs better later in the earthquake cycle. Finer discretization at the fault in the vicinity of the relaxing layer can improve the model performance (e.g. Lundgren et al., 2009), though for higher values of effective viscosity smaller elements strain so much that numerical precision is compromised after several earthquake cycles. Using the mesh shown in Figure 2.3, we obtain a better match to the 20 analytical solution for 3.5x10 Pa s but for higher values (not presented) large element strains prevent convergence of the solution. Our mesh is optimized for the range in Savage parameter values (TC/TM) we explore, and though our postseismic velocity and strain rate estimates are somewhat low for some models, comparisons among different models should not be affected. 14 2.3 Results 2.3.1 Lateral contrasts in plate thickness: Effects on velocity profiles Effects on velocity profiles for the suite of models in which we assess the effects of a plate thickness contrast across the fault, we vary plate thickness contrast D1/D2 (renamed Dr), width over which the contrast occurs (w), and substrate viscosity (η). D1 is 10 km for all of the models. Figure 2.5. Effct of the plate thickness contrast. Modeled velocity profiles are shown for models with three values of D2. η = 3.5 x 1019 Pa s and w = 30 km. (a) velocity profiles in modeled reference frame (i.e., with constant velocities at side boundaries). (b) For D2 = 40 km, velocity profiles relative to a point on the fault. The heavy black line shows the velocity profile for a Savage-Burford style dislocation model with Vo = 40 mm/yr and locking depth (zL) equal to 25 km (the plate thickness at the fault). 19 Figure 2.5 shows three models in which the substrate viscosity is 3.5x10 Pa s, and only Dr is varied. As one would expect, a larger plate thickness contrast results in more asymmetric velocity profiles (Figure 2.5a). Note that when the model side boundary velocities are symmetric about the fault (+Vo/2 and -Vo/2, respectively), the velocity of a point on the fault varies throughout the earthquake cycle. This is necessary to reconcile symmetric coseismic deformation with asymmetric interseismic deformation. Velocity profiles plotted relative to a point on the fault (Figure 2.5b), show that the sense of asymmetry relative to the fault reverses over the earthquake cycle. Initially, the highest velocities 15 (relative to a point on the fault) are on the side with the thin plate. Late in the cycle, this reverses and maximum velocities are greater on the side with the thick plate. Figure 2.6a shows the effect of varying viscosity only, for models with D1= 10 km and D2= 25 km (Dr = 2.5). Asymmetry of the surface deformation is most pronounced in models with low viscosities. Models with high viscosities produce fairly symmetric deformation. Models with the lowest viscosity all show the same sense of asymmetry, with the highest velocities relative to a point on the fault on the thick plate side of the fault (though for t/TC < 0.1, the higher velocities occur on the thin plate side as in Figure 2.5b). Figure 2.6b shows velocity profiles relative to a point on the fault for the intermediate viscosity value (η =3.5x10 19 Pa s). Figure 2.7 illustrates that velocity profiles from models with different values of w (30, 50 and 70 km) look very similar. This finding is consistent with an earlier result from Schmalzle et al. (2006) who found little sensitivity of model results to w Figure 2.6. Effect of the substrate viscosity on models with a plate-thickness contrast. Modeled velocity profiles are (w was varied from 0 to 30 km in an shown for models with D2 = 25 km, w = 30 km, and three values of η. (a) velocity profiles in model reference frame. (b) For the model with η = 3.5 x 1019 Pa s, velocity profiles properties similar to those we present here). relative to a point on the fault. The heavy black line shows the velocity profile for a Savage-Burford style dislocation All model results discussed from here on use model with Vo = 40 mm/yr and locking depth (zL) equal to 17.5 km (the plate thickness at the fault). w = 30 km. earthquake cycle model with dimensions and 16 Figure 2.7. Effect of transition width w on models with a plate-thickness contrast. Modeled velocity profiles are shown for models with D2 = 25 km, η = 3.5 x 1019 Pa s, and three values of w. The heavy black line shows the velocity profile for a Savage-Burford style dislocation model with Vo = 40 mm/yr and locking depth (zL) equal to 17.5 km (the plate thickness at the fault). 2.3.2 Lateral contrasts in viscosity structure Our suite of layered viscosity models (with a uniform plate thickness of 10 km and a viscosity contrast across the fault) show behavior similar to the variable thickness plate models. When velocity profiles are plotted in the model‘s reference frame (with boundary velocities of +Vo/2 and -Vo/2), the higher-viscosity side shows less variability with time than the lower-viscosity side of the fault, and velocities at the fault (relative to the far field) vary with time (Figure 2.8). When the velocity profiles are plotted relative to a point on the fault (Figure 2.8b), it is apparent that the sense of asymmetry reverses with time between earthquakes, as it does for the variable plate thickness models. As expected, the asymmetry is larger for the model with the larger viscosity contrast (75AV13, panels a and b). The layered viscous models (75AV3L13 and 75AV3L13, not plotted) show similar behavior. Figure 2.9 shows modeled velocity profiles for the power-law viscoelastic models. These profiles vary less with time than profiles for all but the stiffest linear viscoelastic models, and the 17 lower-viscosity side for the asymmetric model shows more interseismic variation in velocity profiles. As with the linear models, the sense of velocity-profile asymmetry reverses between t/TC = 0.05 and t/TC = 0.35. We experienced some computational challenges with the power-law models, and show results for just one suite. (A model with uniform AT = 1037 Pan s was also run, though the results are not shown here. After the same number of cycles, it was producing more localized and stationary deformation than the model with AT = 1036 Pan s, but it did not reach a cycle-invariant state in 100 cycles.) For power-law models with smaller values of AT, instantaneous effective viscosities should be lower at all times. In this case the velocity profiles should vary more interseismically, yielding a larger inferred contrast in G but unreasonably low strain rates in the near field (and large inferred locking depths). The opposite should hold for models with larger values of AT. Figure 2.8. Velocity profiles for models with a viscosity contrast. (a) velocity profiles in model reference frame and (b) relative to a point on the fault for model 75AV13. (c) velocity profiles in model reference frame and (d) relative to a point on the fault for model 75AV23. The heavy black line shows the velocity profile for a Savage-Burford style dislocation model with Vo = 40 mm/yr and locking depth (zL) equal to 10 km (the plate thickness D). 18 2.3.3 Inferred elasticity contrasts from forward-modeled profiles Rybicki and Kasahara (1977) present a simple solution for surface velocities due to slip on a strike-slip fault at depth in an elastic halfspace with a contrast in shear modulus across the fault. This solution, which is a modification to the classic arctangent solution of Savage and Burford (1973), may be used to invert our forward-modeled velocity profiles at various times in the earthquake cycle for locking depth and elasticity contrast. This will indicate how a contrast in viscoelastic structure could be misinterpreted as an elasticity contrast. We perform a grid search over locking depth and elasticity contrast, computing residuals at 2-km intervals along a transect that extends 200 km on either side of the fault. Table 2.2 shows the results of this analysis. The weighted residual sum of squares (WRSS) is summed for 200 points along the profile, and we assume a one-sigma GPS velocity uncertainty of 1 mm/yr. The root mean square (RMS) misfit (in mm/yr) is the square root of 1/200 times the WRSS. We also show the ratio of WRSS to the total model variance (WRSSo) on Table 2.2 because in some cases, even for poor fits, the WRSS and RMS misfit are small because the total model variance is small. Table 2.2. Inferred locking depth and elasticity contrast (asymmetric models) Model t/TC inferred z l (km) inferred RG* WRSS/WRSSo WRSS RMS (mm/yr) .15 .45 163 1.5 0.00104 11.84 0.24 .75 260 1.5 0.00150 8.32 0.20 .15 .45 .75 19 43 1.05 1.18 .000287 .000247 16.0 9.44 0.28 0.22 .15 .45 .75 13 15 18 0.98 1.0 1.0 .000075 .000043 .000075 4.32 2.48 4.16 0.15 0.11 0.14 Models with plate thickness contrast 10 and 16 km plate CC111 CC112 CC113 19 Table 2.2. Inferred locking depth and elasticity contrast (asymmetric models), continued. Model t/TC inferred zl (km) inferred RG WRSS/WRSSo WRSS .15 .45 .75 200 300 2.5 2.5 0.00340 0.00406 35.2 22.4 0.42 0.33 .15 .45 .75 16 40 1.0 1.3 .001709 .000266 86.4 8.96 0.66 0.21 .15 .45 .75 17 19 23 0.95 0.98 1.03 .000096 .000042 .000100 5.44 2.24 5.12 0.16 0.11 0.16 .15 .45 .75 30 200 285 4 3.6 3.4 .050377 .003553 .004526 4272 37.12 26.72 4.62 0.43 0.37 CC312 .15 .45 .75 10 36 71 0.65 1.3 1.95 .020000 .000429 .000640 1664 18.72 19.04 2.88 0.31 0.31 CC313 .15 .45 .75 23 26 30 0.92 1.00 1.08 .000133 .000028 .000154 6.880 1.344 6.400 0.19 0.08 0.19 .15 .45 .75 150 240 0.24 0.26 0.0057 0.0048 99.36 37.92 0.70 0.43 .15 .45 .75 18 35 0.8 0.6 4.3e-4 0.0016 24.64 73.12 0.35 0.60 .15 .45 .75 60 150 0.28 0.26 0.0027 0.0027 111.68 44.00 0.74 0.47 .15 .45 .75 65 145 0.38 0.55 0.0027 0.0048 96.00 66.72 0.69 0.58 0.05 0.35 0.75 12 34 50 1.05 0.75 0.75 0.0087 4.3e-04 3.0e-04 584.0 19.04 10.72 1.71 0.31 0.23 10 and 25 km plate CC 211 CC212 CC213 10 and 40 km plate CC311 RMS (mm/yr) Models with viscosity contrast (D1 = D2 = 10 km) 75AV13 75AV23 75AV3L13 75AV3L12 Power-Law Model 100NAV67 ______________________________________________________________________________ *RG is the inferred shear modulus ratio. 20 2.3.3.1 Plate thickness contrast models For the plate thickness contrast models with a low substrate viscosity (3.5x1018 Pa s), we find that early in the interseismic interval, velocity profiles cannot be fit to arctangent functions (so D cannot be estimated very well). Later in the cycle, the ratio of strain rate on opposite sides of the fault is equal to the plate thickness ratio, as suggested by Chery (2008). Inferred locking depths are unreasonably large (> 100 km) for large t/TC, and the inferred elasticity ratio is comparable to the plate thickness ratio. For all of the low-viscosity models, RMS misfit for the best-fitting elastic model decreases with time in the earthquake t/Tc (Table 2.2). At the intermediate viscosity value (3.5x1019 Pa s), inferred locking depth increases substantially with time during the interseismic interval, and the sense of asymmetry reverses. Early in the interseismic interval, velocities relative to a point on the fault are higher on the thick-plate side, giving a locally higher strain rate, so it appears to have a lower rigidity than the thin side. Fit to the elastic solution is poor early in the cycle, so this effect is only shown on Table 2.2 for the model with a factor-of-4 increase in plate thickness (CC312). Later on, the thick plate side has the lower velocities (and local strain rate), giving it the appearance of having a higher rigidity. The inferred elasticity contrasts increase with contrast in plate thickness, to a factor of two late in the earthquake cycle for model CC312. The locking depth increases to over 35 km by 75% of the way through the interseismic interval. As with the low-viscosity models, RMS misfit tends to decrease with increasing t/TC (Table 2.2). At the highest viscosity value (3.5x1020 Pa s), inferred locking depths are fairly stationary and are small enough to be consistent with values inferred from geodetic studies of natural faults (typically 10 to 25 km). Inferred elasticity contrasts from the forward-modeled velocity profiles are well below 10%, even for the model with the most extreme plate thickness contrast. The sense of inferred elasticity asymmetry reverses during the interseismic interval (as described above) for all three of the high-viscosity models. 21 2.3.3.2 Viscosity contrast models Figure 2.8b illustrates that for models with a contrast in effective viscosity across the fault, the sense of asymmetry reverses. Table 2.2 does not show this because at t/TC = 0.1, the velocity profile cannot be fit to an elastic model. As with the plate-thickness contrast models, the largest inferred elasticity contrasts are associated with very large locking depths (e.g., model 75AV13, with a factor-of-100 viscosity contrast, shows a factor-of-four contrast in inferred G but the inferred locking depth exceeds 100 km). Model 75AV23, with a factor-of-10 viscosity contrast, yields more reasonable locking depths and a more modest inferred elasticity contrasts (20 to 40%, Table 2.2). The RMS misfit values decrease with t/TC, but are generally larger than for the plate thickness contrast models. The power-law model gives apparent R values of 1.05, 0.75, and 0.75 at t/TC = 0.05, 0.35, and 0.75, respectively (Figure 2.9 and Table 2.2). The inferred locking depths for the same epochs are 12, 34, and 50 km. Late in the cycle these inferred locking depths exceed observations typical of major strike-slip faults, so these deviations from lateral homogeneity are upper bounds (i.e., effective viscosity contrasts across natural faults are likely smaller). Figure 2.9. Velocity profiles for nonlinear viscoelastic models. (a) velocity profiles in model reference frame and (b) relative to a point on the fault. The heavy black line shows the velocity profile for a Savage-Burford style dislocation model with Vo = 40 mm/yr and locking depth (zL) equal to 10 km (the plate thickness D). 22 2.3.4 Locally thinned plate models Figure 2.10 shows modeled velocity profiles for thinned-plate models with (a) low and (b) high substrate viscosity values, compared with reference models assuming the same viscosities and a uniform-thickness plate. In the reference models, the plate thickness is equal to that of the thinnest part of the plate in the thinned-plate models. For the low viscosity case, the thinnedplate models yield slightly more localized deformation than the reference models. For the model with Savage parameter TC/TM = 40 (η = 3.5x1018 Pa s), late in the earthquake cycle, strain rates at the fault are almost three times their values at about 100 km from the fault (Figure 2.10c). Though this ratio is consistent with the findings of Chery (2008) that the shear strain rate should be proportional to the plate thickness (see discussion), the absolute values of these strain rates during most of the interseismic interval are very low (0.01 to 0.03 microstrains per year). Furthermore, the ratio depends on the distance from the fault at which we define the ―far field‖ strain rate: the strain rate varies with distance from the fault (approaching zero as x approaches infinity in the analytical solution). For models with a higher substrate viscosity (Figure 2.10b) the thinned-plate model yields velocity and strain rate profiles which are nearly stationary throughout the cycle and similar to that for a buried, vertical shear dislocation creeping at a constant rate (Savage and Burford, 1973). Table 2.3 shows how localized plate thinning affects inferred locking depth and slip rate. For higher substrate viscosities, models with a thinned plate yield more localized deformation and hence smaller locking depths late in the cycle, than the layered reference models. If slip rate is allowed to vary, velocity profiles from the thinned plate models with moderate to high viscosities (i.e. 40C20 and 75C21) give a slightly lower estimate of slip rate than the reference models (i.e. CC001 and CC003). For models with a high substrate viscosity, the locally thinned plate model gives somewhat higher inferred locking depths than the reference (thin plate) model, particularly, later in the earthquake cycle (at t/TC = 0.75, inferred D values are 13 km and 10 km, respectively; Table 2.3). For models with a lower substrate viscosity, the locally thinned plate has a greater effect, but the inferred slip rates and locking depths are lower than for reference model CC001 (Table 2.3). For these models, however, inferred locking depths are unreasonably large (i.e. over 10 times D for most of the interseismic interval). 23 To summarize, in thinned-plate models where the effective substrate viscosity is low enough to result in minimal shear tractions on the base of the plate during most of the interseismic interval, shear strain rates are well below values typical of strike-slip faults with a similar slip rate (for example, about 0.5 microstrains per year at the San Andreas Fault, e.g., Thatcher, 1990). Furthermore, these models predict a large variation in the pattern of surface deformation throughout the interseismic interval, which is not seen around major strike-slip faults. Figure 2.10. Velocity and strain rate profiles for thinned-plate models. (a and c) low substrate viscosity case (η = 3.5 x 1018 Pa s) and (b and d) high substrate viscosity case (η = 3.5 x 1020 Pa s). The heavy black line shows the velocity (a and b) or strain rate (c and d) profile for a Savage-Burford style dislocation model with Vo = 40 mm/yr and locking depth (zL) equal to 10 km (the plate thickness at the fault). 24 Table 2.3. Inferred locking depth and slip rate (symmetric models) Model t/TC inferred zl /D 0.15 0.15 0.45 0.45 0.75 0.75 1.0 4.0 13 17 25 36 0.15 0.15 0.45 0.75 0.75 inferred Vslip /Vo WRSS/WRSSo WRSS RMS (mm yr) 1.0 1.9 1.0 1.2 1.0 1.4 0.1333 0.0010 8.85e-4 4.90e-6 9.58e-4 6.52e-5 19099 144.96 12.96 0.072 5.600 0.381 9.77 0.85 0.25 0.02 0.17 0.04 1.0 1.0 1.5 3.5 3.3 1.0 1.18 1.0 1.0 0.98 0.0270 0.0059 2.11e-4 2.14e-4 1.39e-4 2364.2 518.4 12.48 9.120 5.920 3.44 1.61 0.25 0.21 0.17 0.15 0.15 0.45 0.45 0.75 0.75 1.0 1.0 1.0 1.0 1.0 0.99 1.0 1.07 1.0 1.07 1.0 1.05 9.28e-5 9.28e-5 6.14e-5 6.14e-5 1.10e-4 2.94e-5 5.92 5.92 3.84 3.84 6.56 1.76 0.17 0.17 0.14 0.14 0.18 0.09 0.05 0.87 1.13 2.13 2.33 3.13 3.53 1.0 1.15 1.0 1.04 1.0 1.05 0.0117 0.0015 9.9e-04 2.5e-04 4.18e-04 1.55e-05 864.0 109.6 45.44 11.68 15.52 0.587 2.08 0.74 0.48 0.24 0.29 0.05 Layered Models CC001 CC002 CC003 Power-Law Model 100NAV66 0.35 0.75 Models with locally thinned plate* 20C19 0.15 0.15 0.45 0.45 0.75 0.75 1 2.3 19 11 25 12 1 1.5 1 0.7 1 0.6 0.0704 3.58e-4 .0066 .0033 0.003 .0033 7997 40.64 55.68 28.00 17.12 19.20 6.32 0.45 0.53 0.37 0.29 0.31 40C20 0.15 0.15 0.45 0.75 0.75 1 1 1.3 3.5 2.4 1 1.2 1 1 0.9 .0305 0.0022 1.96e-4 .0048 9.11e-4 2806 205.0 11.84 192.6 36.80 3.75 1.01 0.24 0.98 0.43 75C21 0.15 0.15 0.45 0.75 0.75 1 1.1 1.2 1.3 1.3 1 1.01 1 1 0.99 2.03e-4 8.75e-5 6.77e-05 1.74e-4 7.87e-5 12.96 5.600 4.160 10.24 4.640 0.25 0.17 0.14 0.23 0.15 *For the thinned plate models, D = 10 km is assumed when computing ZL/D,. 25 2.4 Discussion and conclusions 2.4.1 Comparison with previous studies Other studies have explored the effects of lateral contrasts in effective viscosity or plate thickness for strike-slip earthquake-cycle models. Assembled from a variety of sources, the results of these studies are consistent with what we have seen here. Malservisi et al. (2001) modeled a uniform elastic plate with a substrate viscosity contrast across the fault, and concluded that asymmetric coseismic displacements must occur to attain zero net strain of the blocks on either side of the fault over an earthquake cycle. We note that for linear and either uniform or layered elasticity, this is not possible. Our finding that the sense of asymmetry of the velocity profile relative to a point on the fault reverses interseismically reconciles symmetric coseismic deformation with asymmetric interseismic deformation. Put another way, a point on our modeled fault trace (at the surface) moves along strike interseismically relative to the far-field (Figure 2.11). Lundgren et al. (2009) also present results of an asymmetric viscoelastic earthquake-cycle model of the southern San Andreas Fault which show interseismic motion of a point on the fault trace. In their model, the plate thickness was uniform and a viscosity contrast across the fault was imposed. We see this effect in both of our categories of asymmetric earthquake cycle models. We also highlight the fact that the resulting interseismic reversal in deformation asymmetry around the fault may make interpretations of interseismic velocity data that invoke elasticity contrasts inappropriate. However, in models that localize deformation to a degree that is consistent with observations (e.g., Reilinger et al., 2006, Wright et al., 2001, Fialko, 2006, Lisowski et al., 1991, Schmalzle et al., 2006, Figure 2.1), viscosity values must be high and deformation asymmetry due to contrasts across the fault is minor. Huang and Johnson (2011) developed layered models with a uniform plate thickness to test the effect of a contrast in the elastic modulus of the plate. Their models yield fairly symmetric velocity profiles, unless very large elasticity contrasts are modeled. However, based on seismic wave velocities in the Earth and laboratory experiments, elasticity contrasts should be modest at depths exceeding a few km. Significant contrasts in G (i.e. > 10%) are probably limited to the 26 upper crust, for example, in sedimentary basins (e.g. Hager et al., 1999), where wide, deep damage zones are present (e.g., at stepovers, Finzi et al., 2010) or where different rock types are juxtaposed (e.g. Schmalzle et al. 2006). Fulton et al. (2010) discuss possible causes of larger crustal elasticity contrasts and suggest that high pore-fluid pressures may act to reduce the shear modulus of rocks in some locations. Figure 2.11. Displacement profiles for an asymmetric plate-thickness model (CC312: D1 = 10 km, D2 = 40 km, and η = 3.5 x 1019 Pa s). The velocity profiles for this model are asymmetric (see Figures 2.4a and 2.4b) but the displacement over an earthquake cycle (i.e. at t/TC = 1) is symmetric (relative to the fault and in the model reference frame). Insets show how a point on the modeled fault moves during the interseismic interval, allowing the sense of velocity asymmetry to reverse, and the displacement at t/TC = 1 to be symmetric. Huang and Johnson (2011) also study a suite of viscoelastic earthquake-cycle models incorporating a plate thickness contrast, developed using the boundary-element technique. Schmalzle et al. (2006) and Fulton et al. (2010) have also developed smaller suites of such models, specifically for the San Andreas Fault at the Carrizo Plain, using the finite-element method. In the Huang and Johnson (2011) and Fulton et al. (2010) models, the thickness contrast 27 is abrupt (i.e., w = 0). Given the apparent insensitivity of our results and those of Schmalzle et al. (2006) to variations in w, the Huang and Johnson (2011) and Fulton et al. (2010) results should be broadly comparable to ours. Like us, Huang and Johnson (2011) find that velocity profiles vary less with time on the thick-plate side than on the thin-plate side of the fault. Fulton et al. (2010) and Schmalzle et al. (2006) find that the contrast in plate thickness across the Carrizo Plain is not in itself an explanation of the asymmetric velocity profile across the San Andreaas Fault in that region, and that a finite-width region of low-G rocks must be present in the upper crust northeast of the San Andreas Fault. Chery (2008) presented a model of shear deformation for an elastic plate of variable thickness, and no basal tractions, and he showed that for this model, the strain rate is proportional to the plate thickness. Chery then suggested that GPS profiles of fault-parallel velocities across faults tell us about lateral variations in the elastic plate thickness but nothing about fault slip rates. In the Earth, shear tractions at the base of the elastic plate are not zero and, even if they are quite small, they vary interseismically. For models with very low substrate viscosities (which would most closely approximate a shear traction-free condition on the base of the plate), enormous postseismic velocity transients will occur. Interseismic velocities and strain rates for most of the cycle for such models are quite small (e.g. Figure 2.10c). Hence, the Chery (2008) model does not apply to interseismic deformation around natural faults, though for large regions and long timescales, time-averaged strain rates should correlate with plate thickness as he suggests. Our models with a substrate viscosity of 3.5 x 1019 Pa s or less (20C19 and 40C20) suggest that a locally thinned plate could bias inferred locking depth downward by amplifying strain rates where the plate is thinnest. However, in contrast with Chery (2008), we find that surface strain rate does not directly correlate with plate thickness. In these models, earthquake cycle effects swamp the effect of steady-state elastic stress and strain rate concentration through variations in plate thickness, and late in the cycle, strain rates at the fault are very low (Figure 2.10c). As shown by Savage and Prescott (1978), models with very high substrate viscosities should approach the Savage and Burford (1973) solution, and due to large basal tractions are not analogous to those of Chery (2008). 28 Cohen and Kramer (1984) developed a set of thinned plate models with a viscous channel occupying a finite-width notch in the plate, rather than a continuous plate thickness variation as modeled here. In those models, sensitivity of modeled strain rates to surface velocity profiles was systematically evaluated for a range of channel widths and depth intervals, for Savage parameter (TC/TM) = 5. For models in which the depth interval of the channel was fixed and the channel width exceeded the plate thickness, near-field results were insensitive to the channel width and essentially equivalent to viscoelastic halfspace models. 2.4.2 Ramifications: Likely causes of deformation asymmetry Models incorporating moderate substrate viscosities and a large contrast in effective viscosity across the fault suggest surface deformation asymmetry which should be easily detectable with GPS and inSAR. However, these models also produce large apparent locking depths (i.e., hundreds of km), and dramatic interseismic changes to both apparent locking depth and the sense of asymmetry. Marked variations in deformation throughout the interseismic interval are not seen around major strike-slip faults (for example, the San Andreas and the North Anatolian faults, Thatcher, 1983; Ayhan et al., 2002). Localized strain around faults is often noted late in the earthquake cycle (e.g., Fialko, 2006). This also suggests fairly stationary deformation through most of the interseismic interval: elastic rebound theory requires that the integrated interseismic strain plus the coseismic strain must equal zero (on average, or over several cycles). Time-dependent earthquake cycle effects, though likely present, appear to have a limited influence on the interseismic surface velocity field. Our models with η = 3.5x1020 Pa s below the plate produce fairly stationary and localized deformation, but show very little deformation rate asymmetry. Hence a large asymmetry in surface deformation rate is unlikely to be due to a contrast in plate thickness or effective viscosity. Seismic studies of continental crust also show that lateral elasticity contrasts at length scales comparable to the thickness of the upper crust are limited to a few tens of percent or less. This suggests that geometric effects (a creeping fault zone at depth which is offset from surface trace or a dipping fault) or local elasticity contrasts (e.g. a damage zone) are more likely explanations for velocity asymmetry around major strike-slip faults (e.g. Fialko, 2006; Jolivet et al., 2008; and 29 Finzi, 2009). For extreme asymmetry, a geometrical explanation (dipping fault zone or shear zone offset at depth) is required. On the other hand, where observations of tremor in the lower crust outline a deep extension of the fault zone directly below the surface trace (Shelley, 2010), deformation asymmetry must be explained in terms of a lateral contrast in elasticity, perhaps resulting from asymmetric damage, deep sediments on one side of the fault, or high pore-fluid pressures (Fulton et al., 2010). We note that taken together, asymmetric and localized deformation throughout the interseismic interval are consistent with a thick-lithosphere model of the earthquake cycle. 2.5 Data and resources No data were used in this paper, with the exception of the GPS velocities shown on Figure 1. The North Anatolian Fault Zone GPS velocity data are from Reilinger et al. (2006). The Carrizo Plain GPS velocities were obtained on November 12, 2009 from the following URL: http://ncwebmenlo.wr.usgs.gov/research/deformation/gps/auto/CentralCalifornia/CentralCalifornia.cleaned/C entralCalifornia 30 Chapter 3 Earthquake cycle models of the Pacific-North America plate boundary at Point Reyes, California 3.1 Introduction North of the San Francisco Bay, relative motion between the Pacific and North American plates is accommodated by several sub-parallel, right-lateral strike-slip faults, all located within about 80 km of the Pacific coast (Figure 3.1). From the coast inland, these are the San Andreas Fault, the Rodgers Creek Fault, the Napa Fault, and the Green Valley Fault. Together, these faults accommodate about 36 mm/yr of relative plate motion (e. g., Prescott and Yu, 1986). I have developed a suite of earthquake-cycle models to investigate locking depths and the viscosity of the lower lithosphere in this region. Using these models I also assess whether current ideas about faults (slip rates, dips, and recurrence intervals) are consistent with GPS velocity data from the region. 3.2 Data I use three main categories of data to develop and constrain models of earthquake-cycle deformation in this region. The first includes information on slip rate, interseismic interval, and timing of earthquakes for each fault, based principally on geologic studies (though geodetic slip rate estimates exist for all four faults). The second data category relates to material properties of the subsurface and fault geometry. We base these on the US Geological Survey (USGS) 3D seismic velocity model for the San Francisco Bay Area (Aagaard et al., 2010). GPS measurements of current surface velocities are used to calibrate the model. We also make use of a catalogue of precisely located earthquake hypocenters (Waldhauser and Schaff, 2008), to investigate whether fault dips might differ from those in the USGS block model, potentially affecting inferences of fault slip rates from the GPS data. 31 Figure 3.1 Study area, showing active faults red indicates faults with historic ruptures, orange indicates faults with Holocene slip, green indicates faults with late Quaternary slip, and purple indicates undifferentiated Quaternary faults. The yellow star shows Vedanta site. The black line indicates the model profile line. Faults in the figure are labeled by their abbreviation as: San Andreas Fault (SAF), Rodgers Creek Fault (RCF), Napa Fault (NF) and Green Valley Fault (GVF). 3.2.1 Slip rates, creep, interseismic intervals, and timing of earthquakes Historic and geologic information provide some constraints on slip rate, recurrence interval, and timing of the most recent large earthquake (MRE) for each fault. Table 3.1 shows these parameters for the San Andreas Fault, the Rodgers Creek Fault, and the Green Valley Fault. These parameter values are used for most of the models, which incorporate three faults. Models with four faults also include the Napa fault. In these models, the fault parameters for the San Andreas Fault and Rodgers Creek Fault Rodgers Creek Fault given in Table 3.1 are used. The recurrence interval and timing of the most recent earthquake for the Napa Fault and the Green Valley Fault are assumed to be the same, and the summed slip rate for both faults is 9 mm/yr. I investigate the partitioning of this 9 mm/yr between the two faults in this study. 32 (a) San Andreas Fault The last major earthquake on the Marin segment of the San Andreas Fault was in 1906, and prior to that, in the mid 1600‘s (Goldfinger et al., 2003 and 2007). The Working Group on California Earthquake Probabilities, or WGCEP (2007) gives a preferred recurrence interval of 248 years for the northern San Andreas Fault at the Vedanta site (LAT 38.032 LON -122.789; Figure 3.1) based on event ages from Zhang et al. (2006). For the North Coast section (in the vicinity of the Vedanta site), the preferred recurrence interval is 288 years (Kelson, 2006), and based on Noyo Canyon turbidite data (Kelson et al., 2006) the preferred recurrence interval is 199 years. This is the same as the ~200-year interval estimated by Goldfinger et al. (2007). Slip in the most recent large event is about 5.3 +- 0.3 m (WGCEP, 2007; based on an average of slip for this segment from Thatcher [1997]). Goldfinger et al. (2007) suggest that at least 8 of the last 10 turbidite-generating earthquakes ruptured the 320 km distance from the Mendocino triple junction to San Francisco. Based on a coseismic rupture depth of 10 km (Thatcher et al. 1997) and scaling relationships given by WGCEP (2007), this suggests a M 7.6 or greater earthquake, with an average slip of at least 3 m (based on Hanks and Kanamori, 1979). The geologically-determined slip rate on the North Coast segment of the San Andreas Fault, based on Niemi and Hall (1992) and Prentice, et al. (1991), is 24 +- 3 mm/yr (WGCEP, 2007). This rate drops to 17 mm/yr for the Peninsula segment (i.e. at San Francisco, just south of the modeled profile). The San Andreas Fault slip rate for the Marin segment (i.e. the southernmost part of what WGCEP calls the North Coast segment, and the specific location of our profile) slips at a geodetically-determined rate of about 20 mm/yr (d‘Alessio et al., 2005). To match this slip rate, we model earthquakes every 250 years with 5 m of slip. Given that the last large earthquake was in 1906, the ratio of time since the last earthquake to recurrence interval (t/Tc) was 104/250, or 0.42, in 2010. A range from 200 to 288 years is reported for mean recurrence interval, so t/TC may range from 0.36 to 0.52. The Marin segment of the San Andreas Fault does not appear to creep at the surface (WGCEP 2007). Aseismic slip rates inferred at the surface are less than 0.4 mm/yr (Galehouse & Lienkaemper, 2003, Burford & Harsh, 1980, and other references cited by WGCEP (2007). 33 (b) Green Valley Fault WGNCEP (1996) provides a preferred recurrence interval of 330 years for the Green Valley Fault. More recently, the recurrence interval for the Green Valley Fault has been estimated at 220 plus or minus 160 years, with the most recent event dated at AD 1600 +- 130 years (Lienkaemper et al., 2008). These estimates are from studies of a single trench. Another trenching study yields an interevent time of 865 years, based on dates for three earthquakes (Baldwin et al., 2008) but the authors report a 200-year interevent time from another trench (Mason Road) and suggest that a combined analysis of both sites could yield a better estimate of interseismic interval for the Green Valley Fault. Baldwin et al. (2008) also find that the most recent surfacerupturing earthquake occurred between 1573 and 1799 (95% confidence range). We use the more recent recurrence interval estimate of 220 years in our modeling, and assume that we are currently at the end of that interval. According to D‘Alessio (2005), the geodetically derived slip rate on the Green Valley Fault is 9 mm/yr. To achieve this slip rate, we model earthquakes every 220 years with 2.0 m of slip. For the purpose of modeling, we assume that the last earthquake was 220 years ago and that t/T is just under 1. The minimum late Holocene dextral slip rate is 3.8 mm/yr to 4.8 mm/yr (Baldwin and Lienkaemper, 1999. Galehouse & Lienkaemper (2003) give a creep rate of 4 mm/yr at the surface. According to J. Lienkamper (pers comm., 2010) this creep rate is based on 25 years of observations at one site and a shorter timespan of observations at several other sites (c) Rodgers Creek Fault Paleoseismic observations on the Rodgers Creek Fault (Budding and others, 1991; Schwartz et al., 1993) show three surface-rupturing earthquakes between about AD 1000 and 1776. These three events produced 5.1 to 7.2 m of offset, with slip during the most recent event of 1.8 to 2.3 meters. The recurrence interval is estimated as 131 to 370 years with a preferred value of 230 years (Schwartz et al., 1993). The most recent large earthquake on the Rodgers Creek Fault occurred between 1690 and 1776, most likely between 1715 and 1776 (Hecker et al., 2005). WGCEP (2007, Appendix B) estimates a preferred recurrence time of 305 years (with a minimum of 220 years and maximum of 390 years), based on data from Budding et al. (1991) and the most recent earthquake date from Hecker et al. (2005). WGCEP (2007) also estimates a 34 geologic slip rate of 9 2 mm/yr for the Rodgers Creek Fault, which is an average of values from Schwartz, et al. (1992) and a slip rate from Hayward fault (Lienkaemper and Borchardt, 1996). To be consistent with the geodetically-determined slip rate (7 mm/yr, D‘Alessio et al., 2005), and the most recent geologic estimate of recurrence interval, we model 2.1 m of slip every 300 years. The most recent quake is assumed to have been in 1745, putting us 265 years into the earthquake cycle (t/TC = 0.88). For the shorter recurrence times (e.g., from Schwartz et al., 1993), t/TC is greater than 1, and the minimum possible value based on a 390-year recurrence interval estimate is about 0.8. The Rodgers Creek Fault creeps at the surface at a rate of about 3 mm/yr (Hecker, pers. comm., 2010. According to Galehouse & Lienkaemper (2003) the rate is 0.4 to 1.6 mm/yr and this is also cited by WGCEP (2007). Table 3.1. Slip rate and earthquake timing values used for modeling of earthquakes along the San Andreas Fault, the Rodgers Creek Fault and the Green Valley Fault. (d) Time since Last quake Fault Slip rate (mm/yr) Period (Yr) San Andreas 20 250 104 Rogers Creek 7 300 265 Green Valley 9 220 220 (Yr) Napa Fault At first our FE model included just the San Andreas Fault, the Rodgers Creek Fault and the Green Valley Fault. The Napa Fault, between the Rodgers Creek Fault and the Green Valley Fault, was assumed to be inactive, as no geologically-determined, Holocene slip rates were available for this fault (Bryant, 2000). The Napa Fault was added to the model after I had trouble fitting the GPS velocity profile with just three faults. This addition was supported by recent publications suggesting that Calaveras fault slip (6 mm/yr; Simpson et al., 1999) is partitioned into both the Napa Fault and Green Valley Fault. Evidence from detailed geologic mapping studies (Brossy et al., 2010; Wesling and Hanson, 2008) and block models (d‘Alessio et al., 2005) support an active Napa fault. Brossy et al. (2010) argue that the slip from the northern 35 Calaveras fault transfers north via the Contra Costa Shear Zone to the Napa Fault, rather than to the Concord Fault, as had previously been assumed. Based on detailed mapping, and on observations from a trench at the Napa Airport and one other surface exposure of the Napa Fault, Wesling and Hansen (2008) argue for Holocene activity of the Napa Fault, though they cannot provide a slip rate or event dates (other than ―several late Quaternary events‖ in the trench). Kinematic block models constrained to GPS data (D‘Alessio et al., 2005) also suggest that slip is partitioned between the Napa Fault and the Green Valley Fault; their models estimate slip rate of 3-4 mm/yr for the Napa Fault. In models with four faults, I model the Napa Fault using the same earthquake timing as the Green Valley Fault. Slips per event for both faults are set to partition a total of 9 mm/yr of slip (i.e. the previously inferred Green Valley Fault rate) between them. 3.2.2 Material properties and fault geometry Distinct lithologic units within the modeled profile are based on the USGS 3D geological model (Aagaard et al., 2010), and are shown on figure 3.2. For each lithologic unit, I have estimated the mean P-wave velocity using values from the USGS Bay Area Velocity model, version 08.3.0 (Aagaard et. al. 2010). This model provides a three-dimensional view of the geologic structure and physical properties of the region down to the depth of 45 km. The database of the velocity model is located in two files (―USGSBayAreaVM-08.3.0.etree.gz‖ and ―USGSBayAreaVMExt-08.3.0.etree.gz‖) downloaded (http://earthquake.usgs.gov/regional/nca/3Dgeologic/). from the USGS model website Model query software is required to query the seismic velocity model for physical properties. This software provides P wave velocity and density of layers in a particular group of cells. The average of the values for cells comparison a block is used to compute the elastic moduli of that block. Figure 3.2 shows the value of density in kg/m3 and P wave velocity in m/s for the blocks. Note that depths to the upper crust - lower crust boundary and the Moho vary along the profile. To get the elastic modulus of rock units, P wave velocity and density are used to compute elastic moduli for cells by the following relation. 36 vp (k 43 G) v s G (3.1) Where p and s are respectively P and S wave velocity, is density, G and k are respectively shear and bulk moduli. (Poissons ratio) is obtained from: 2 1 2 v s v p 2 v 2 2 s v p (3.2) Whenever we assume Poisson ratio equal to 0.25 (a typical value for the crust) we will have, v p 3 vs (3.3) and G 1 2 3 vp , k 5 2 9 vp (3.4) Fault geometries are based on the USGS 3D geological model (Figure 3.2, Aagaard et al., 2010; http://earthquake.usgs.gov/regional/nca/3Dgeologic/). According to this source, the Napa Fault is a minor feature that does not extend completely through the crust, and the San Andreas Fault and Green Valley Fault are vertical. The Rodgers Creek Fault dips 80 degrees to the east. The San Andreas Fault is the main boundary between the Pacific and North America plates. The upper crust at the west side of San Andreas Fault is granite (the Salinian block) and at the east side is Franciscan, metamorphosed seafloor rocks. To determine the location of the faults on the profile and their distance from San Andreas Fault, I use the fault activity map of California and the adjacent area (Jennings, 1994) which does not shows the Napa Fault. In my revised model with four faults, I use the new version of the Jennings map (2010) which show the Napa Fault. 37 Figure 3.2. Rock types and material property of the blocks in the central part of the model. Density of blocks is in kg/m3 and p-wave velocity is in m/s. SAF stands for San Andreas Fault, RGF stands for Rodgers Creek Fault, NF stands for Napa Fault and GVF stands for Green Valley Fault. 3.2.3 GPS data Modeled interseismic velocities at points on the surface are compared with data from GPS networks in the area. Of 158 continuously operating GPS sites (Table 3.3) in northern California, 14 are within 45 km of the profile (Figure 3.3). GPS velocity data from continuous GPS sites for the San Francisco Bay Area are available to the public via the USGS website: http://earthquake.wr. usgs.gov/monitoring/gps/data/SFBayArea/cleaned_velocity_file. 38 Figure 3.3. GPS velocities north of the San Francisco Bay at sites within 45 km of the modeled profile. The San Andreas Fault is indicated in red and the black line shows the modeled profile. Velocities are relative to North America reference frame. In addition to the continuous GPS site velocity data, I use data from campaign-mode GPS sites. From 24 GPS campaign sites in the profile area (Table 3.4), 21 are within 45 km of the profile line shown on Figure 1. The North Bay campaign data have been folded into a larger set of campaign stations (142 stations). Velocity data from these sites are available from the following USGS website: http://earthquake.wr.usgs.gov/monitoring/gps/data/Pillsbury/cleaned velocity file The latitude-longitude coordinates of GPS sites along the modeled profile are transferred to a local Cartesian coordinate system with its origin located at the intersection of the San Andreas Fault and the profile line. The local GPS site coordinate and GPS velocity vectors are rotated 125.5 degrees to obtain respectively the sites‘ distance from the San Andreas Fault and the GPS velocity resolved in the San Andreas Fault strike-parallel direction. To compute the variance of the GPS velocity in the strike-parallel direction, the following relation is used (Schofield and Breach, 2007). 39 a2 E2 sin 2 (a) 2 EN sin(a) cos(a) N2 cos 2 (a) (3.5) Where is the strike azimuth (positive clockwise angle from the north), E is easting error, N is northing error, and EN is the covariance of easting and northing. Figure 3.4 shows the San Andreas Fault strike-parallel velocity components for the GPS sites, with their uncertainties, along the profile line. Blue points indicate continuous and green points indicate campaign GPS sites. The locations of the faults on the profile line are also shown by cyan line. These data are summarized in Tables 3.3 and 3.4. Figure 3.4. Profile of plate-boundary parallel GPS site velocities in a stable North America reference frame. Error bars show 95% confidence limits. 3.2.4 Seismicity data To compare the models with another source of data and more precisely determine the faults‘ location and geometry, I used high-precision hypocenter locations for 445,493 earthquakes that were recorded between January 1984 and December 2008 by stations of the Northern California Seismic Network (NCSN) (Waldhauser and Schaff, 2008 and Schaff and Waldhauser, 2005). 40 1813 epicenters from this catalogue are within 45 km of our profile (Figure 3.5) and are used in this study. Figure 3.5. Earthquake epicenters within 45 km of the model profile. 1813 earthquakes were recorded in this area between 1984 and 2008. All magnitudes are less than 5. Blue line indicates the San Andreas Fault. These hypocenters are inspected in plan view to and projected onto the profile plane to determine fault dips and surface coordinates. Figure 3.5 shows the epicenters in plan view. My analysis of fault geometries based on these hypocenters is deferred until the Discussion section because the initial model design (and most of our models) do not make use of the refined revised fault geometry based on these hypocenters. 3.3 Method To investigate plate-boundary deformation along this profile, I have developed 2D profile viscoelastic finite-element models (Appendix A). These models incorporate an elastic upper crust underlain by a heterogeneous viscoelastic half space representing the lower crust and 41 mantle. The locked zone of the faults extends from the surface to the base of the upper crust. Viscoelastic shear zones extend each fault zone through the lower crust and upper mantle (Figure 3.6 A). To avoid the effect of mesh boundary conditions on the modeled velocities, the model mesh extends to 500 km west and east of the San Andreas Fault, and 500 km in depth. The model is discretized to very fine elements, less than 2 km, along the faults in the upper crust. Toward the boundaries of the model, elements size is increased gradually to a maximum dimension of 50 km. Figure 3.6 (B) shows the mesh of the model around the faults. These models are discretized to linear hexahedron elements. Each model has about 10000 elements. The top and bottom surfaces of the model mesh are stress-free boundaries. The front and back faces of the model mesh are free to displace in the fault-parallel direction. A relative velocity of 36 mm/yr is imposed using constant velocity boundary conditions at the east and west sides of the model mesh, in the upper crust. Below this, the side boundaries are traction-free. The "split node" technique (Melosh and Raefsky, 1981) is used to represent the faults, and the depth of the interseismically locked interval is varied. For most simulations, the depth of the interseismically locked (and coseismically slipping) interval is fixed at 12 km. A viscous shear zone extends each fault zone downward into the mantle. These shear zones are 4 km (2 elements) wide, and extend to a depth of 70 km (Figure 3.6 A). The code Geophysical Finite Element Simulation Tool (GeoFEST; Parker et al., 2003) is used for this modeling. The integrated slip rate at this part of the Pacific-North America plate boundary is 36 to 38 mm/yr (d’Alessio et al., 2005), which is partitioned between the Green Valley Fault, the San Andreas Fault, and the Rodgers Creek Fault. Three faults with different earthquake timing schedules (different MRE dates and recurrence intervals) result in a system-wide periodicity of 16500 year. For models incorporating the Napa Fault, I assume that the earthquake period and last quake are the same as for the Green Valley Fault, so the system-wide period does not change. WRSS d i mi 2 i2 (1.1) where d is data, m is value of model and d is the standard deviation of data. Values belonging to individual GPS sites (site number i) are summed. 42 Figure 3.6. A) Model schematic showing three faults, and shear zones extending the faults into the mantle. Dark red surfaces are coseismic rupture planes which are interseismically locked. The boundaries between the upper crust, and the lower crust and the Moho, are shown. B) Model mesh in the vicinity of the faults. 43 Smaller WRSS misfit means a more robust model. The value is strongly dependent on the data set; by increasing the number of data points or decreasing the standard deviation of them the value of WRSS increases, so that for the same model with different data set WRSS can change dramatically. The data set to be used in this study is from 30 GPS sites on the profile band (Tables 3.3 and 3.4).To recognize how WRSS is sensitive with the used data set; WRSS for two models are compared. First one is a dislocation model for four parallel faults at the data region; WRSS for this model is 1155. The second model is same as first one with replaced fault 5 km to oceanic side; WRSS for the model is 1597. For comparison the summed variance, (minimum WRSS for a model with no change in velocity across the profile) is 48290. The GPS velocity data are in a North America reference frame, but the model is in a local coordinate system. To compare data and a model, I need to transfer data to the local coordinate system which the model is in. In a comparison of data and model the most important issue is relative velocity of the points. Because all velocity of GPS sites have uncertainties, we can‘t take any of them as a reference point and compare other points that one, unless we assume one of the sites with a small velocity error is the reference point (and it has a negligible velocity error). The most reasonable method is that we find a local datum for which the weighted residual sum of square is a minimum. So we need to shift the velocities by some constant offset, T, to place them in the new datum. To find the value of T we should minimize the WRSS which is computed in new datum: WRSS d i T mi 2 i2 (1.2) The derivative of Equation (1.3) with respect to T should be set to zero, to find the T for the minimum WRSS. T d i mi i2 1 (1.3) 2 j 44 3.4 Results 3.4.1 Locking depth The first parameter I vary in this study is locking depth. Locking depth defines the depth interval on a fault which slips in characteristic earthquake and is ―locked‖ interseismically, and provides an upper limit on seismic moment (and hence moment magnitude). Seismic moment is shear modulus times rupture area times fault offset. Since the rupture surface area depends on the locking depth interval, and rupture area also correlates with offset, locking depth provides constraints on maximum earthquake magnitude. Here, I run models in which locking depth is varied from 5 km to the depth of upper crust at each fault (up to 25 km). For all these models, viscosity of the lower crust and mantle is assumed to be 5x1019 Pa s and viscosity of the shear zones extending the faults through the lower crust is 1x1019 Pa s. Figure 3.5 shows results from these models. Models with a 25km locking depth have the maximum misfit, 2492. WRSS misfit is minimized when locking depth is 13 km. This result is consistent with maximum seismicity depths in the region, as well as the results of elastic block models (D‘Alessio et al., 2005). Figure 3.7. WRSS of the model with respect to locking depth. Mantle viscosity is 5x1019 Pa s and viscosity of the shear zone is 1019 Pa s. 45 Although misfit is minimized at 13 km, it is insensitive to variation in locking depth between 10 and 13 km. Models with a locking depth of 12 km have practically the same misfit as those with a locking depth of 13 km. Denser and more precise GPS data could improve this estimate only somewhat; the elastic upper crust tends to filter short-wavelength strain rate features at the surface. 3.4.2 Viscosity of shear zone Viscous shear zones, where reduced grain size and elevated volatile content may reduce the effective viscosity relative to the host rock, extend fault zones downward into the crust. The depths to which these shear zones extend are uncertain and are difficult to resolve with surface deformation data. To place limits on acceptable shear zone viscosity values, I have run several sets of models with different mantle viscosity values. First of all it is obvious that a shear zone has a lower viscosity than its surroundings, so in all studies the maximum viscosity for the shear zone is equal to the minimum viscosity of adjacent rock. Generally, lower viscosity improves the fit to the data and there is no minimum for the shear zone viscosity–misfit curve (Figure 3.6). Practically any shear zone viscosity is acceptable when the contrast in viscosity with the host rock is over about 100. Figure 3.6 shows misfit plotted against shear zone viscosity for a model with a 10 km locking depth and 5x1019 Pa s lower crust and upper mantle viscosity. A viscosity of 1017 Pa s is acceptable for the shear zone since misfit does not change meaningfully for lower viscosity. Furthermore, a high contrast in the viscosity of adjacent zones dramatically slows the simulation and sometimes causes numerical problems. Similar results are obtained for models with a higher viscosity for the lower crust and mantle (1021 Pa s) and more complicated models with different viscosities for the lower crust and upper mantle (1021 and 5x1019). 46 Figure 3.8. WRSS of models with different shear zone viscosities. A 10 km locking depth and 5x1019 Pa s viscosity of lower crust and mantle are assumed. 3.4.3 Viscosity of lower crust and upper mantle To calibrate viscosity of the lower crust and mantle and to check the probability of a contrast in viscosity across the San Andreas Fault, several models with different viscosities from 1018 to 5x1021 Pa s are developed. Misfit for the lower viscosity value is very high. The minimum misfit is for model with 1020 Pa s west of the San Andreas Fault and 1021 Pa s for the east. I find that 47 lower crust and mantle in the area of study have a high viscosity. After additional models were run with moderately high viscosities, I found that viscosity values of 1 to 3x1020 Pa s for the west and 5x1020 to 1021 Pa s for the east side of the San Andreas Fault were best. Figure 3.9. Contour plot of WRSS for different viscosities of the lower crust and mantle. Viscosities east and west of the San Andreas Fault are varied independently. Blue points indicate individual models. Green line indicates laterally homogeneous models (no viscosity contrast). For all the models in this figure (/w) viscosity of the shear zone is 1017 Pa s/m and locking depth is 12 km. The minimum misfit region is offset somewhat from this line, illustrating that the effective viscosity of the lower crust and upper mantle west of the San Andreas Fault is perhaps five times smaller than to the east. Uniform models with an effective viscosity of about 1020 to 1021 Pa s are also admissible. 48 3.4.4 West Napa Fault Inspection of misfit along the model profile (for all of our models) shows large residuals just west of the Green Valley Fault (Figure 3.10). There are several data points in this location which show lower velocities relative to the model. At first glance, the data seem to suggest that the Green Valley Fault is located west of its assumed (modeled) location. To test this idea, I designed a new model in which the Green Valley Fault is moved 10 km toward the San Andreas Fault. The new model reduced the WRSS by a factor of two (Figure 3.10 b). Figure 3.10. Model velocity profiles, showing sensitivity of WRSS to location of the Green Valley Fault. 49 On the other hand, the location of the Green Valley Fault‘s surface trace (relative to the San Andreas Fault) is correct. The West Napa Fault, which is 10 km west of the Green Valley Fault, is shown to be inactive on the 1994 Jennings map of active California faults, but recent block modeling and geologic studies suggest some level of activity on this structure (D‘Alessio et al., 2005, Wesling and Hanson, 2008). In particular, D‘Alessio et al. (2005) suggest that the 9 mm slip rate of Calaveras fault is partitioned between the Napa Fault and Green Valley Fault in our study area. To find how 9 mm/yr of slip is partitioned between these two faults, I ran a set of models with a range of slip rates for each, summing to 9 mm/yr. In all of these models locking depth is 12 km and beneath it, shear zones extend to 70 km depth with viscosity of Pa s. The lower crust and upper mantle viscosity west of the San Andreas Fault is 5x1020 Pa s and east of the San Andreas Fault is 1021 Pa s. Slip rates of the first model for the Napa Fault and the Green Valley Fault are respectively 0 and 9 mm/yr. On the next models slip rate gradually increases for the Napa Fault and decrease for the Green Valley Fault, and for the last model, these values are 9 and 0 mm/yr. Figure (3.11) shows the results of the models. Figure 3.11. Sensitivity of WRSS to relative slip rate for the Napa Fault and the Green Valley Fault. Note that by decreasing the Green Valley Fault and increasing the Napa Fault slip rate, misfit between data and model is reduced. 50 WRSS misfit is continuously improved by increasing slip rate on the Napa Fault and reducing it on the Green Valley Fault. It means that new model can‘t remedy the models misfit to the GPS data in the area and there must be another explanation for it. In contrast, D‘Alessio et al. (2005) find a slip rate of just 4 mm/yr along the Napa Fault and geological studies suggest very low rates (Wesling and Hanson, 2008). 3.5 Discussion Modeled surface velocities are highly sensitive to mantle and lower crust viscosity, and models with high viscosity values provide the best fit to the GPS velocity profile. The models and GPS velocities also suggest a viscosity contrast across the San Andreas Fault, with a somewhat lower viscosity west of the San Andreas Fault, but this finding is less robust. The GPS velocity profiles show localized strain around the faults, all of which are late in their earthquake cycles, and most of the relative plate motion is accommodated within about 80 km of the plate boundary (e.g. Prescott and Yu, 1986). This is compatible with high-viscosity models (e.g., Savage and Prescott, 1978). In the low-viscosity models, plate-boundary strain is distributed across hundreds of kilometers, so the relative motions across the profile are far too small. Boosting the fault slip rates to unreasonably high values would improve the fit of the lowviscosity models to the GPS velocity profiles, but this would violate geological constraints on fault slip rates. We did not investigate models with effective viscosity values exceeding 3 x 10 21 Pa s. For models with high effective viscosity values to reach cycle invariance, we need to run the models for a longer period of time. This requires not only longer run times, but also special codes for modeling large strains. Over long simulation times, large strains accrue in model elements near the faults (specifically, in the viscous shear zones and in the mantle just below these) and either re-meshing of the modeled domain or other techniques (e.g. an EulerianLagrangian formulation) is required to prevent inaccuracy or numerical instabilities. Figure 3.12 compares the surface velocities from the best viscoelastic finite element model with a dislocation model which has just the faults and shear zones which creep at a constant rate, embedded in an elastic half space. GPS velocities in the figure are shifted 28.74 mm/yr to 51 transform from the North America frame to the local frame. All other parameters for both of the models are same. There are 4 faults in the model: the San Andreas Fault, Rodgers Creek Fault, Napa Fault and Green Valley Fault, and their slip rates respectively are 20, 7, 3 and 6 mm/yr. Locking depth is 12 km for all of the faults. For this particular case WRSS for viscoelastic finite element model is 672 which is comparable to than 802, the WRSS of the dislocation model. Visual inspection shows that the viscoelastic layer localizes strain around the faults. Figure 3.12. A comparison of elastic dislocation and viscoelastic finite element models incorporating all four faults. In both models slip rates for the San Andreas Fault, Rodgers Creek Fault, Napa Fault and Green Valley Fault are respectively 20, 7, 3, and 6 mm/yr and locking depth is 12 km. WRSS for the finite element and dislocation models are 672 and 802, respectively. Models with different slip rates split among the Napa Fault and the Green Valley Fault show that a higher slip rate on the Napa Fault results in a smaller WRSS. However, as mentioned in the Results section, block models and geologic studies suggest more slip on the Green Valley Fault than the Napa Fault. The contradiction between previous studies and these models indicates that there is some other mechanism in the region which causes large residuals west of the San Andreas Fault. One possibility is that the modeled fault positions or dips are incorrect, or that it is inappropriate to use a 2D profile model because of changes in fault strike north or south 52 of the modeled profile. Using high-precision hypocenter locations (Waldhauser and Schaff, 2008), I investigate these possibilities. The plan view of epicenters shows that the Rodgers Creek Fault and Napa Fault are approximately parallel to the San Andreas Fault (Figure 3.13). At first view, the Green Valley Fault does not look parallel to the San Andreas Fault, but when the focal depth of the earthquakes are considered, it is obvious that south of the model profile (in Region 2 on Figure 3.13) hypocenters are generally deeper. The dip angle of the Green Valley Fault causes epicenters of deeper earthquakes south of the profile to appear to be closer to the San Andreas Fault. Although south of the profile the Green Valley Fault does deflect toward the southwest, it can be accepted as approximately parallel to the San Andreas Fault across the modeled profile. In Figure 3.13, the 90-km swath centered on the profile is divided to two zones. Zone 1 is north of the modeled profile line, and zone 2 is to the south. All hypocenter data from these two zones are projected on to the profile plane (Figure 3.14) for comparison. The distribution of hypocenters in both zones in Figure 3.13 shows that the surface trace of the Green Valley Fault in the model is in the right place with respect to the San Andreas Fault. Simple linear regressions are used to find the best fitting fault line to the hypocenter positions along the profile. For the Rodgers Creek Fault, most of the hypocenter data are from zone 1; the regression shows that the surface trace of the Rodgers Creek Fault is 35 km from the San Andreas Fault. There is dense hypocenter data for the Napa Fault in zone 1, however the result of the regression for zone 2 is compatible with the result from zone 1. The majority of the earthquake hypocenters on the Napa Fault are located at depth between 6 and 9 km, (Figure 3.13) and only a small number of them are located at shallower depths. The location of the surface trace of the Napa Fault confirms the result and compensates for the lack of shallow data. The result of the regression shows that the surface trace of the Napa Fault on the profile is 54 km east of the San Andreas Fault. This distance is 61 km in the previous models. Shifting the Napa Fault to the west improves the fit of the model to the GPS data result, but does not completely account for the misfit. For the Green Valley Fault the results of both zones indicate the fault on the surface is located 71 km from the San Andreas Fault. The study also shows that the local coordinate of GPS sites on profile direction is separated 2.6 km toward west from the local coordinate of the model. The effect of replacement on the result is negligible and WRSS misfit improves less than 0.1%. 53 Figure 3.13. Epicenters of earthquakes within 45 km of the modeled profile. Colors of points indicate the focal depth. The profile line, in the middle, divides the region to 2 zones. The numbers in the white squares indicate the zone numbers. Epicenters are from Waldhauser and Schaff (2008). The dip angles of the Napa Fault and the Green Valley Fault, which have been assumed to be vertical, may explain the GPS-model residuals in their vicinity. If these faults dip to the west and intersect the lower crust west of their surface traces, areas of high interseismic strain rates will be deflected to the west, which will improve the fit of the models to the GPS data. To assess fault dip, all the hypocenters of Zones 1 and 2 on Figure 3.13 are projected onto the profile plane. The result of linear regression shows that the Napa Fault and Green Valley Fault are not vertical, but are dipping to the west (Figure 3.14). The best-fitting line through the Napa Fault hypocenters has a dip angle of 118 degrees (that is, 62 degrees to the west). Separate regressions for the hypocenters of Green Valley Fault earthquakes in Zones 1 and 2 give, respectively, dip angles of 109 and 115 degrees (that is, 81 and 75 degrees to the west). When hypocenter data from both zones are analyzed together, the dip angle is 117 degrees (73 degrees to the west). Since the fault trace is bending (trending north rather than normal to the model profile) in Zone 2, the result for Zone 1 is more reliable. 54 Dip angles of the Napa Fault and the Green Valley Fault can probably explain the too-low modeled surface velocities and strains west of the Green Valley Fault. The pattern of interseismic strain depends on where shear zone creep and concentrated viscoelastic lower crust and mantle relaxation is occurring (e.g. Segall, 2010), and with the dipping faults, these are offset several kilometers to the west of the fault traces at the surface. Several hypocenters, mostly in Zone 2 and at a depth of 5 to 10 km, fall along a plane that intersects the surface 63 km east of the San Andreas Fault. There is small amount of hypocenter data for this fault in zone 2 which gives a similar result. If this plane is an active fault, strain associated with it may also contribute to the surface deformation west of the Green Valley Fault surface trace. This fault has a 114 degree dip angle (72 degrees to the west) and its surface trace is located at about 10 km west of Green Valley Fault. This fault, shown with a dashed line, and its corresponding hypocenters, are shown in Figure 3.14. In Table 3.2 I refer to this structure as ―Unknown Fault‖ (UF) . Figure 3.14. Earthquake hypocenters (Waldhauser and Schaff, 2008) projected on to the profile plane. Red circles indicate Zone 1 in figure 3.12 (North West of the profile line) and blue circles indicate Zone 2. Faults in the figure are labeled by their abbreviation as: San Andreas Fault (SAF), Rodgers Creek Fault (RCF), Napa Fault (NF), Unknown Fault (UF) and Green Valley Fault (GVF). Table 3.2 shows the result of the regression with hypocenter coordinates from Zone 1, Zone 2, and the whole dataset for the faults in Figure 3.13 and 3.14. For some faults, hypocenter data is 55 mainly located in one of the zones. For instance, there are few Rodgers Creek Fault hypocenters in Zone 2 and few Unknown Fault hypocenters in Zone 1, so the results from these zones for these faults are not robust. When the inferred fault dip angles for both zones are in agreement (as for Unknown Fault) the dip estimate is more robust. Table 3.2. Fault parameters (distance from the San Andreas Fault and dip angle) obtained from hypocenter data. Zones 1 and 2 are shown in Figure 3.13. Data of Zones 1 and 2 are shown in Figure 3.14 with red and blue circles respectively. Unknown Fault (UF) is also shown in Figure 3.14 (dashed line). Fault RCF 3.6 Zone 1 hypocenters Distance from Dip . angle SAF 35 km 83˚ Zone 2 hypocenters Distance from Dip angle SAF 32 km 83˚ Both zones Distance from Dip angle SAF 36 km 89˚ NF 54 km 118˚ 53 km 116˚ 54 km 118˚ UF GVF 64 km 71 km 124˚ 109˚ 62 km 71 km 114˚ 115˚ 63 km 71 km 119˚ 117˚ Conclusions I use viscoelastic earthquake cycle models to model Pacific-North America plate boundary deformation at Point Reyes, California. In this area, about 36 mm/yr of right-lateral relative motion is accommodated by the San Andreas, Rodgers Creek, Napa and Green Valley faults. The models are calibrated to plate-boundary parallel GPS site velocities, and take into account the timing and recurrence intervals of earthquakes. A 12-km locking depth is inferred from the models and is consistent with local seismic data (e.g. D‘Alessio et al, 2005). Shear zones are found to have a viscosity of 1018 Pa s or less. In order to fit the high velocity gradient across the system of faults, a high effective viscosity for the lower crust and mantle (5x1020 Pa s west of the San Andreas Fault side and 1021 Pa s to the east) is required. The data support a modest contrast in effective viscosity of the lower crust and upper mantle across the San Andreas Fault, with lower viscosity values to the west. In the region between the Rodgers Creek Fault and the Green Valley Fault, GPS data suggest higher strain rates than the models can explain. Even after adding the Napa Fault and allocating all the slip rate of Calaveras fault, 9 mm, to the Napa Fault, there is still a higher strain rate in the region than my models can reproduce. By studying local earthquake hypocenter data, I find that 56 both the Green Valley Fault and the Napa Fault dip toward the west. The hypocenter data also suggest the presence of another fault zone between the Napa Fault and the Green Valley Fault. The westward dip angles of these faults explain why high strain rates are offset to the west of surface traces of the Green Valley Fault and the Napa Fault. Future models incorporating dipping faults should be developed to explore how the slip rate of the Calaveras fault is partitioned between the Green Valley Fault, the Napa Fault, and the other fault (UF) between them. In addition, slip rates estimated from block models (particularly those with small block dimensions) may be affected if even modestly dipping faults are assumed to be vertical. Since dips of 5 to 10 degrees off vertical are often missed in seismic studies (e.g. Fuis et al., 2008), creeping shear zones at depth may be offset by several kilometers from the surface trace, leading to incorrect inferences of material asymmetry (chapter 2 and Vaghri and Hearn, 2011) and possibly, errors in the attribution of slip rates to closely spaced, active faults. 57 Table 3.3. Velocities of continues GPS sites, is latitude, is longitude, vE is easting velocity, vN is northing velocity, E is easting error, N is northing error, and EN is the covariance of easting and northing. Site vE vN E N EN 208P 39.109302 -122.303869 -9.45 6.36 0.09 0.22 -0.0333 229P 37.749435 -121.977956 -14.21 15.28 0.31 0.23 -0.0245 BRIB 37.919405 -122.152552 -14.13 16.73 0.28 0.28 -0.0038 CHAB 37.724116 -122.119310 -15.78 19.29 0.58 0.36 -0.0063 CHO1 39.432639 -121.664962 -9.46 6.98 0.06 0.09 0.0099 CMBB 38.034176 -120.386040 -9.85 7.98 0.31 0.06 -0.1080 CNDR 37.896405 -121.278495 -11.84 7.85 0.69 0.05 0.0254 DIAB 37.878575 -121.915627 -11.41 11.30 0.19 0.35 -0.0367 EBMD 37.815008 -122.283805 0.00 0.00 0.00 0.00 -0.0449 FARB 37.697208 -123.000766 -26.85 37.25 0.04 0.20 -0.0090 HOPB 38.995184 -123.074727 -18.10 21.79 0.26 0.34 0.0103 LNC1 38.846512 -121.350235 -10.12 5.69 0.44 0.36 -0.0442 LUTZ 37.286851 -121.865225 -19.22 23.59 0.21 0.19 0.0030 MHCB 37.341533 -121.642578 -11.87 11.53 0.05 0.32 0.0037 MOLA 37.946574 -122.419923 -19.08 23.01 0.48 0.39 0.0013 MONB 37.485326 -121.866865 -15.25 16.55 0.13 0.10 0.0070 MUSB 37.169941 -119.309350 -9.96 8.69 0.26 0.15 -0.0158 OHLN 38.006253 -122.272988 -14.13 17.43 0.20 0.43 -0.0017 OXMT 37.499363 -122.424315 -24.23 32.32 0.12 0.26 0.0430 P059 38.928346 -123.726197 -22.83 35.31 0.39 0.31 -0.0486 P171 36.485524 -121.792518 -28.04 35.88 0.11 0.93 -0.0328 P176 37.471773 -122.357139 -23.03 30.76 0.50 0.41 -0.0114 P177 37.528168 -122.495052 -30.36 31.47 0.61 0.53 -0.0499 P178 37.534519 -122.332363 -22.68 27.39 0.47 0.38 -0.0521 P181 37.914545 -122.376755 -18.36 24.21 0.30 0.23 -0.0230 P182 38.495015 -123.181245 -20.12 31.62 0.39 0.31 -0.0940 P183 38.313663 -123.068887 -21.33 32.96 0.36 0.28 -0.0368 P184 39.117165 -123.708937 -20.77 29.52 0.45 0.36 -0.0366 P185 39.261304 -123.749335 -19.89 28.22 0.40 0.31 -0.0350 P186 39.150179 -123.518152 -19.45 23.40 0.45 0.36 -0.0995 P187 39.352477 -123.602541 -19.67 23.47 0.52 0.30 0.0334 P188 38.667856 -123.229563 -18.63 27.86 0.40 0.32 -0.1047 P189 38.987449 -123.348428 -19.21 25.92 0.64 0.21 -0.0219 P190 39.241955 -123.204047 -15.87 17.46 0.15 0.43 -0.0286 P192 39.319728 -123.105189 -13.64 14.51 0.41 0.58 -0.0229 58 Table 3.3, continued site vE vN E N EN P194 38.185716 -122.816257 -18.69 28.07 0.48 0.40 -0.0337 P195 38.664665 -122.959307 -15.34 25.16 0.45 0.37 0.0024 P196 38.298144 -122.742641 -16.75 25.13 0.12 0.19 -0.0231 P197 38.428561 -122.767385 -16.16 24.77 0.65 0.15 -0.0216 P198 38.259874 -122.607450 -17.24 21.91 0.79 0.91 -0.0195 P199 38.263692 -122.503438 -12.62 20.62 0.33 0.26 -0.0275 P200 38.239831 -122.451701 -13.96 19.20 0.24 0.14 -0.0252 P201 38.559806 -122.658434 -11.97 17.82 0.58 0.48 -0.0317 P202 38.423581 -122.496000 -12.66 17.24 0.48 0.39 -0.0733 P203 38.866113 -122.917004 -9.92 14.84 0.46 0.37 -0.0384 P204 38.666500 -122.710533 -11.59 18.62 0.47 0.38 -0.0184 P205 39.398124 -122.963064 -11.87 12.06 0.33 1.02 -0.0362 P206 38.777817 -122.575793 -11.48 14.78 0.36 0.20 -0.0270 P207 39.260045 -122.718567 -8.36 7.48 0.66 0.58 -0.0350 P209 37.069255 -122.126714 -25.73 32.76 0.47 0.38 0.0271 P210 36.816137 -121.731841 -27.36 32.83 0.11 0.34 -0.0244 P211 36.879176 -121.698037 -25.49 32.53 0.44 0.35 -0.0095 P212 36.962010 -121.862733 -25.79 32.56 0.14 0.05 -0.0536 P213 37.201712 -121.990836 -20.96 27.95 0.33 0.26 -0.0314 P214 37.001020 -121.796568 -23.11 35.19 1.51 1.51 0.0108 P215 37.048785 -121.762949 -21.29 27.67 0.44 0.35 -0.0287 P216 37.002426 -121.726207 -22.03 28.01 0.43 0.35 -0.0630 P217 37.104497 -121.650632 -19.45 22.67 0.32 0.25 -0.0534 P218 37.203514 -121.713959 -20.17 22.92 0.34 0.27 -0.0134 P219 37.342490 -122.284820 -19.26 30.18 0.69 0.61 -0.0281 P220 37.329887 -122.214284 -22.65 30.13 0.52 0.43 -0.0299 P221 37.336952 -122.099049 -20.11 26.82 0.52 0.44 -0.0312 P222 37.539239 -122.083263 -19.51 24.03 0.34 0.27 -0.0156 P223 37.722054 -122.099779 -15.45 17.72 0.50 0.41 -0.0383 P224 37.863896 -122.219059 -14.68 18.53 0.36 0.12 -0.0258 P225 37.713866 -122.058327 -14.59 17.26 0.30 0.23 -0.0205 P226 37.336776 -121.825585 -20.08 23.55 0.39 0.32 -0.0200 P227 37.532975 -121.789597 -11.82 11.71 0.38 0.31 -0.0246 P228 37.601836 -121.686936 -9.80 10.78 0.36 0.29 -0.0180 P230 37.818965 -121.786399 -9.69 9.74 0.30 0.23 0.0254 P231 36.621677 -121.905409 -28.20 35.46 0.40 0.32 -0.0430 59 Table 3.3, continued Site vE vN E N EN P232 36.724020 -121.579041 -27.05 34.27 0.44 0.35 -0.0440 P233 36.800426 -121.420283 -17.31 24.61 0.38 0.29 -0.0308 P234 36.858526 -121.591225 -26.43 31.50 0.37 0.29 -0.0305 P235 36.814272 -121.541540 -26.16 33.98 0.46 0.37 -0.0248 P236 36.903545 -121.554460 -21.13 25.50 0.12 0.42 -0.0238 P237 36.637023 -121.386829 -26.45 33.65 0.43 0.34 -0.0637 P238 36.849080 -121.452768 -17.05 23.65 0.67 1.37 -0.0267 P239 36.962462 -121.547786 -15.00 23.81 0.61 0.52 -0.0289 P240 37.007807 -121.542032 -18.83 23.59 0.07 0.07 -0.0331 P241 37.213011 -121.574259 -11.22 11.76 0.50 0.42 -0.0294 P242 36.953932 -121.463182 -18.89 22.12 0.41 0.44 -0.0183 P243 36.918189 -121.335152 -10.78 9.43 0.45 0.36 -0.0476 P244 37.010824 -121.354517 -12.46 9.71 0.15 0.11 -0.0335 P245 37.713116 -119.706121 -9.52 5.34 0.60 0.51 -0.0364 P248 37.975609 -121.868697 -9.73 8.37 0.49 0.40 -0.0252 P250 36.950036 -121.268437 -11.19 10.26 0.45 0.36 -0.0356 P251 36.811450 -121.347949 -11.78 13.44 0.14 0.06 -0.0416 P252 37.169556 -121.057729 -10.05 8.41 0.12 0.10 -0.0257 P253 37.478450 -121.653009 -11.61 8.90 0.53 0.44 -0.0207 P254 37.489633 -121.469195 -10.02 9.42 0.49 0.42 -0.0252 P255 37.581885 -121.324852 -10.02 8.49 0.12 0.13 -0.0277 P256 37.931965 -121.604839 -10.75 7.94 0.11 0.19 -0.0231 P257 37.755291 -121.464034 -10.61 8.01 0.13 0.04 -0.0307 P258 37.385389 -121.283272 -9.48 8.39 0.49 0.41 -0.0454 P259 37.433020 -121.100579 -9.72 8.10 0.32 0.24 -0.0225 P260 37.711221 -121.083952 -10.47 7.80 0.11 0.71 -0.0163 P261 38.152960 -122.217539 -9.29 12.83 0.72 0.90 -0.0219 P262 38.025148 -122.096143 -12.98 14.09 0.70 0.42 -0.0475 P263 38.577694 -122.429176 -11.43 13.75 0.60 0.52 -0.0496 P264 38.444215 -122.195329 -10.97 10.54 0.16 0.14 -0.0315 P265 38.530186 -121.954194 -9.59 7.59 0.05 0.05 -0.0234 P266 38.183968 -121.843526 -10.42 8.11 0.13 0.05 -0.0287 P267 38.380336 -121.823234 -8.46 7.88 0.77 0.41 -0.0183 P268 38.473526 -121.646411 -11.43 7.34 0.76 0.14 -0.0145 P269 38.999525 -122.354553 -6.90 6.04 0.62 0.52 0.0333 P270 39.243770 -122.055213 -8.84 6.31 0.57 0.19 -0.0235 60 Table 3.3, continued site vE vN E N EN P271 38.657351 -121.714549 -11.08 6.60 0.37 0.10 -0.0206 P272 39.145482 -121.943062 -9.26 5.70 0.05 0.05 -0.0281 P273 38.115813 -121.388077 -10.46 10.16 0.33 0.25 -0.0182 P274 38.283124 -121.460512 -13.11 4.92 0.97 1.72 -0.0079 P275 38.321529 -121.214586 -10.42 6.66 0.37 0.29 -0.0204 P276 38.645308 -121.095165 -10.35 7.10 0.05 0.05 -0.0272 P277 37.192372 -122.366883 -26.10 35.54 0.47 0.38 -0.0420 P303 37.054383 -120.705295 -9.76 8.85 0.12 0.11 -0.0269 P305 37.352214 -120.196758 -9.22 8.04 0.44 0.36 -0.0337 P306 37.795167 -120.644458 -10.24 7.45 0.44 0.36 -0.0332 P309 38.089991 -120.951237 -9.85 7.74 0.59 0.10 -0.0286 P312 39.529179 -123.698374 -19.97 23.08 0.37 0.29 -0.1053 P313 39.554289 -123.564544 -19.93 20.65 0.37 0.29 -0.0383 P314 39.685664 -123.581847 -20.76 19.08 0.42 0.19 -0.0216 P315 39.863583 -123.716896 -19.87 19.79 0.17 0.24 -0.0421 P317 39.905670 -123.551850 -15.94 13.24 0.37 0.28 -0.0014 P318 39.452370 -123.371853 -17.33 14.59 0.55 0.50 -0.0427 P319 39.707097 -123.295078 -15.52 8.41 0.52 0.48 -0.0344 P320 39.492827 -123.186785 -15.71 12.57 0.81 0.76 -0.0160 P321 39.897221 -123.394947 -14.02 15.87 0.69 0.60 -0.0087 P322 39.962144 -123.182481 -12.29 4.91 0.72 0.61 -0.0613 P323 39.845748 -122.963442 -9.27 9.79 2.27 2.34 -0.0316 P333 39.621381 -122.975753 -9.71 7.58 0.47 0.39 -0.0287 P334 39.493592 -122.735871 -10.43 6.01 0.84 0.78 -0.0314 P335 39.726190 -122.873637 -6.38 6.49 0.65 0.57 -0.0057 P336 39.528080 -122.430483 -7.19 5.36 0.48 0.39 -0.0419 P339 40.034111 -122.668249 -7.70 6.09 0.45 0.36 -0.0325 P340 39.409339 -123.049822 -13.30 8.39 0.60 0.51 -0.1665 P344 39.929123 -122.027971 -8.12 5.78 0.43 0.34 -0.0441 P345 40.271235 -122.270806 -6.92 5.03 0.43 0.34 -0.0479 P534 37.061231 -122.237610 -26.26 34.75 0.47 0.38 -0.0395 PBL1 37.853054 -122.418946 -21.07 22.80 0.32 0.25 0.0212 POTB 38.202625 -121.935386 -10.12 8.87 0.29 0.22 0.0174 PPT1 37.187085 -122.389951 -27.59 35.61 0.57 0.14 0.0126 PTRB 37.996175 -123.018718 -25.27 36.35 0.05 0.08 0.0113 S300 37.666507 -121.558267 -11.01 8.76 0.07 0.05 0.0203 61 Table 3.3, continued site vE vN E N EN SAOB 36.765303 -121.447184 -27.98 34.62 0.17 0.29 0.0046 SBRN 37.686219 -122.410441 -20.15 26.74 0.09 0.09 -0.0624 SLAC 37.416517 -122.204263 -21.69 28.30 0.26 0.20 -0.0096 SODB 37.166404 -121.925518 -21.47 26.17 0.17 0.25 0.0212 SRB1 37.874358 -122.266962 -16.14 25.54 0.83 0.76 -0.0228 SUAA 37.426905 -122.173288 -19.84 26.75 0.59 0.50 0.0296 SUTB 39.205836 -121.820598 -9.67 6.78 0.03 0.11 0.0179 SVIN 38.033179 -122.526318 -18.87 24.16 0.32 0.16 -0.0179 THAL 37.351486 -121.935488 -18.30 23.96 0.32 0.24 0.0191 TIBB 37.890875 -122.447601 -19.16 24.58 0.23 0.14 -0.0088 UCD1 38.536240 -121.751231 -9.78 8.09 0.24 0.23 0.0107 UCSF 37.762968 -122.458150 -19.52 28.63 0.79 0.72 -0.0411 WIN2 37.652652 -122.140628 -19.54 21.79 0.76 0.68 -0.0355 WINT 37.652643 -122.140567 -19.09 23.40 0.14 0.15 0.0195 P193 38.122937 -122.908141 -22.01 32.78 0.45 0.37 -0.0291 62 Table 3.4. Campaign GPS site velocities along the profile, is latitude, is longitude, vE is easting velocity, vN is northing velocity, E is easting error, N is northing error, and EN is the covariance of easting and northing. site vE vN E N EN 1395 38.087308 -122.812645 -22.38 30.48 0.40 0.34 -0.0190 2271 38.532743 -121.949879 -11.84 9.09 1.29 0.98 -0.0128 ADOO 38.236422 -122.527469 -15.90 22.06 0.35 0.31 0.0095 AIRR 38.223184 -122.455673 -14.78 19.58 0.44 0.36 0.0102 CAML 38.416679 -121.994710 -8.63 9.09 0.43 0.37 -0.0160 CORD 38.185992 -122.595287 -17.26 24.14 0.31 0.28 0.0043 CURR 38.475086 -121.822568 -13.28 8.37 1.27 0.98 -0.0211 DEAL 38.257812 -122.337977 -13.96 17.18 0.52 0.47 0.0425 FARB 37.697207 -123.000765 -27.75 38.29 0.33 0.32 0.0088 UCD1 38.536239 -121.751230 -10.36 8.89 0.32 0.31 0.0133 GAME 38.350685 -122.175246 -10.35 11.19 0.35 0.31 0.0325 GORR 38.331178 -122.114607 -10.94 10.56 0.46 0.40 -0.0211 HAGG 38.323878 -122.259202 -11.98 15.44 0.42 0.36 0.0030 HENN 38.283205 -122.361752 -11.70 17.50 2.62 1.50 -0.0789 MADI 38.312607 -122.203192 -12.65 14.91 0.39 0.34 0.0212 NICC 38.092654 -122.736543 -20.45 28.50 0.36 0.31 -0.0026 PRH2 38.079711 -122.868729 -23.67 33.42 0.33 0.29 0.0167 PRH3 37.996405 -123.014906 -25.39 37.02 0.38 0.33 -0.0159 PRNC 38.103525 -122.936598 -31.80 31.54 1.94 1.36 -0.0462 PTRB 37.996174 -123.018718 -26.28 36.84 0.35 0.33 0.0268 VAC3 38.397877 -122.102923 -10.84 9.45 0.46 0.40 0.0507 P261 P267 38.152960 38.380336 -122.217539 -121.823234 -10.05 -9.96 12.77 7.99 0.68 0.68 0.65 0.65 0.0021 0.0073 P198 38.259874 -122.607450 -18.10 23.47 0.68 0.65 0.0020 63 Chapter 4 Earthquake cycle models of the Cascadia subduction zone 4.1 Introduction Subduction zones are convergent plate margins where one tectonic plate overrides another. Typically, oceanic lithosphere (composed of the crust and uppermost mantle) slips beneath the continental lithosphere (or more buoyant oceanic lithosphere), and deeper into the mantle (Figure 4.1). Great earthquakes at subduction zones are the largest earthquakes worldwide, and often cause disastrous tsunamis. In this chapter, I describe models of Cascadia subduction zone earthquake cycles. One goal of such modeling is to understand the dynamics of subduction zones and to estimate absolute stresses in the lithosphere (something existing models still cannot do, for reasons given below). Understanding stresses acting along subduction zone faults is required for us to understand the mechanics of slow earthquakes and non-volcanic tremors (Dragert et al., 2001) and the extent to which these events (which occur about every 400 days) might enhance the probability of subduction zone earthquakes. Another objective is to predict surface motions and compare these with GPS velocity data, as a means to rule out candidate models (e.g. fault zone and mantle rheologies) and to bracket the likely rupture area for future, large subduction zone earthquakes. 4.2 Subduction zones Depending on the age of the subducting plate, subduction zones may be divided in to extensional and compressional subclasses. In general, old thick plates subduct with a large dip angle, causing the continental plate to advance seaward and resulting in back-arc extension. Subduction zones with young, thin, and buoyant oceanic plates subduct with a small dip angle, and the continental plate boundary region is compressed. Subduction zones are categorized into seven subclasses; the end-members are type one (strongly extensional with weak coupling between plates) and type seven (strongly compressional; Stern, 2002). 64 The rheology of the subduction interface is controlled mainly by temperature (Figure 4.1). The shallowest part of the interface, where the accretionary prism (poorly consolidated sediments) overrides the downgoing plate, is aseismic though large ruptures may propagate into this region. The clay minerals which allow stable sliding dehydrate at greater depths and change to frictionally unstable (potentially seismogenic) materials at about 150° C (Hyndman and Wang 1995), so the part of the interface which is located where the temperature is between 150° and 350° is thought to be locked between earthquakes. This part of the subduction zone interface is characterized by velocity-weakening friction. Earthquakes nucleate in this interval though once started, a rupture may propagate both updip and downdip. In general, the zone where the temperature is between 350° and 450° C (called the ―transition zone‖) is conditionally stable; earthquakes may propagate well into this area if the dynamic stress causes enough of a velocity jump. Figure 4.1. Schematic illustration of a subduction zone, showing terms used in the text. The red line indicates the seismogenic zone and the orange dashed line shows the transition zone. (modified from Fluck et al. 1997) Subducted oceanic crust is dominantly basalt, and its most ductile mineral phase is plagioclase feldspar. Feldspar becomes ductile at about 450° C, so the part of the subduction interface with a higher temperature creeps viscously. If an earthquake rupture propagates into 65 this area it will stop very fast because the background differential stress is well below the friction threshold and viscous resistance to sliding is proportional to slip velocity (strain rate; e.g., Scholz 1989). Depending on the thermal gradient of the subduction zone, the depths with the mentioned temperatures (and frictional stability thresholds) above will vary. Episodic tremors and slip events (Dragert et al. 2001) occurs in the conditionally stable depth interval. 4.2.1 The Cascadia subduction zone The Cascadia Subduction Zone extends from the Mendocino Triple Junction offshore northern California to the Explorer Triple Junction south of the Queen Charlotte Islands (Figure 4.2). The Juan De Fuca plate is subducting beneath the North America plate with an azimuth of 69 degrees and a speed of 42 mm per year relative to North America (Flück et. al., 1997). The Juan De Fuca plate is young and thin (~10 km) with a subduction dip angle smaller than 10 degrees (Wang et al. 1995). This indicates that the Cascadia Subduction Zone is a compressional subduction zone. The plate interface has a complicated geometry (Figure 4.2). In Washington beneath the Olympic peninsula its dip angle is less than in southern British Columbia beneath Vancouver Island. Under the Vancouver Island the locked zone extends to 18 km depth and the transition zone extends from 18 to 22 km (Dragert et al. 2001). In this zone earthquakes happen more frequently but take weeks to slip. These events are Figure 4.2. : Geometry of the Cascadia Subduction accompanied by non-volcanic tremor and Zone (modified from Fluck et al. 1997). Contour intervals show depth to the modeled subduction zone are now called ―episodic tremor and slip‖, fault. Red indicates where the Juan De Fuca and North America plates are locked and orange shows the 66 transition zone. or ETS events. Cascadia episodic tremor and slip happen about every 400 days and last for about 2 weeks in any location. Relatively steady aseismic slip occurs at depths greater than 22 km. episodic tremor and slip propagates along strike (to the north and south) at about 6 km/day (Dragert et al. 2001), with more rapid updip propagation at particular locations along strike (Rubinstein et al, 2010). A 2D finite-element thermal model for Vancouver Island and northern Oregon and Washington (Wang et al. 1995) indicates little frictional heating on the Cascadia Subduction Zone Cascadia Subduction Zone. Wang et al. conclude that the Cascadia Subduction Zone fault must be very weak (that is, it must have a low friction coefficient). The result also shows that temperature of the Juan De Fuca plate at the deformation front is more than 200°C and it reaches to 350°C and 450°C respectively at about 14 km and 18 km depth. Therefore this study suggests that the seismic zone in southern British Columbia is from the trench to a depth of about 14 km and that the transition zone extends to about 18 km. This locking depth estimate is shallower than that of Dragert et al. (2001). 4.2.2 Deformation models of subduction zone In his dislocation model of a subduction zone, J. C. Savage (1983) expressed strain accumulation and release at a subduction zone as a result of two kinds of slip. The first kind is periodic slip on the main thrust zone, and the second is steady aseismic slip on the remainder of the plate interface. To model the aseismic slip, he defines a model with slip at the plate convergence rate along whole length of interface, and adds another model with a reverse slip (―backslip‖) at the plate convergence rate along the main thrust zone, that is, the locked part. The sum of these is deformation due to steady slip at depth and a locked (non-slipping) seismogenic interval, that is, interseismic deformation. The sum of the interseismic deformation and the coseismic deformation is deformation over a complete earthquake cycle. This ―backslip‖ model produces no deformation (or uplift) of the continental plate after one complete earthquake cycle. Over long times (many earthquake cycles), interseismic and coseismic deformation cancel and there is no net uplift (or subsidence). This is because the model is purely elastic and the geometry of the subduction interface does not warp. 67 Savage [1983] suggests that the subduction earthquake cycle can be described adequately by an elastic backslip model because the records of interseismic deformation show approximately constant interseismic velocities (at least, long after the last earthquake). One drawback of the model is that it is physically unrealistic, that is, the fault has to be extended to a depth of about 100 km, which is very unlikely, and the Earth is assumed to be a uniform elastic halfspace. This is also a purely kinematic model, that is, slip is specified by the modeler and does not result from forces acting on a plate interface with a realistic creeping rheology. This model also cannot predict absolute stress levels in the lithosphere - just perturbations. Flück et al. [1997] developed a 3-D dislocation model of the Cascadia Subduction Zone based on Okada‗s solution (Okada, 1985). The model is based on the surface deformation due to shear faulting in an elastic half space. It contains curved fault geometry and nonuniform locking zone. The model is a uniform elastic half space with a stress-free, flat surface. This model calculates the surface deformation for a rupture with variable width along the margin of Cascadia Subduction Zone; and it has a linear down dip transition zone between locked and free slip zones. The locked and transition zone are based on the result of previous 2D model of Hyndman and Wang (1995). The results are a narrower locked zone off central Oregon and slightly wider off southern Vancouver Island compared to previous models (Flück, et. al. 1997). Wang (2003) developed a 3-D dislocation model for Cascadia subduction zone to model interseismic deformation. Like Flück‘s model, this model is a uniform elastic half space with a flat surface; however it has an exponential decay in slip in the downdip transition zone to prevent over predicting of coastal and underpredicting of inland GPS velocities (relative to North America). In this zone, slip from the locked to the free slip zone (from 0 mm/yr to the relative plate rate) increases exponentially with distance. Also in this model the effect of secular motion of the central and southern Cascadia fore arc is subtracted to get the effective convergence amount. The method of calculation is same as Flück‘s model (Wang, et. al. 2003). The problem with these elastic half space models with kinematic dislocations is that they require a wide effective transition zone to approximate the stress relaxation in the mantle wedge. Because this relaxation effect is large after the earthquake and gets smaller later on, modeling of the subduction zone at different times in the earthquake cycle with an elastic model requires a transition zone which widens over time since the last large earthquake. This means that an elastic 68 dislocation does not give a physical explanation for the form of the effective transition zone. Another issue is that these models cannot predict the likely width of a coseismic rupture. Both give maximum values because both are based on deformation data from late in the earthquake cycle. Wang (2007) reviews other reasons for caution in the use of backslip models for subduction zone modeling, in particular, the fact that ―reasonable-looking‖ backslip distributions may represent physically unrealistic patterns of actual interseismic slip (for example, a high rate of interseismic slip in a small area which is completely surrounded by a locked interface). Williams and McCaffrey [1991] developed a 2D, elastic finite element model of the Cascadia Subduction Zone, in which the subduction interface was represented with a constant stress boundary condition and only the continental plate was modeled. This model yields estimates of stresses and surface velocities but, like the models above, it assumes a constant rate of interseismic deformation. Viscoelastic models of the subduction zone earthquake cycle account for time-dependent deformation between earthquakes due to relaxation of stresses in viscoelastic mantle and lower crust (e.g., Thatcher and Rundle, 1979; Thatcher, 1984). Wang et al., (2001) developed a 3-D viscoelastic finite element model for the Cascadia Subduction Zone and compared the results with horizontal geodetic deformation observations and uplift data. This model was developed to study the variations of interseismic deformation, which cannot be seen with the elastic models described above. Viscoelastic properties of the mantle asthenosphere were incorporated in the model, unlike the previous dislocation models. Also, a thin viscoelastic layer was used along the down dip transition zone of the Cascadia Subduction Zone interface in order to approximate aseismic creep. Variable convergence direction and speed along strike is used to contain change in convergence direction from normal to trench in the north British Columbia to oblique to the trench in the southern Cascadia. The results show that the surface deformation rate decreases through the interseismic period, and the direction of local maximum contraction is much less oblique than plate convergence. An update to the Wang et al. (2001) model was recently developed by Yan Hu at the University of Victoria (Hu, 2011). This is a 3D, viscoelastic earthquake cycle model but it also makes use of the backslip technique to model interseismic deformation along the Cascadia Subduction Zone. Hence, although it gives time-dependent interseismic stress and velocity 69 estimates, stresses along the Cascadia Subduction Zone are likely inaccurate because the physics of stress-driven afterslip is not being modeled. Liu and Rice (2005) developed a 3D elastic Cascadia Subduction Zone earthquake cycle model which has uniform elastic structure and simplified subduction geometry. However, this model incorporates a sophisticated version of rate-and-state friction, and spontaneously generates both great Cascadia Subduction Zone earthquakes and aseismic slip events. Like the Wang et al. (2001) model, this model yields time-varying interseismic deformation, but it would require more complex elastic structure and viscoelasticity to be more realistic. The earthquake cycle model I propose below will have some of the elements of the Liu and Rice model and the Wang viscoelastic model, as well as a more realistic, heterogeneous elastic and rheologic structure. All of the above models have one thing in common: they do not model stress-driven afterslip along the subduction interface. Rather, the afterslip is assigned kinematically. This is done to preserve the geometry of the subduction zone interface and to prevent gradual accumulation of stress and uplift which occurs in models of subduction zones which are driven by far-field relative motions of plates. 4.3 Data We need two kinds of data for modeling the Cascadia Subduction Zone. The first type is data required to build the model and the second is required to calibrate the model. Seismic wave velocity profiles are needed in order to compute elastic moduli for the different rock units in the finite element model, and to define where important contrasts (like the subduction interface) are located. Gravity data are needed to constrain densities of rock units. Even heat flow data can be used indirectly through thermal models. Thermal models are needed to estimate temperatures and from this rheology of rocks (e.g. parameters governing stable frictional slip and viscoelastic relaxation). Geodetic data are the most important for constraining the subduction zone model. A wide range of these data can be used, for example, leveling data, data from classical geodetic networks (triangulation and trilateration), data from strain arrays, tide gauge records and velocity fields 70 from Global Positioning system (GPS) networks. Of these data, GPS networks provide the most important and precise measurements of ―instantaneous‖ deformation. Longer-term geological studies of uplift rates are also used to constrain the average uplift over many earthquake cycles. 4.3.1 Seismic wave velocity and density data The seismic wave velocity structure for the southern Canadian cordillera has been obtained from refraction and wide angle reflection profiles (Clowes et al., 1995). Data from LITHOPROBE and Vancouver Island seismic experiments (Drew and Clowes 1990) in southern British Columbia were used to achieve a two-dimensional P wave velocity model. The closest part of refraction and wide angle reflection lines to the model profile are chosen and projected onto the profile in order to get the P-wave velocity structure. Figure 4.4 shows the position of the finite element model profile across the Canadian Cordillera; the purple straight line indicates the profile. Other lines on the figure are refraction and wide angle reflection lines; the colored parts of these lines were chosen to get the P wave velocity structure. The P wave velocity in the accretionary prism is 4 km/s, and for other shallow layers it varies from 5.8 km/s to 6.1 km/s. With increasing temperature and pressure (toward the lower crust) P wave velocities increase to the range of 6.3 km/s to 6.8 km/s. For the mantle the typical value of 8 km/s is assumed. Figure 4.3. Figure illustrates the densities in kg/m3 and magnetic susceptibilities in SI units along the profile. (Clowes et al., 1996). 71 Densities of crustal layers from the oceanic basin to the Cascadia volcanic arc are also required to compute elastic moduli. A density model was obtained with the use of potential field (gravity) and seismic reflection data (Clowes et al., 1997). The density of blocks depends on the type of rocks. In the accretion zone the density is about 2000 upper crust it varies from 2290 to 2930 and in other parts on the . In general, for the same kind of rock, density increases with depth. Density estimates are available just for 450 km landward distance from trench Figure 4.3). I extend the end blocks of density to the end of the model. This part of the model is far from the subduction interface, so we assume that the approximation is adequate. Figure 4.4. Red line indicates the project profile across the Canadian Southwest Cordillera. The other colored lines show the parts of the seismic refraction lines which are used for computing the elastic parameters in this project (modified from Clowes et al., 1997). 72 4.3.2 Calibration data: GPS velocities and long-term uplift rate The Geological Survey of Canada (GSC) established the Western Canada Deformation Array (WCDA), a permanent GPS network, to monitor crustal deformation in southwestern British Columbia and study seismic hazard in the region. The Geological Survey of Canada, in cooperation with the Geodetic Survey of Canada, installed the first GPS site (DRAO at Penticton, B.C.) in 1991. Since then 11 additional GPS sites have been added to address several crustal deformation in the area (http://gsc.nrcan.gc.ca/geodyn/wcda/index_e.php). Figure 4.5 shows the locations of the Western Canada Deformation Array GPS sites in Vancouver Island and south-west British Columbia. Figure 4.5. Locations of GPS sites on Vancouver Island and southwest British Columbia are shown in blue squares and red circles. Red circles indicate campaign-mode GPS sites and blue squares are continuous GPS sites. The black line shows profile line of the finite-element model. 73 There are also plenty of campaign GPS sites in and around Vancouver Island (Mazzotti et al. 2003. 4.4 Finite element modeling of the Cascadia subduction zone In this project I developed 2D elastic and viscoelastic models of the Cascadia Subduction Zone. The modeled profile is perpendicular to deformation in southern British Columbia. It starts 200 km from the deformation front offshore (that is, where the Cascadia Subduction Zone fault intersects the seafloor), passes across Vancouver Island and continues 1000 km along northeast direction (black line on Figure 4.5). The model profile crosses the Insular, Coastal, Intermountain, Ominica and Foreland belts. To avoid boundary condition effects we have to continue the model to 800 km depth or use infinite element at the deep end of 200 km-deep model. I used GeoFEST (Parker et al. 2003), a finite element simulation code for the modeling the Cascadia Subduction Zone. It is a 2 or 3 dimensional finite element package for modeling of stress and strain in geophysics and other continuum domain applications. I also used TECTON, a popular finite element package for geophysics (Williams and Richardson 1991), for comparing and checking the results of some calculations. The principles of finite-element modeling for continuum mechanics problems (as implemented in both of these codes) are described in Appendix A. 4.4.1 Elastic structure The first ingredients required for finite element modeling of an elastic body are its elastic moduli, so we need to have the elastic moduli of lithosphere across the subduction zone. Because direct access to crustal layers to measure elastic moduli is impossible, we have to get them from seismic wave velocity and density of the layers (or rock units for heterogeneous models). 74 To get the elastic modulus of rock units, I used density and P wave velocity models and partitioned the profile so that each partition has an approximately constant P wave velocity and density. Then these values were used to compute elastic moduli for cells by the (2.4) relation. Figure 4.6. : Bulk modulus values for the modeled profile, based on Lithoprobe seismic velocity and density data (see text). Partitioning the model by the change of elastic moduli plays a significant part in the process, because it may affect the deformation dramatically. Figure 4.6 shows the difference of vertical surface deformation due to a dislocation, for uniform and complex elastic structure. In this test model, the modeled slip is 24 m – representing the total of Cascadia Subduction Zone convergence in an average cycle of 500 years - and the dislocation geometry is as shown on Figure 4.6 (see next section for details). The red dashed line is the result for the uniform elastic structure. In this case density and P-wave velocity are 2700 and 6.4 km/s, respectively. The blue solid line is the result for the model with complex elastic structure (shown in Figure 4.6). Modeled subsidence of Vancouver Island for the heterogeneous elastic case is almost 2 m more than for the uniform elastic model. Masterlark (2003) has also shown that elastic heterogeneity affects coseismic (and postseismic) surface deformation for a general test model of a subduction zone. 75 Figure 4.7. Figure shows modeled coseismic surface subsidence for two elastic finite element models. Red dashed line belongs to the uniform elastic model and blue solid line indicates the result of the heterogeneous model. 4.4.2 Geometry of the Cascadia subduction interface The Cascadia Subduction Zone has a complex geometry. Along strike, the subduction interface is shallowest in Washington State, dipping more steeply in southern British Columbia to the north and in Oregon to the south. At the edge of the continental plate, there is a lowseismic velocity and low-density sediment layer (accretionary prism) above the subducting Juan de Fuca plate. This sediment is scraped from the subducting oceanic plate and added to the continental plate. It extends about 100 km inland from the trench, and the dip angle at this part of the subduction zone is less than 10° (Figure 4.6). Further inland, there is debate about the location of the subduction interface. Two highly reflective regions, the C and E layers (Green et al., 1986; Hyndman et al., 1990) are interpreted either as shear zones or horizons with locally high pore fluid pressures. These layers are separated by a 5- to 10-km thick interval of apparently stronger rock (either underplated ocean crust or lower continental crust; Hyndman et al., 1990). The E layer, which is 25 to 45 km below Vancouver Island, is a zone with a prominent low Swave velocity, a high inferred Poisson‘s ratio, and an estimated thickness of about 5-10 km (Audet et al., 2009). Audet et al. (2009) conclude that it is a zone of dehydrating ocean crust with high pore fluid pressures, and that it represents the subduction interface. Earlier studies (Green et al., 1986) also infer that the E layer is the subduction interface. A deeper interface (the F reflector) has also been interpreted as the active subduction interface (e.g., Drew and Clowes, 76 1990) and for this study; I assume that this is the case. At 30 km depth, the apparent subduction interface dip angle increases to more than 17° (Clowes et al. 1995). In contrast to the geometry of the subduction interface in the crust, in the upper mantle the down going oceanic slab dips more steeply beneath northern Washington than beneath British Columbia. The dip angle of the slab in the British Columbia upper mantle is on average 50°, however it is 60° beneath northern Washington (Bostock et al., 1994). 4.4.3 Meshing the model When developing a finite element mesh, the following guidelines are important: 1) It should be possible to mesh different parts of the model (zones) independently, with a pre-determined element size. 2) Element size should increase gradually with distance from fault. 3) Elements should be as close to equant as possible. High aspect ratios may result in computational instability after even modest element strains. I wrote the software ―Mesh‖ to automatically mesh the model domain, while satisfying the above conditions. This software contains 34 files, which are in total more than 8000 lines of MATLAB code. This software can mesh a 2D model using triangle or quadrilateral elements, and can mesh a 3D model using hexahedral elements. Mesh has graphical outputs which display the meshed model, zones and their properties, boundary conditions, split and slippery nodes, and output nodes. It also shows modeled element stresses and nodal displacements, which are calculated by the finite element models (GeoFEST or TECTON). The software also generates input files for GeofEST and TECTON. It enables me to change the geometry and parameters of huge models easily. To make such a change without this software sometimes is a very tough job and requires several hours or even days. (See Appendix B for a full description of Mesh, and Appendix D for input and output file formats.) 77 Figure 4.8. Part of the unstructured Cascadia profile mesh. Colors indicate zones with different elastic modulus values. Note that the elements are smallest (1 to 2 km in dimension) at the subduction interface. Also note that the surface topography is represented. 4.4.4 2D Cascadia elastic model: Mesh testing and benchmarking The code Mesh requires information about the crust and mantle structure, including the geometry of the subduction interface. The model profile is meshed using both triangle and quadrilateral elements, and the element dimensions increase gradually with distance from the fault. The subduction interface is modeled using the ―split node‖ technique (Melosh and Raefsky, 1981). To assess whether the model mesh is fine enough for the problem I am addressing, I have varied average element size (from a coarse mesh to a fine mesh) and compared modeled displacements and element stresses for the elastic case. The element size is refined until modeled 78 displacements become insensitive to element size. I also set a minimum element size so the modeled strain does not excessively distort the elements; otherwise the elements will collapse or the results will be inaccurate. To eliminate the problem of mesh distortion, the model should be re-meshed during long simulations. I did not have such a condition in these models. These models are run for 10 to 25 earthquake cycle and for each cycle the largest displacement is 20m, so for the whole process the displacement is 200 to 500 m. The finest element of the models is about 2 km, so the maximum strain is 10% to 25% which is acceptable. Specified velocity boundary conditions at the edges of the model domain can affect the solutions. For the elastic models, I assume that the side boundaries are fixed and that the top boundary is ―free‖, that is, it is a zero-stress boundary. To avoid undesirable boundary effects on the solution, different-sized models with simple fault geometry were developed. Their results were compared with each other and with the analytical solution to find the right size of the model domain. The test models incorporate uniform elastic properties and a straight, dipping fault (rather than the curved, Cascadia geometry) so the model results can be compared with the analytical solution. The modeled fault is straight, with a dip of 10 degree. It extends from the surface to a depth of 18 km. 24 meter of slip are modeled. I began with a 200-km deep mesh, and found that the finite element model failed to yield the surface deformation predicted by the analytical solution for a dipping dislocation in a uniform elastic halfspace (Freund and Barnett, 1976). The same problem occurs for both TECTON and GeoFEST. To better emulate an infinite half space I had to extend the model domain to a depth of 800 km or use infinite elements at the bottom of the model. Figure 4.9 illustrates the results of modeling for surface vertical and horizontal displacement. The TECTON model is 200 km deep, but with an infinite element at the bottom and sides. The GeoFEST model domain is 800 km deep (infinite elements are not implemented in GeoFEST). Both of the finite-element models match the analytical solution, and the large model domain dimensions required by the GeoFEST model were incorporated in my Cascadia models. 79 Figure 4.9. A comparison of numerical and analytical solutions for surface displacements, for an elastic model with a straight dislocation. This comparison shows that both GeoFEST and Tecton can replicate the analytical solution, provided that the model side and bottom boundaries are sufficiently far from the subduction zone. 4.4.5 Incorporating viscoelastic relaxation and stress-driven slip Accurate geodetic deformation data are available only for a short period of time in comparison with the subduction zone earthquake cycle. In order to use these data to calibrate a subduction zone model, we actually need to model the entire earthquake cycle. To prepare such a model, viscoelastic stress relaxation and stress-driven interseismic slip (for example, via ratedependent friction) should be included, and deformation should be driven by remotely applied plate-driving forces. Such an ideal model for subduction zones is not available yet, though models addressing some of these three requirements have been developed. The finite element method has been used to develop viscoelastic earthquake cycle models of subduction zones (e.g., Wang et al., 2003; Hashimoto et al., 2003; Taylor et al., 1996), but all of these models make use 80 of backslip to represent interseismic creep of the Cascadia fault zone. The problem with backslip is that it is a kinematic approach; hence, when it is used, stresses in the plate and along the subduction interface are not modeled correctly. In this section I describe how viscoelasticity and stress-driven fault slip are incorporated into the FE model. In the subsequent section on tectonic loading and displacement of the subduction interface (4.4.6), I describe in detail my efforts to model the earthquake cycle without resorting to backslip, to get closer to an ideal, simultaneous treatment of the subduction interface and tectonic loading. Though several approaches were tried, none were ultimately successful for viscoelastic models. 4.4.5.1 Viscoelastic Relaxation Viscoelastic materials exhibit both viscous and elastic characteristics when undergoing deformation. In viscous materials, shear stress is proportional to strain rate. Elastic materials strain instantaneously when stressed, and instantly return to their original state once the stress is removed. Viscoelastic materials have elements of both these properties and exhibit time dependent strain. When stress is applied to such material it starts to strain, and if strain is held constant, stresses relax exponentially with time. The upper mantle of the Earth is viscoelastic. When an earthquake happens and stress changes in the upper mantle, stresses in the mantle gradually decrease with time. This causes deformation in the overlying crust. So the surface deformation is in part due to mantle viscoelastic relaxation. The other part of surface deformation is a result of locking of the seismic and transition zone. If we ignore the viscoelastic effect, the model will incorrectly determine the Cascadia Subduction Zone interface stresses, and potentially the limits of the locked and transition zones. In subduction zones, viscoelastic deformation is clearly responsible for a significant part of the observed interseismic surface deformation (e.g., Thatcher and Rundle, 1984). The effective viscosity of the Cascadia lower crust and upper mantle has been estimated from modeling studies and inferred based on estimated temperatures, volatile content, and composition. In most earthquake cycle modeling studies, a Maxwell viscoelastic rheology is assumed. Based on post-glacial rebound, the mantle viscosity in the Cascadia Subduction Zone 81 region is 5 to 50 x 1018 Pa s (James et al., 2000). This value falls within the range given by experimental studies on ―wet‖ olivine with temperature, pressure, differential stress and volatile content characteristic of the upper mantle at a subduction zone (Dixon et al., 2004; Currie et al., 2004; Bürgmann and Dresen, 2008). 4.4.5.2 Stress-Dependent Slip ―Split nodes‖ may be used to represent the fault interface in finite-element models, and may be used to model both seismic and aseismic slip [Melosh and Raefsky, 1981]. Split nodes require that the slip per timestep be provided, either by the modeler or from the previous time step. (―Slippery‖ nodes, the other commonly-used fault representation in finite element models, merely slip by the amount that is consistent with zero shear stress on the fault.) The amplitude of split-node displacement for aseismic slip is equal to the coseismic slip which accommodates the total plate convergence during an earthquake cycle. To avoid a stress singularity and have a better physical representation of the rupture tips, coseismic slip is tapered at the both ends of the seismic zone. Along the transition zone, coseismic slip decreases from full slip at the seismic end to zero at the aseismic end. Far downdip, aseismic slip occurs at the full plate convergence rate every time step. In between, the slip rate depends on stresses and varies with time in the earthquake cycle, but the total integrated slip sums to the full convergence rate times the interseismic interval. The transition zone is the most controversial part of the Cascadia Subduction Zone. This part is likely velocity-strengthening; an earthquake cannot nucleate in it, but a rupture may propagate into it (e.g., Tse and Rice, 1986). It is often assumed that slip tapers linearly or exponentially from full slip to zero across the zone (Hyndman and Wang 1995, Wang et al. 2001, Flück et al. 1997, Wang et al. 2003, Wang, 2007), but to gain a physical understanding of the Cascadia Subduction Zone, we have to apply stress-related slip rate to the transition zone. For this purpose we have to obtain stresses at the split node locations, compute the slip rate based on these stresses, and apply the new value of integrated slip to the split node at the next time step. Stress at a node is the average of stresses of all the elements around the node. To get the shear stress at the split node we must rotate the stress tensor such that one of the coordinate axes is parallel to 82 the fault plane, so the shear stress is the off-diagonal component of the symmetric, 2x2 stress tensor. To get the shear stress for a node we need to take an average for all elements connected to the node. To compute the stress-dependent slip rate, we assume either of two conceptual models of the fault zone. For the first, we can assume that there is a thin layer of low-viscosity material along the interface (similar to the idealization of Wang et al. [2003]). For the viscous shear zone we have: (4.1) Where is stress, is viscosity of the shear zone and is strain rate. If we use shear stress we will have: xy xy xy w s t s w xy t (4.2) t is time interval, s is slip amount in t and w is the width of the viscoelastic layer; for a time unit we will have s xy w (4.3) For the second approach, we assume that aseismic slip follows a velocity-strengthening friction law. We use the ‗slowness law‘, which agrees with experimental data better than other rate/state-variable constitutive laws (Scholz, 1998), to determine slip rate. The friction at steady state is: V 0 (a b) ln w V0 (4.4) Where ’ is shear stress over width of shear zone which we adjust it on our model, is effective normal stress, V is slip velocity and V0 is a reference velocity, a and b are material 83 property and is the steady-state friction at V=V0. By substituting V=s/t and V0 = s 0 /t where s and s0 respectively are slip and reference slip, we can get slip rate as: 0 V V0 exp ( a b ) ( a b) 0 s s0 exp ( a b) ( a b ) (4.5) The FE modeling can be adapted to model stress-dependent, aseismic fault slip using the approach outlined above. A similar approach has been used for modeling postseismic deformation following the 1989 Loma Prieta, the 1999 Izmit, and the 2004 Parkfield earthquakes (Linker and Rice, 1997; Hearn and Bürgmann, 2002; Johnson et al., 2006), as well as earthquake-cycle deformation along straight, strike-slip faults (Hearn, pers. communication). 4.4.6 The problem of tectonic loading and displacement of the subduction interface There are several forces acting on a subducting oceanic plate. Together, these forces cause the oceanic plate to move toward the trench, bend, and slip beneath the continental plate into the mantle at an angle that in the case of the Cascadia Subduction Zone, is fairly steady. The first force is acting under the oceanic plate because of viscoelastic coupling between the plate and the mantle asthenosphere. The second force is like the first one, but between the continental plate and the asthenosphere. The third is at the mid ocean ridge where the plate is pushed away by the gravitational potential energy gradient (―ridge push‖), and the fourth is a resisting force of transform faults at the mid ocean ridge. At the subduction zone there are at least three other forces: the frictional force between continental and oceanic plate because of relative motion between them, the negative buoyancy force of the down going oceanic plate (or ―plate pull), and viscous resistance of the mesosphere when the plate plunges into it (Forsyth et. al. 1975). The balance of these and other forces determine the subduction interface geometry, as well as whether the subduction zone is compressional or extensional (rolling back) and whether permanent deformation of the continent (averaged over many earthquake cycles) should occur. Along the subduction interface, the downgoing plate bends and the subduction interface is curved. We do not explicitly model all of the forces that cause this flexure, but the result is that 84 subduction is accomplished along this curved geometry without offset of the subduction interface or distortion of the overriding plate. The permanent uplift rate of western North America along the Cascadia Subduction Zone for the most part is less than 0.25 mm/yr (Kelsey et al. 1994), so compared with interseismic deformation it is negligible. It means that the subduction interface geometry is not shifting. If there were not a locked part of Cascadia interface, there would not be significant deformation of the continent. Since there is locked zone, it causes strain accumulation during the earthquake cycle. Ideally, in an earthquake, this strain is released and there is no net deformation. In order to model a subduction zone dynamically, we have to apply a far-field tectonic load as a horizontal velocity boundary condition. But continues horizontal motion of the oceanic plate in the trench will compress the continental plate of the model and push the subduction interface toward the continent. As the geometry of the interface distorts, the continental plate in the model over the subduction zone will compress. Superposition of this deformation over repeated earthquake cycles produces significant uplift of the back arc over geological time. However, there is no evidence of such uplift (e.g. Kelsey et al., 1994). This deformation will occur in my model if I do not require the Cascadia Subduction Zone interface to be stable. I will call this steady (non earthquake-cycle) deformation ―net deformation‖ in subsequent sections. To solve this problem in my finite element model, I have to improve the model so that the geometry of the subduction interface does not change. We don‘t have enough information about all of the forces acting on a subducting plate to model them explicitly. Even if they were known, it would require a very huge, complicated model. Therefore, we have to apply velocity boundary conditions at the sides of the model, and use constraint equations in the model to obtain a stationary subduction zone interface. This requires editing the finite-element code by adding correction terms to the nodal displacement equations. The correction terms are applied to either the nodal force vector or directly to the nodal displacements after each time step (see below). Note that I require the plate interface to be stationary over an entire earthquake cycle. Within the seismic cycle, or coseismically, the subduction interface may displace. Besides the net deformation, the interface of the subduction zone deforms elastically by an amount equal to the total plate convergence along the interface after a whole earthquake cycle. I will call this recoverable deformation ―elastic displacement‖ in subsequent sections. Figure 4.10 shows 85 schematically the trajectory of each fault node along the earthquake cycle. P and P′ are respectively position of a node at the beginning and at the end of the earthquake cycle on the down going plate. The curve and line between P and P′ show the trajectory of the fault node during the earthquake cycle. Figure 4.10. Schematic trajectory of split node displacement over a seismic cycle. P represents the mean position of the split node. The total interseismic displacement includes a recoverable elastic component and a net (permanent) component. For the complete seismic cycle, the elastic component is zero normal to the subduction interface. My goal in applying correction terms to the FE equations is to eliminate the permanent component. After one complete cycle, the displacement of the fault node includes a component of motion along the subduction interface and a net motion normal to the interface which should be removed. Because the plate convergence rate is constant, we assume that the net deformation is approximately linear with respect to time. So we have to apply a constraint equation to fault nodes every time step and force the node to move back by the net deformation amount for the time step. This will not ensure that the subduction interface is stationary at all times, merely that it is stationary over an entire cycle. Figure 4.11 schematically illustrates the total displacement of a split node (U), its net (permanent) displacement (UM) and its elastic displacement (UE) for each time step. The constraint equations are 86 Figure 4.11. Vectors indicating elastic displacement, model-net (permanent) displacement and total displacement vector for a time step. These vectors are exaggerated in scale. Elastic displacement of the split node occurs solely along a tangent to the subduction interface. The goal of corrections to the FE equations is to remove the model-net displacement component (red vector), which continuously accrues over multiple earthquake cycles (and is inconsistent with the lack of significant long-term uplift of Vancouver Island). U x XC g1 (u ) U x X C U z Z C g 2 (u ) U Z Z C (4.6) XC and ZC are respectively x and z direction component of UM , UX and UZ are displacements of the node respectively in x and z direction. g1(u) and g2(u) are constraint equations. Now we can use the Lagrange multiplier method to restrict the nodes to move only along the Cascadia Subduction Zone interface tangent direction (see Appendix C). In addition, a small fault-normal displacement equal to the coseismic Cascadia Subduction Zone interface-normal displacement times (time step / interevent time) must be added to bring the node back to the Cascadia Subduction Zone interface after a complete earthquake cycle. 87 4.4.7 Corrections to the FE model equations to preserve the geometry of the subduction interface I developed several numerical methods to apply a correction term to the nodal displacement equations, which would preserve the long-term geometry of the subduction interface. All of the refinements described here were made to the Geophysical Finite Element Simulation Tool (GeoFEST 4.7), as a separate module. I tried seven different methods in seven versions of the subduction module. The modules need extra input data like the node numbers for the split nodes along the subduction interface, nodes or elements defining the adjoining ocean plate (slab), nodes and elements of the adjoining continental plate, and boundary conditions for the ideal slip condition on subduction zone, it means at which direction and how much the nodes move when the slab sliding regularly and continually under the continental plate without locking zone. I have two approaches for applying the correction. In the first approach, the correction is applied in terms of displacement. In the second, it is applied in terms of force. For each correction category, several methods are used. Here, just the most important methods are described. Appendix A outlines how the system of equations we use in finite-element modeling is derived and what it means. As noted there, we can write the global system of equations for the finite element method as follows: k11 k 21 k 31 k m1 k12 k 22 k 32 k m2 k13 k 23 k 33 k m3 k1n U 1 F1 k 2 n U 2 F2 k 3n U 3 F3 k mn U n Fm (4.7) The matrix k is the stiffness matrix, U is the displacement vector, and F is the force vector. This system of linear equations is an algebraic representation of Hooke‘s law of elasticity (a boundary-value differential equation) for the entire modeled domain. Each row n represents a degree of freedom for the nodal displacements, that is, for my model, horizontal and vertical displacements for each node whose velocity is not specified. k incorporates information about the 88 nodal coordinates, their spatial derivatives, and the elastic properties of the material, which may vary within model elements (Appendix A). F is zero for most elements, except at nodes falling along model boundaries with imposed velocities or stresses, or where faults are modeled (Appendix A). If gravitation is modeled, the gravitational force associated with each model node is added to the force vector. If F is zero for all nodes (for example, all boundaries are tractionfree and are free to displace) then there is no deformation (and an infinite number of solutions, as the model domain may displace or rotate arbitrarily). 4.4.7.1 Corrections to nodal displacements To apply the correction to the displacement vector U, it is required to compute and correct displacement of all model nodes. The correction vector is separated to two parts, first one belongs to the nodes of slab and the second one belongs to the nodes of continental plate. So to compute the correction vector the whole model is divided to two parts (Figure 4.11). For the continental model, the correction term is negative amount of displacement vector of steady sliding model (total displacement), because these nodes should be stationary. To compute the correction vector for the slab, we need to find the elastic deformation for the nodes and subtract the total displacement. To obtain the elastic deformation of the subducting slab, it is modeled separately Figure 4.12. Figure shows the model when it is separated to two parts. 89 To apply this method, the first step is to run the model with steady slip. This gives the total displacement (U). Then the model for the slab with the boundary condition explained above (fig 4.11) is run. This model gives the elastic displacement (Ue) for the slab‘s nodes. Finally the correction vector in terms of displacement is computed in two parts as follow: { (4.8) These two parts are assembled to form the correction vector which is added to the displacement vector after each time step. Although this method looks very clear, the result for it is one of the worst I tried, because it causes serious problems for the calculation of the relaxation force (right hand side) at the next time step. It is like a solving a set of equation by a numerical repeating method, in which after each solution you change the some of the value so the process will not converge thus not lead an answer. Figure 4.13 shows the vertical displacement for the model with displacement correction term. Figure 4.13. Vertical displacement of the model with displacement correction. 90 4.4.7.2 Corrections to nodal forces The force correction method involves adding terms to the force vector F in equation 4.7, in order to prevent the net displacement of the nodes representing the subduction interface. (a) First force correction method The first method I developed for computing the force vector correction term involved modeling steady sliding along the subduction interface, computing the forces at the split nodes, and using these forces to adjust the force vector in the earthquake-cycle model. If the displacement vector is the sum of the model-net-displacement vector and the elastic displacement vector in equation (4.7), we will have: k11 k 21 k 31 k m1 k12 k 22 k 32 k m2 k13 k 23 k 33 k m3 k1n Ue1 Un1 F1 k 2 n Ue2 Un2 F2 k 3n Ue3 Un3 F3 k mn Uen Unn Fm (4.9) Ue is elastic displacement and Un is model-net-displacement, so the correction force FCor will be: k11 k 21 k 31 k m1 k12 k 22 k 32 k m2 k13 k 23 k 33 k m3 k1n Un1 F1cor k 2 n Un2 F2cor k 3n Un3 F3cor k mn Unn Fmcor (4.10) The final equation which gives elastic displacement, Ue, will be: k11 k 21 k 31 k m1 k12 k 22 k 32 k m2 k13 k 23 k 33 k m3 k1n Ue1 F1 F1cor k 2 n Ue2 F2 F2cor k 3n Ue3 F3 F3cor cor k mn Uen Fm Fm (4.11) 91 As the first step, a model with slip at a constant rate over the modeled Cascadia Subduction Zone interface is run. The displacement that comes out from this model is the sum of model-netdisplacement and the elastic displacement. We know that the elastic displacement of the split nodes for this steady sliding model is sliding down dip along the interface at half of the relative plate rate, so by computing the elastic displacement at each point on the interface and subtracting that from total displacement, obtained from initial model, the model-net-displacement may be computed. By replacing the model-net-displacement in equation (4.13), the correction force may be obtained. Figure 4.14 shows vertical uplift for this method. Figure 4.14 Vertical uplift of the model with force correction term on interface nodes. (b) Second force correction method In the second force vector correction method, the correction force is applied to all the nodes of the element which touch the interface, rather than just the split nodes themselves. In a triangular element like element number 11 in figure (4.12) which touches the interface at node 2, the local equation is: 92 k11 k12 k 21 k 22 k 31 k 32 k13 U 1 F1 k 23 U 2 F2 k 33 U 3 F3 (4.12) This equation is simplified shape of general equation and there is just one line for each node. There should be 2 equations for each node (in x and y direction). As with the first method if the displacement is divided to 2 parts – model-net-displacement and elastic displacement- the equation for the correction displacement, which is negative amount of model-net-displacement, is: k11 k12 k 21 k 22 k 31 k 32 k13 0 F1Cor k 23 Un2 F2Cor k 33 0 F3Cor (4.13) Correction is applied just for the interface node, so for node number 2 and 3 it is zero. Finally correction force for the nodes of the element which touches interface (element number 11 in figure 4.12) is: F1Cor k12Un2 F2Cor k 22Un2 F3Cor k 32Un2 (4.14) Figure 4.15. Figure shows elements which are touching each other on the interface and have a common nodes over the interface. 93 Final correction force of each node is the sum of correction force for the node in all the elements which the node belong to. For example the global correction force of the node number 2 in element 11 is the sum of all correction force of the node in elements number 10, 11, 12, 47 and 48 (figure 4.13). Beside the interface nodes at this method there are correction force for the node 1 and 3, on the other word correction force is applied for the all nodes of the elements which touch interface. As with the first method, an initial model with steady sliding along the interface is run. The displacement from the initial model is the sum of model-net-displacement and elastic displacement. Elastic displacement is known, the elastic displacement for this steady sliding model is sliding rate at the interface direction toward down, so by computing the elastic displacement for each nodes on the interface and subtracting that from total displacement, model-net-displacement for interface nodes are obtained. Replacing them in equation 4.12 and following the procedure provides the correction force. Figure 4.16 shows the vertical displacement of the method. Figure 4.16 Vertical displacement of model with force correction term on the all nodes of the elements which touch the interface. 94 (c) Third force correction method At this group the method which works better than others is the one which uses slab model to compute the correction force. This method separates the Cascadia Subduction Zone model into two models: one of the overriding plate and one of the slab. The model-net-displacement for this method is computed using equation 4.11 and equation 4.14 is used to compute the correction force. Figure 4.17 shows the vertical displacement of result of this method. Figure 4.17 vertical displacement for the model with force correction term computed by modeling slab with boundary condition on the interface nodes. 4.5 Recommended correction method and future research The most acceptable approach for applying the correction terms, from the view of this study, is applying a reaction force to the Cascadia Subduction Zone interface nodes. To compute the reaction force, the model-net-displacement is needed. This is obtained as in the first method of section. The most important part is the second step, when the model-net-displacement as boundary condition and split node must be applied to the interface nodes at the same time. If the finite element software does not have the ability to do this, it must be implemented. Applying the boundary condition and split nodes at the same time gives the right condition to compute the 95 reaction force which should cancel the model-net displacement. The displacement vector which is obtained from the second step (UI) is used to compute the reaction force as: k j1 k j2 k j 3 k jn U 1I I U 2 U 3I F jR U I n (4.15) j indicates the interface nodes; so k is formed from rows of stiffness matrix which is relevant to the interface nodes j and FRj is the reaction force relevant to the interface nodes j. The correction force is the negative value of the reaction force. F Cor F R (4.16) In this new method, tectonic loading is applied and net long term deformation from the model is eliminated. Since earthquake cycle models may be run to the point of cycle invariance, stresses may be estimated, so relaxation of nonlinear viscoelastic layers may be modeled. Stressdependent fault zone creep may also be implemented. This would allow modelers to explore several scientific questions. 4.5.1 Physics of aseismic deformation along the subduction zone The seismogenic zone is often treated in models as a single locked fault patch, though the locked zone may be a combination of locked asperities separated by aseismic fault patches (e.g., Wang 2007). This is the case when the seismic coupling coefficient is less than 1, that is, when the seismic moment release in earthquakes is less than the moment release expected from the relative plate motion. Seismic moment is directly related to rupture dimension. The relation between the log of the seismic moment and the log of rupture area is linear (Kanamori and Anderson 1975). If the seismogenic zone contains aseismic patches, the computed or predicted seismic moment with seismic zone size will be greater than the true moment. Modeling the subduction zone with stress-driven slip will make it possible to compare the models with few 96 asperities along the fault surface, separated by aseismic fault patches, with a model comprising one big asperity, on the locked part of the subduction zone. The method can evaluate whether GPS velocity data can discern whether the interface is fully locked, and if there is surface deformation data on the land above the locked interval. The recent M9 earthquake in Japan occurred where geodetic and seismic studies had given conflicting results on the estimated sizes of future great earthquakes (Mazzotti et al., 2000; Kawasaki et al., 2001, Heki et al., 1997, Petersen and Seno, 1984; Pacheco, 1993). Thus, understanding the physics of stress-driven creep along subduction zones is a first-order question. Quantifying the physics of stress-driven creep along subduction zones is fundamental to understanding the physics of episodic tremor and slip (ETS), though other evidence (e.g. tidal triggering) has illuminated much about the stress state along the subduction interface where these events occurs. 4.5.2 Evolution of stresses throughout the earthquake cycle, and seismic hazard To understand absolute stresses on faults in tectonically active areas, it is not enough to model stress perturbations from earthquakes. Deformation must be modeled until (in theory) a cycleinvariant solution is obtained; that is, stresses have built up to the point where a steady solution (from one cycle to the next) has been obtained. Once this is achieved, stresses and stressing rates may be computed along the Cascadia Subduction Zone interface in order to analyze how these vary with time between large subduction zone earthquakes. The evolution of Coulomb stresses (shear stress plus friction coefficient times normal stress) on active faults in the overriding plate may be used to infer how seismicity rates and earthquake probabilities should evolve over time between great earthquakes(Stein, 1999; Taylor et al., 1997 [for faults in the ocean plate]). This could shed light on where and when potentially dangerous shallow earthquakes will be promoted by great subduction zone earthquakes. 97 4.5.3 Inferred locations of episodic tremor and slip events Rheological structure may affect inferred locking along the Cascadia Subduction Zone, as well as inferred locations of aseismic slip events. Masterlark [2003] analyzed the assumption of a homogeneous, isotropic, Poisson-solid half-space (HIPSHS), for subduction zone models. He found that static deformation is sensitive to the assumptions of a Poisson solid ( = 0.25), isotropy, and homogeneity, but that it is insensitive to topography (i.e., the assumption of a flat, free surface). His findings are consistent with my finding that surface deformation is sensitive to elastic structure (Figure 4.7). Shear zone creep and viscoelastic relaxation in the mantle might also affect geodetically inferred locations of slow slip. Some models of postseismic deformation following great subduction zone earthquakes invoke a Burgers (biviscous) rheology for the mantle with a very low initial effective viscosity (e.g. 5 x 1017 Pa s for the 2004 Sumatra earthquake; Pollitz et al., 2006). If this represents the mantle wedge rheology, then over the ~1 week timescales of slow slip events, assuming a high mantle effective shear modulus based on seismic data could be inappropriate. Earthquake-cycle models could show whether such a mantle wedge rheology would be at least consistent with current GPS observations from Cascadia. 98 Chapter 5 Conclusion Key contributions of my research are summarized in this chapter. I also describe the limitations of my findings and make suggestions for future research opportunities. 5.1 Contributions The main topic of this thesis is characterization of the effect of lower crust and mantle viscosity on the surface deformation around plate boundaries. I make use of the finite element method to model plate-boundary faults and study the effect of substrate viscosity and other parameters on the surface deformation. For my first project (Chapter 2), I developed a suite of earthquake-cycle models to investigate whether asymmetric deformation often noted around strike-slip faults could be due to rheological contrasts at depth. For my second study (Chapter 3), I modeled a continental transform plate boundary where the strain is accommodated by several parallel faults. For the final study (Chapter 4), I modeled subduction zone earthquake cycles. My main scientific findings are summarized below. For all three projects, I developed meshing software in Matlab (―Mesh‖, Appendices B and D). This software automatically meshes a 2D or 3D model domain, which may contain subdomains (partitions) with different material properties, and faults. It is a very flexible code, and makes it easy to change material parameters, fault properties, and the geometries and properties of partitions. It can discretize a complicated model in a few minutes, whereas it would take a long time to do it by hand. By automatically preparing input files for the geophysical finite-element codes like GeoFEST and Tecton, Mesh makes discretizing a model fast and errorfree. 5.1.1 Asymmetric deformation around strike-slip faults is not due to material contrasts For my first project (Chapter 2), I developed a suite of earthquake cycle models of an infinite, vertical strike-slip fault, incorporating a contrast in elastic plate thickness in which plate 99 thickness changes gradually from one side of the fault to the other. I find that that the sensitivity of surface velocity to elastic plate thickness contrast depends on substrate viscosity. Low viscosity (1018 Pa s) causes more asymmetry and high viscosity (1021 Pa s) occasions trivial asymmetry. The width over the zone across the plate transitions from thick to thin has little effect on surface deformation. I also developed test models with a viscosity contrast across the fault (and a uniform elastic plate), and these show the same result. All of the models may cause asymmetric interseismic velocity profiles, but only for models with modest viscosities on one side of the fault, which produce diffuse strain around the plate boundary late in the earthquake cycle I conclude that a contrast in substrate viscosity does not explain the asymmetric deformation seen around some strike-slip faults (e.g., Le Pichon et al., 2005). All of my models that produce this effect require modest viscosities on one side of the fault in contrast with the other side, resulting in a broad and time-varying deformation around the fault at different times in the earthquake cycle. This is not seen around natural faults, where deformation appears to be localized and invariant through most of the interseismic interval. My models of the Point Reyes region (Chapter 3) and other models of different strike-slip faults (e.g., Hearn et al., 2009) indicate that this requires a high viscosity for the lower lithosphere (5x1020 – 1021 Pa s). Hence viscosity contrasts cannot account for the asymmetry. I find that deformation asymmetry is more likely caused by an offset between the surface trace of a fault and its position at depth, where it is creeping aseismically. Although other modelers have developed models of asymmetric properties around particular strike-slip faults, mine is the first published study to systematically investigate a suite of asymmetric models to address the general question of asymmetric deformation profiles. 5.1.2 Incorrect assumptions about fault geometry may bias slip rate estimates from geodetic data My study of the Point Reyes area indicates that some plate-boundary strain is occurring along the West Napa fault, which has been assumed to be inactive during the Holocene (Jennings, 1994). This result is in agreement with a recent elastic block model of the region (D‗Alessio et 100 al., 2005) as well as new geological studies (e.g. Wesling and Hanson, 2008). My analysis of new, double-difference hypocenter data from the region also shows that both the Green Valley and the Napa faults are dipping toward the southwest. The West Napa fault dip angle is 62 degrees and the Green Valley‘s is 75 degrees to the southwest. This means that their creeping extensions are offset several km from their surface traces, leading to an apparent asymmetry in the GPS velocity field around these faults. Accounting for this offset dramatically improves the fit of my models to GPS velocity data. This finding links the Point Reyes project to Chapter 2: only high substrate viscosities and creeping faults at depth offset from their surface traces were able to explain the GPS velocity field. This study shows that in general, modest departures from a vertical fault plane may complicate interpretations of GPS data for slip rates (for example, via elastic block modeling). 5.1.3 Subduction zone earthquake cycle models without backslip: still an unresolved problem For my third project, I modeled the Cascadia subduction zone with a realistic, curved interface and its interseismic slip driven by relative motion of the far-field model boundaries. The results show unrealistic, continuous uplift along the North America plate boundary. Because of this problem, models of earthquake cycle deformation for subduction zones make use of a kinematic approach (―backslip‖), in which the upper and lower plates are assumed to be nondeforming, allowing the modeler to represent interseismic deformation with specified, backward slip along the seismogenic part of the fault. The problem with this approach is that it cannot be used to model earthquake-cycle deformation with stress-driven interseismic fault zone creep. Efforts to get around this problem with specified velocities at points in the subducting plate (to force its rotation) resulted in downward motion of the slab, and widespread, steady subsidence. I developed several numerical approaches to this problem, applying a correction to the model in terms of either forces or displacements, at nodes along the subduction interface or throughout the model domain. This correction term is essentially meant to represent forces which drive the flexure of the subducting plate. The final result of this study shows that the correction term should be computed as a reaction force in the initial model when both the half-slip rate, as a 101 boundary condition, and a split-node condition are simultaneously enforced at each interface node. Currently, we know of no finite-element codes in which this is implemented. I began to develop some modules to give GeoFEST this capability, but more work is required. 5.2 Limitations The main limitations encountered throughout the course of this thesis can be attributed to: The coverage and quality of the GPS dataset. Computational restrictions. For instance, with reference to Chapter 3, there are only two campaign-mode GPS sites in the vicinity of the Napa fault and there are no continuous GPS sites nearby. There is also a big gap of data coverage to the west of the Napa fault. This is the area where we need more data with high precision to be able to locally calibrate the model. This part of the profile is strongly affected by the slip of Napa fault, so precise GPS velocity data would lead to a better understanding of this fault. There is also a data coverage problem in southwest British Columbia. There are only 12 continuous GPS sites in the area, and all but four are located on Vancouver Island, which covers less than 100 km of my model profile line. Campaign-mode GPS sites are also mainly located on Vancouver Island and just a few of them are on the Lower Mainland. There is a huge gap in the GPS surface velocity field further inland from Vancouver. Although our modeling never advanced to this point, the sparse GPS data would have made it difficult to calibrate the model. Computational restrictions refer to two issues. First, there are hardware limitations which limit the size of the models I can run, and the number of timesteps. This is usually a matter of patience and in most cases was not an issue for me (see below). The second computational restriction is the lack of software which can model the processes of interest. Both kinds of computational restrictions are experienced in this study. When running a viscoelastic finite element model, depending on the viscosity, we need to reduce time steps to as short as day (or in some nonlinear viscoelastic models, a minute) during 102 the immediate postseismic period. This lengthens the computation time to weeks in a good situation and several months for some of the models. For earthquake cycle models incorporating power-law viscous flow, simulations were very slow due to this dynamic time stepping and I ended up running just a few of such models. Despite adjusting the model and code parameters to optimize the calculation, the shortest computation time for my nonlinear models was more than two weeks. The next limitation for this study was the lack of finite element codes for geophysical purposes which could apply a particular kind of boundary condition (specified stress and/or displacement simultaneous with a split node condition at the same node). In addition, existing codes do not output the reaction force of the nodes which is required to compute the correction force at the subduction interface. 5.3 Future work Although Chapter 2 of this thesis explores a large suite of asymmetric deformation models, models incorporating more complex linear and nonlinear viscoelastic rheologies should be explored in the future. Such rheologies have been proposed to explain deformation around major faults where sufficient data exist to justify detailed postseismic and interseismic modeling (e.g., Hearn et al., 2009; Hetland and Hager, 2005). In Chapter 3 of this thesis I point out the dipping geometry of faults in the vicinity of Napa, California, and demonstrate that the West Napa fault is active. This result shows the necessity of establishing a denser network of GPS sites in the region, and more sophisticated modeling, in order to better understand earthquake hazard in this area. In general, where block models are being used to infer slip rates on faults, they need to account for non-vertical dips of strike-slip faults in the upper crust. This is more important than incorporating heterogeneous viscoelastic rheology. Earthquake cycle models for strike-slip faults are developed much more than those for subduction zones. After about 30 years of models making use of the backslip approach, it is time to model subduction zone earthquake cycles (with stress-driven fault zone creep) more 103 realistically. This requires the development and use of a correction force which is discussed completely and explicitly in Chapter 4. 104 Bibliography Aagaard, B.T., R.W. Graves, A. Rodgers, T.M. Brocher, R.W. Simpson, D. Dreger, N.A. Petersson, S.C. Larsen, S. Ma and R.C. Jachens, 2010, Ground-Motion Modeling of Hayward Fault Scenario Earthquakes, Part II: Simulation of Long-Period and Broadband Ground Motions, Bull. Seis. Soc. Am., 100, 2945-2977, doi: 10.1785/0120090379. Audet, P., M. Bostock, N. Christensen and S. Peacock, 2009, Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing, Nature, 457, 76-78. Baldwin, J.N. and J. J. Lienkaemper, 1999, Paleoseismic investigations along the Green Valley fault, Solano County, California, Unpublished report, Bay Area Paleoseismological Experiment Contract No. 98WRCN1012, 18 pp. Baldwin, J., R. Turner, and J. Lienkaemper, 2008, Earthquake record of the Green Valley Fault, Lopes Ranch, Solano County, California, Final Technical Report, U. S. Geological Survey National Earthquake Hazards Reduction Program Award Number 06HQGR0144. Bennett, R., W. Rodi and R.E. Reilinger, 1996, Global PositioningSystem constraints on fault slip rates in southern California and northern Baja, Mexico, J. Geophys. Res., 101, 21943-21960. Bostock, M. G. and J. C. Vandecar, 1995, Upper-Mantle Structure of the Northern Cascadia Subduction Zone, Canadian J. Earth Sci., 32, no. 1, 1-12. Brocher, T., 2005, Empirical relations between elastic wavespeeds and density in the Earth‘s crust, Bull. Seis. Soc. Am., 95, 2081-2092. Brossy, C., K. Kelson and M. Ticci, 2010, Final Technical Report: Digital compilation of data for the Contra Costa Shear Zone for the northern California Quaternary fault map database: Collaborative research with William Lettis and Associates and the U. S. 105 Geological Survey, USGS Earthquake Hazards Program External Research Support (NEHRP) Award Number 07HQGR0063. Bryant, W.A. (compiler), 2000, Fault number 36a, West Napa fault, Browns Valley section, in Quaternary fault and fold database of the United States: U.S. Geological Survey website, http://earthquakes.usgs.gov/regional/qfaults, accessed 04/22/2011 at 12:13 PM. Budding, K.E., Schwartz, D.P. and Oppenheimer, D.H., 1991, Slip Rate, Earthquake Recurrence, and Seismological Potential of the Rodgers Creek fault zone, Northern California: Geophys. Res. Lett., 18, 447-450. Burford, R. and P. Harsh, 1980, Slip on the San Andreas Fault in central California from alinement array surveys, Bull. Seis. Soc. Am., 70, 1233-1261. Bürgmann, 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. Chery, J., 2008, Geodetic Strain across the San Andreas Fault reflects elastic plate thickness variations (rather than fault slip rate), Earth Plan. Sci. Lett., 269, 352-365. Clowes, R. M., C. A. Zelt, J. R. Amor, and R. M. Ellis. 1995. Lithospheric structure in the southern Canadian Cordillera from a network of seismic refraction lines. Canadian J. Earth Sci., 32, no. 10, 1485-1513. Clowes, R. M., D. J. Baird and S. A. Dehler, 1997, Crustal structure of the Cascadia subduction zone, southwestern British Columbia, from potential field and seismic studies, Canadian J. Earth Sci. 34, no. 3, 317-335. Cohen, S. and M. Kramer, 1984, Crustal deformation, the earthquake cycle, and models of viscoelastic flow in the asthenosphere, Geophys. J. Royal Ast. Soc., 78, 735-750. 106 Currie, C., K. Wang and R. Hyndman, 2004, The thermal effects of steady-state slab-driven mantle flow above a subducting plate: The Cascadia subduction zone and backarc, Earth Plan. Sci. Lett., 223, 35-48. D‘Alessio, M., I. Johanson, R. Bürgmann, D. Schmidt and M. Murray, 2005, Slicing up the San Francisco Bay Area: Block kinematics and fault slip rates from GPS-derived surface velocities, J. Geophys. Res., 110, doi:10.1029/2004JB003496. Dieterich, J. H. and B. Kilgore, 1996. Implications of fault constitutive properties for earthquake prediction, Proc. Nat. Acad. Sci. USA, 93, no. 9, 3787-3794. Dixon, T., E. Norabuenal, L. Hotaling, 2003, Paleoseismology and global positioning system: earthquake-cycle effects and geodetic versus geologic fault slip rates in the Eastern California shear zone, Geology., doi:10.1130/0091-7613(2003)031. Dixon, J., T. Dixon, D. Bell and R. Malservisi, 2004, Lateral variation in upper mantle viscosity: Role of water, Earth Plan. Sci. Lett., 222, 451-467. Dragert, H., K. L. Wang and T. S. James, 2001, A silent slip event on the deeper Cascadia subduction interface. Science, 292, 1525-1528. Drew, J. J. and Clowes, R. M., 1990, A re-interpretation of the seismic structure across the active subduction zone of western Canada, In Green, A. G., ed., Studies of Laterally Heterogeneous Structures Using Seismic Refraction and Reflection Data, Geol. Surv. Can., Paper 89-13, 115-132. Fialko, Y., 2006, Interseismic strain accumulation and the earthquake potential on the southern San Andreas fault system, Nature, 441, 968-971. Finzi, Y., E. Hearn, Y. Ben Zion and V. Lyakhovsky, 2010, Structural properties and deformation patterns of evolving strike-slip faults: Numerical simulations incorporating damage rheology, PAGEOPH Topical Volumes, 1537-1573, DOI: 10.1007/978-3-0346-0138-2_2. 107 Finzi, Y., 2009, Strike-slip fault structure and fault system evolution: A numerical study applying damage rheology, UBC PhD Dissertation. Flück, P., R. D. Hyndman and K. Wang. 1997, Three-dimensional dislocation model for great earthquakes of the Cascadia subduction zone. J. Geophys. Res. 102, 20539-20550. Forsyth, D. and S. Uyeda, 1975, Relative importance of driving forces of plate motion, Geophysical Journal of the Royal Astronomical Society 43, no. 1, 163-200. Freund, L. B., and D. M. Barnett, 1976, 2-Dimensional Analysis of Surface Deformation due to Dip-Slip Faulting, Bull. Seis. Soc. Am., 66, 667-675. Fuis, G., V. Langenheim and M. Kohler, 2008, The San Andreas Fault In Southern California Has a ―Propeller‖ Shape—Implications for Tectonics and Seismic Hazard, Joint Meeting of The Geological Society of America, Soil Science Society of America, American Society of Agronomy, Crop Science Society of America, Gulf Coast Association of Geological Societies with the Gulf Coast Section of SEPM, paper 233-4. Fulton, P., G. Schmalzle, R. Harris and T. Dixon, 2010, Reconciling patterns of interseismic strain accumulation with thermal observations across the Carrizo segment of the San Andreas Fault, Earth Plan. Sci. Lett., doi:10.1016/j.epsl.2010.10.024. Gordon, R. G., 1998, The plate tectonic approximation: Plate nonrigidity, diffuse plate boundaries, and global plate reconstructions, Ann. Rev. Earth Plan. Sci., 26, 615-642. Goldfinger, C., C. H. Nelson, J. E. Johnson and the Shipboard Scientific Party, 2003, Holocene earthquake records from the Cascadia Subduction Zone and Northern San Andreas Fault based on precise dating of offshore turbidites, Ann. Rev. Earth Plan. Sci., 31, 555-577. Goldfinger, C, Morey, A., Hans Nelson, C., Gutiérrez-Pastor, J., Johnson, J.E., Karabanov, E., Chaytor, J., Eriksson, A. and the Shipboard Scientific Party, 2007, Rupture Lengths and Temporal History of Significant Earthquakes on the Offshore and North Coast Segments 108 of the Northern San Andreas Fault Based on Turbidite Stratigraphy, Earth Plan. Sci. Lett., v. 254, p. 9-27. Green, A. G., R. Clowes, C. Yorath, C. Spencer, E. Kanasewich, M. Brandon and A. Sutherland Brown, 1986, Seismic reflection imaging of the subducting Juan de Fuca Plate, Nature, 319, 210-213. Hager, B., G. Lyzenga, A. Donnellan and D. Dong, 1999, Reconciling rapid strain accumulation with deep seismogenic fault planes in the Ventura Basin, California, J. Geophys. Res., 104, 25207-25219. Hanks, T. and H. Kanamori, 1979, A Moment Magnitude Scale, J. Geophys. Res., 84, 23482350. Hashimoto, C., K. Fukui and M. Matsu‘ura, 2003, 3-D Modelling of plate interfaces and numerical simulation of long-term crustal deformation in and around Japan, Pure. Appl. Geophys., 161, 9-10, 2053-2068. Hearn, E. H., R. Burgmann and R. E. Reilinger, 2002, Dynamics of Izmit earthquake postseismic deformation and loading of the Düzce earthquake hypocenter, Bull. Seis. Soc. Am., 92, 172-193. Hearn, E., S. McClusky, S. Ergintav and R. Reilinger, Izmit earthquake postseismic deformation and the dynamics of the North Anatolian Fault Zone, J. Geophys. Res., 114, doi:10.1029/2008JB006026, 2009. Hearn, E. H., B. H. Hager, and R. E. Reilinger, 2002, Viscoelastic deformation from North Anatolian Fault Zone earthquakes and the eastern Mediterranean GPS velocity field, Geophys. Res. Lett. 29, 1549. 109 Hecker, S., D. Pantosti, D. Schwartz, J. Hamilton, L. Reidy and T. Powers, 2005, The most recent large earthquake on the Rodgers Creek Fault, San Francisco Bay Area, Bull. Seis. Soc. Am., 95, 844-860. Herring, T., 1999, Geodetic application of GPS, Proceeding of IEEE, Vol. 87, No 1, 0018-9219. Hetland, E. and B. Hager, 2005, Postseismic and interseismic displacements near a strike-slip fault: A two-dimensional theory for general linear viscoelastic rheologies, J. Geophys. Res., 110, B10401, doi:10.1029/2005JB003689. Hu, Y., 2011, Deformation processes in great subduction zone earthquake cycles, University of Victoria PhD Dissertation. Huang, W. and K. Johnson, in preparation, 2010. Hyndman, R., C. Yorath, R. Clowes and E. Davis, 1990, The northern Cascadia subduction zone at Vancouver Island: seismic structure and tectonic history, Canadian J. Earth Sci., 27, 313-329. Hyndman, R. D. and K. Wang, 1995, The rupture zone of Cascadia great earthquakes from current deformation and the thermal regime, J. Geophys. Res., 100, 22133-22154. James, T., J. Clague, K. Wang and I. Hutchinson, 2000, Postglacial rebound at the northern Cascadia subduction zone, Quat. Sci. Rev., 19, 1527-1541. Jennings C. W. and G. J. Saucedo, 1994, Fault activity map of the California and adjacent area, California Geologic Data Map Series, Map No. 6, Calif. Division of Mines and Geology. Jennings, C. W. and W. A. Bryant, 2010, 2010 Fault Activity Map of California, California Geological Survey, Geologic Data Map No. 6, Calif. Division of Mines and Geology. 110 Johnson, K., R. Burgmann and K. Larson, 2006, Frictional properties on the San Andreas Fault near Parkfield, California, inferred from models of afterslip following the 2004 earthquake, Bull. Seis. Soc. Am., 96, 5321-5338. Jolivet, R., R. Cattin, N. Chamot-Rooke, C. Laserre and G. Peltzer, 2008, Thin-plate modeling of interseismic deformation and asymmetry across the Altyn Taugh fault zone, Geophys. Res. Lett., 35, L02309, doi:10.1029/2007GL031511. Kanamori, H. and D. L. Anderson, 1975, Theoretical basis of some empirical relations in seismology, Bull. Seis. Soc. Am., 65, 1073-1095. Karato, S. and P. Wu, 1993, Rheology of the upper mantle - a synthesis, Science, 260, 771-778. Kelsey, H. M., D. C. Engebretson, C. E. Mitchell and R. L. Ticknor, 1994, Topographic form of the Coast Ranges of the Cascadia margin in relation to coastal uplift rates and plate subduction, J. Geophys. Res., 99, 12245-12255. Kelson, K. I., Simpson, G. D., Lettis, W. R. and Haraden, C. C., 1996, Holocene slip rate and recurrence of the northern Calaveras fault at Leyden Creek, eastern San Francisco Bay region, J. Geophys. Res., 101, 5961-5975. Kelson, K. I., A. R. Streig and R. D. Koehler, 2006, Timing of late Holocene paleoearthquakes on the northern San Andreas fault at the Fort Ross orchard site, Sonoma county, California, Bull. Seis. Soc. Am., 96, 1012-1028; DOI: 10.1785/0120050123. Le Pichon, X., Kreemer, C. 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/2004/B003343. Lienkaemper, J. J. and Borchardt, G., 1996, Holocene slip rate of the Hayward fault at Union City, California, J. Geophys. Res., 101, 6099-6108. 111 Lienkaemper, J., Sickler, R., Brown, Baldwin, J.N., and Turner, R., Jr., 2008, The Green Valley fault: Geologic record of four large earthquakes in the past millennium, Third Conference on Earthquake Hazards in the Eastern San Francisco Bay Area, October 22-24, 2008, abstract, Conference Program with Abstracts, p. 67. Linker, M. and J. Rice, 1997, Models of postseismic deformation and stress transfer associated with the Loma Prieta earthquake, in U.S. Geological Survey Paper 1550-D: The Loma Prieta, California Earthquake of October 17, 1989: Aftershocks and Postseismic Effects, D253–D275. Lisowski, M., J. Savage and W. Prescott, 1991, The velocity field along the San Andreas Fault in Central and Southern California, J. Geophys. Res., 96, 8369-8389. Liu, Y. J. and J. R. Rice, 2005, Aseismic slip transients emerge spontaneously in threedimensional rate and state modeling of subduction earthquake sequences, J. Geophys. Res. 110, B08307, doi:10.1029/2004JB003424. Lundgren, P., E. Hetland, Z. Liu and E. Fielding, 2006, Southern San Andreas-San Jacinto fault system slip rates estimated earthquake cycle models constrained by GPS and interferometric synthetic aperture radar observations, J. Geophys. Res., 114, B02403,doi:10.1029/2008JB005996. Malservisi, R., K. Furlong and T. Dixon, 2001, Influence of the earthquake cycle and lithospheric rheology on the dynamics of the Eastern California Shear Zone, Geophys. Res. Lett., 28, 2731-2734. Masterlark, T., 2003, Finite element model predictions of static deformation from dislocation sources in a subduction zone: Sensitivities to homogeneous, isotropic, Poisson-solid, and half-space assumptions, J. Geophys. Res., 108, 2540, doi:10.1029/2002JB002296. Mazzotti, S., X. Le Pichon, P. Henry, and S. Miyazaki, 2000, Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: An analysis of uniform elastic strain 112 accumulation in Japan constrained by permanent GPS, J. Geophys. Res., 105, 1315913177. Mazzotti, S., H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, and M. Craymer, 2003, Current tectonics of northern Cascadia from a decade of GPS measurements. J. Geophys. Res., 108, 2554, doi:10.1029/2003JB002653. Meade, B., B. Hager, S. McClusky, R. Reilinger, S. Ergintav, O. 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. Seis. Soc. Am., 92, 208-215. Meade, B., B. Hager, 2005, Block models of crustal motion in southern California constrained by GPS measurements, J. Geophys. Res., doi:10.1029/2004JB003209. Melosh, H. J. and A. Raefsky, 1981, A simple and efficient method for introducing faults into finite-element computations, Bull. Seis. Soc. Am., 71, no. 5, 1391-1400. Nicholson, T., M. Bostock, and J. F. Cassidy, 2005.,New constraints on subduction zone structure in northern Cascadia. Geophys. J. Int., 161, no. 3, 849-859. Niemi, T. M. and Hall, N. T., 1992, Late Holocene slip rate and recurrence of great earthquakes on the San Andreas fault in northern California, Geology, 20, no. 3, p. 196-198. Okada, Y., 1985, Surface deformation due to shear and tensile faults in a half-space, Bull. Seis. Soc. Am., 75, no. 4, 1135-1154. Parker, J., A. Donnellan, G. Lyzenga, J. Rundle and T. Tullis, 2003, Performance modeling codes for the QuakeSim problem solving environment, Computational Science - ICCS 2003, Lecture Notes in Computer Science, 2659, 855-862. 113 Pollitz, F. F., R. Bürgmann, and P. Bannerjee, 2006, Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth, Geophys. Res. Int., 167, 397-420. Prentice, C., Niemi, T. N. and Hall, N. T., 1991, Quaternary tectonics of the northern San Andreas fault, San Francisco Peninsula, Point Reyes, and Point Arena, California [field trip guide], California Division of Mines and Geology Special Publication, v. 109, p. 2534. Prescott, W. H., and S.-B. Yu, 1986, Geodetic Measurement of Horizontal Deformation in the Northern San Francisco Bay Region, California, J. Geophys. Res., 91(B7), 7475–7484. Reilinger, R. E., S. McClusky, M. Oral, W. King, M. TOksoz, 1997, Global positioning system measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone, J. Geophys. Res., Vol. 102, No. B5, pages 9983-9999. Reilinger, R. E., S. McClusky, P. Vernant, S. Lawrence, S. Ergintav, R. Cakmak, H. Ozener, F. Kadirov, I. Guliev, R. Stepanyan, M. Nadariya, G. Hahubia, S. Mahmoud, K. Sakr, A. ArRajehi, D. Paradissis, A. Al-Aydrus, M. Prilepin, T. Guseva, E. Evren, A. Dmitrotsa, S. V. Filikov, F. Gomez, R. Al-Ghazzi and G. Karam, 2006, GPS constraints on continental deformation in the Africa-Arabia- Eurasia continental collision zone and implications for the dynamics of plate interactions, J. Geophys. Res., 111, B05411, doi:10.1029/2005JB004051. Romanyyuk, T. V., R. Bladely, W. D. Mooney, 1998, The Cascadia subduction zone: two contrasting models of lithospheric structure, Phys. Chem. Earth, Vol. 23 pp.297-301 Rubinstein, J., D. R. Shelly and W. Ellsworth, 2010, Non-volcanic tremor: A window into the roots of fault zones, in New Frontiers in Integrated Solid Earth Sciences, International Year of Planet Earth, Cloetingh, S. and J. Negendank, eds., Springer Science+Business Media. 114 Rybicki, K., and Kasahara, K., 1977, A strike slip fault in a laterally inhomogeneous medium, Tectonophysics, 42, 127-138. Sato, K., N. Minagawa, M. Hyodo, T. Baba, T. Hori, and Y. Kaneda, 2007, Effect of elastic inhomogeneity on the surface displacements in the northeastern Japan: Based on threedimensional numerical modeling, Earth Planets Space, 59, 1083-1093. Savage, J. and R. Burford, 1973, Geodetic determination of relative plate motion in central California, J. Geophys. Res., 78, 822-845. Savage, J. and W. Prescott, 1978, Asthenosphere readjustment and the earthquake cycle, J. Geophys. Res., 83, 3369-3376. Savage, J. C. 1983, A dislocation model of strain accumulation and release at a subduction zone, J. Geophys. Res. 88, 4984-4996. Savage, J., 1990, Equivalent strike-slip earthquake cycle models in half-space and lithosphereasthenosphere Earth models, J. Geophys. Res., 95, 4873-4879. Schaff, D. P. and F. Waldhauser, 2005, Waveform cross-correlation-based differential traveltime measurements at the Northern California Seismic Network, Bull. Seism. Soc. Am., 95, 2446-2461. Schmalzle, G., T. Dixon, R. Malservisi and R. Govers, 2006, Strain accumulation across the Carrizo segment of the SanAndreas fault, California: impact of laterally varying crustal properties, J. Geophys. Res., 111 (B05403), doi:10.1029/2005/JB003843. Schofield, W and M. Breach, 2007, Engineering Surveying, sixth edition, Elsevier Ltd. Scholz, C. H., 1998, Earthquakes and friction laws, Nature, 391, 37-42. 115 Schwartz, D. P., Pantosti, D., Hecker, S., Okamura, K., Budding, K. E. and Powers, T., 1993, Late Holocene behavior and seismogenic potential of the Rodgers Creek fault zone, Sonoma County, California, California Division of Mines and Geology Special Publication 113, pp. 393-398. Segall, P., Earthquake and volcano deformation, Princeton University Press, 2010. Shelley, D., 2010, Migrating tremors illuminate complex deformation beneath the seismogenic San Andreas Fault, Nature, 463, 648-652. Simpson, G., J. Baldwin, K. Kelson and W. Lettis, 1999, Late Holocene slip rate and earthquake history for the northern Calaveras Fault at Welch Creek, eastern San Francisco Bay area, California, Bull. Seis. Soc. Am., 89, 1250-1263. Stein, R. S., 1999, The role of stress transfer in earthquake occurrence, Nature, 402, no. 6762, 605-609. Stern, R. J., 2002, Subduction zones, Rev. Geophys., 40(4), 1012, doi:10.1029/2001RG000108 Taylor, M. A. J., G. Zheng, J. Rice, W. Stuart, and R. Dmowska, 1996, Cyclic stressing and seismicity at strongly coupled subduction zones. J. Geophys. Res., 101, 8363-8381. Thatcher, W., 1983, Nonlinear strain buildup and the earthquake cycle on the San Andreas Fault, J. Geophys. Res., 88, 5893-5902. Thatcher, W. and J. B. Rundle, 1984, A viscoelastic coupling model for the cyclic deformation due to periodically repeated earthquakes at subduction zones, J. Geophys. Res., 89, 76317640. Thatcher, W., 1995, Microplate versus continuum descriptions of active tectonic deformation, J. Geophys. Res., 100, 3885-3894. 116 Thatcher, W., 2009, How the continents deform: the evidence from tectonic geodesy, Annu. Rev. Earth Planet. Sci., doi: 101146/annurev.earth.031208.100035. Thatcher, W., G. Marshall and M. Lisowski, 1997, Resolution of fault slip along the 470- kmlong rupture of the great 1906 San Francisco earthquake, J. Geophys Res., 102, 53535367. Tse, S. T. and J. R. Rice, 1986, Crustal earthquake instability in relation to the depth variation of frictional slip properties. J. Geophys. Res., 91, 9452-9472. UN reports, 24 January 2011, UN: 2010 among deadliest years for disasters, urges better preparedness. UN News Center. Waldhauser, F. and D. P. Schaff, 2008, Large-scale relocation of two decades of Northern California seismicity using cross-correlation and double-difference methods, J. Geophys. Res., 113, B08311, doi:10.1029/2007JB005479. Wang, K., 2007, Elastic and viscoelastic models of crustal deformation in subduction earthquake cycles, in The Seismogenic Zone of Subduction Thrust Faults, T. Dixon and J. C. Moore, eds., Coumbia University Press. Wang, K. and J. He, 2008, Effects of friction behaviour and geometry of subduction fault on coseismic seafloor deformation, Bull. Seis. Soc. Am., 98, 571-579. Wang, K., R. Wells, S. Mazzotti, R. Hyndman and T. Sagiya, 2003, A revised dislocation model of interseismic deformation of the Cascadia subduction zone. J. Geophys. Res., 108, 2026, doi:10.1029/2001JB001227. Wang, K. L., 1995, Coupling of tectonic loading and earthquake fault slips at subduction zones. Pure Applied Geophys., 145, no. 3-4, 537-559. Wang, K. L., H. Dragert and H. J. Melosh, 1994. Finite-element study of uplift and strain across Vancouver Island, Canadian J. Earth Sci., 31, no. 10, 1510-1522. 117 Wang, K. L., J. H. He, H. Dragert and T. S. James, 2001. Three-dimensional viscoelastic interseismic deformation model for the Cascadia subduction zone. Earth Planets Space, 53, no. 4, 295-306. Wesling, J. R. and K. L. Hanson, 2008, Final Technical Report: Mapping of the West Napa Fault Zone for input into the northern California Quaternary Fault Database, USGS Earthquake Hazards Program External Research Support (NEHRP) Award Number 05HQAG0002. Working Group on Northern California Earthquake Potential (WGNCEP), 1996, Database of potential sources for earthquakes larger than magnitude 6 in Northern California, USGS Open File Report 96-705. Working Group on California Earthquake Probabilities (WGCEP), 2007, The Uniform California Earthquake Rupture Forecast, Version 2, USGS Open File Report 2007-1437. Vaghri, A. and E. Hearn, 2011, Can Lateral Viscosity Contrasts Explain Asymmetric Interseismic Deformation around Strike-Slip faults?, accepted manuscript, Bull. Seis. Soc. Am. Williams, C. and R. McCaffrey, 2001, Stress rates in the central Cascadia subduction zone inferred from an elastic plate model, Geophys. Res. Lett., 28, 2125-2128. Zhang, H., Niemi, T. and Fumal T., 2006, A 3000-year record of earthquakes on the northern San Andreas fault at the Vedanta marsh site, Olema, California, Seis. Res. Lett., 77, p. 176 118 Appendix A Finite element modeling The finite element method (FEM), or finite element analysis (FEA), was first developed in 1943 by R. Courant (Peter Widas, 1997). This method was very expensive in its early decades and was limited to aeronautics, automotive, defense, and nuclear industries. Since the rapid decline in the cost and increase in the power of computers after the 1970's, the FEM has developed into an incredibly precise numerical method and has become a very common tool for engineering and science. The finite element method is broadly defined as a group of numerical methods for approximating the governing equations of any continuous system. It is suitable and easily adaptable for solving problems with general geometries and boundary conditions as well as different material models. However, depending on the size of the model, developing finite element programs and building finite element models can consume a considerable amount of time. We may summarize the basic steps in a linear finite element structural analysis as follows: 1) Modeling and discretization (pre-processing): In this step the structure is divided into small elements that will be easier to study. 2) Element characteristic equation: In this step each element is studied separately to derive its equilibrium of characteristic equations. 3) Assembly of element equations: In this step, all element equations are first assembled to get the global system equations. 4) Solution of global equations: In this step, we first apply the boundary conditions of the global structure then solve the global equations. 5) Post-processing (calculation of other solution variables): In this step we compute other parameters that we need but that do not come out directly from solving global equations. For example, we obtain nodal displacements but may need displacements and velocities at particular locations (not just nodal coordinates), or stresses (obtained from spatial derivatives of modeled displacements). 119 In the following simple example of axial members (Figure A.1), there are two members which are can be considered as spring elements (here, represented with a node on each end). The elements are just moving in one direction so they have one degree of freedom (DOF) per node. A tensional force is applied at the left end, and the right end is fixed. Figure A.1. Simple model with two axial members. In the first step we need to discretize the modeled object into two elements. Element 1 on the left has a larger stiffness (k1), and has nodes 1 and 2 on its left and right ends. A tensional force is applied at node 1. Element 2 has a smaller stiffness (k2) and is defined with nodes 2 and 3, at its left and right ends. Node 3 is fixed and does not move. Figure A.2. Elements are separated to show the nodes of each elements and forces on each node of the elements. To obtain the characteristic equation for element 1, we have: 120 AE u 2 u1 k1 u 2 u1 L F11 F21 k1 u1 u 2 F21 (A.1) where F is force, E is Young's modulus, L is the length of element, A is the cross-sectional area of the element, u2 and u1 are the displacements at nodes 2 and 1, and k1 is the stiffness of the element 1. In this one-dimensional case, the stiffness of the element depends on the shape of the element and it is a function of element‘s length and cross-section. The superscript and subscript of the F respectively indicate the element and node numbers. In matrix form, equation (A.1) can be written as: k1 k 1 k1 u1 F11 k1 u 2 F21 (A.2) For element 2 we have this characteristic equation: k2 k 2 k 2 u 2 F22 k 2 u 3 F32 (A.3) To assemble the element equations, we need to know the global DOF. In this example we have 3 nodes with one DOF per node, so the global DOF is 3. This means that the global stiffness matrix k is a 3x3 matrix k 33u31 F31 (A.4) The right hand side (RHS) of the global equation is made up of the global forces on each node, it means F1=-P1, F2=0 and F3=R3 (Figure A.2). U is the displacement of each node, and the ij component of k corresponds to the uj in equation/equations belong to Fi so the global equation k1 k 1 0 k1 k1 k 2 k2 0 u1 P1 k 2 u 2 0 k 2 u 3 R3 (A.5) 121 To solve the global system of equations, we need to apply boundary conditions. In this example, the displacement of node 3 (u3) is known so the equation corresponding to it is omitted. k1 k 1 k1 u1 P1 k1 k 2 u 2 k 2 u 3 (A.6) Because node 3 is fixed (u3=0), the final system of equations may be written as: k1 k 1 k1 u1 P1 k1 k 2 u 2 0 (A.7) Now, by finding the inverse of the stiffness matrix (k), the equation is solved. Fk 1 u (A.8) Finally, in the post processing step we compute the parameters which we need, but which do not come out of the global solution. For instance, if we need tension at node 3 we need to compute: R3 k 2 u 2 (A.9) Generally, in any kind of FEM problem, the most important task is to establish the element characteristic equation (step 2 above) and the stiffness matrix (k). k e ue f e (A.10) e means these values belong to the nodes of one element. For solid mechanics problems, in the first step we need to express the displacement inside the element in terms of nodal displacements. This is made possible by shape functions. To find shape functions, we assume the displacement field within element is a linear combination of xiyjzk where i, j and k can be 0,1,… depending on the kind of element. is coefficient for the each part of linear combination ux, y m x i y j z k (A.11) 122 By replacing nodal displacements in equation (A.12), ‘s are obtained in terms of nodal displacement. Replacing ‘s in equation (A.12) gives shape functions, u( x, y, z ) 3,1 N ( x, y, z ) 3,1 u ej1 (A.12) ‗i‘ is the dimension of the model and ‗j‘ is dimension of model times number of nodes in the element. The next step is to find strain in terms of nodal displacements: u ( x, y, z ) N ( x, y, z ) e u x x u ( x, y, z ) N ( x, y, z ) e u y y u ( x, y, z ) N ( x, y, z ) e u z z u ( x, y, z ) u ( x, y, z ) N ( x, y, z ) N ( x, y, z ) e u x y x y xx yy zz xy xz u ( x, y, z ) u ( x, y, z ) N ( x, y, z ) N ( x, y, z ) e u x z x z yz u ( x, y, z ) u ( x, y, z ) N ( x, y, z ) N ( x, y, z ) e u y z y z (A.13) In the matrix equation we will have: x, y, z 6,1 Bx, y, z 6, j u ej,1 (A.14) We use the constitutive equation (Hooke‘s Law) to get stress ( in terms of nodal displacements. C x, y, z 6,1 C6,6 Bx, y, z 6, j u ej,1 (A.15) C is the elastic stiffness tensor, which is different from the stiffness matrix for finite element modeling k. Finally we use a governing equation to get the element characteristic equation. Here at this point we use the virtual work equation for infinitesimal deformations, which is: 123 dv u bdv u tds T T v T v (A.16) s b is body force and t is traction. By substituting values from equations (A.12), (A.14) and (A.15), equation (A.16) may be written as: u B C B u dv u N tds u N bdv e T e T e e e v T T e s T v T (A.17) ue is nodal displacement and ue is constant, so it is omitted from both sides of the equation: B C B dvu N tds N bdv e T v e C B dv : k B e T T e e T s v Generalized stiffness matrix v f b N bdv : Body T force vector (A.18) v f N tds : Traction T t force vector s [k ]{u e } { f t } { f b } The following example is for a triangular element, for a plane strain problem. For a triangular element we assume the displacement field within the element is: u x, y 11 21 x 31 y vx, y 12 22 x 23 y (A.19) u is displacement in the x direction and v is displacement in the y direction. Replacing nodal displacements in the equation we will have a set of 6 equations. We will have the values of in terms of nodal displacements. Replacing these values in equations (A.19) we will get shape functions: 124 0 N 2 x, y 0 N 3 x, y u x, y N1 x, y N 1 x, y 0 N 2 x, y 0 vx, y 0 N i x, y u1 v 1 0 u 2 N 3 x, y v 2 u 3 v3 1 ai bi x ci y 2 (A.20) (A.21) Where, a1 x 2 y3 x3 y 2 b1 y 2 y3 c1 x3 x 2 a 2 x3 y1 x1 y3 b2 y3 y1 c 2 x1 x3 a3 x1 y 2 x 2 y1 b3 y1 y 2 c3 x 2 x1 (A.22) And area of the triangular element 1 x1 1 1 x2 2 1 x3 y1 y 2 area of the triangula r element y3 (A.23) In the next step we need to compute the B matrix. For this example we need xx, yy and xy from equation (A.13). We get B with the use of equation (A.21): b1 0 b2 1 B 0 c1 0 2 c1 b1 c 2 e 0 c2 b2 b3 0 c3 0 c3 b3 (A.24) For plain strain, the elastic stiffness tensor is 3x3: 125 1 0 C E 2 1 0 1 1 0 0 2 (A.25) By substituting C, B and N matrix in equation (A.18) we get the stiffness matrix and force vectors. k BT C B ds BT C B s (A.26) Because this a 2D example, the body integral will change to a surface and it is the area of the element. Depending on how the force is applied to the element, we can obtain force vector from equations (A.18). To model time-dependent deformation, U or F may be updated for each time step and the system then re-solved. For instance, to model deformation in response to viscoelastic relaxation or afterslip, creep per time step is modeled and displacements U are updated. This results in updated element stresses, which are used to compute creep or aseismic fault slip rates for the subsequent time step. For afterslip, the split node displacement must be updated (and the force vector changed) following each time step. 126 Appendix B Mesh software for automatic mesh generation The first step of finite element modeling is to define the geometry of the model and mesh it; that is, to discretize the model into small elements. Finite element software for geophysical applications requires geometry (nodal coordinates and a mapping of nodes to model elements) as well as element material properties, boundary conditions, fault information, time-stepping and output specifications, and some other information. Only some commercial FE software can fully automatically mesh the model domain, but none can do so for 3D models with subdomains, or for 3D model domains cut by split-node faults. This is challenging because the nodal coordinates must match at the boundaries of adjoining model subdomains (or along a split-node fault surface). For a small model with a limited number of elements, it is usually straightforward to generate elements and other information for input files by hand. For models with a large number of nodes and elements, or with a complex 3D geometry, we need to mesh the model automatically. Otherwise the process is time consuming and it is difficult to find meshing mistakes (which are unavoidable). The most important concepts of meshing are: Each zone should be meshed independently with a pre-determined element size range Element size should increase gradually with distance from the fault Aspect ratios of elements should not be large: meshes with equant elements yield the most numerically stable system of equations. I wrote the program ―Mesh‖ in Matlab to generate meshes and generate input files for GeoFEST and/or Tecton, while satisfying the above conditions. The software is flexible and can generate input files for other FE codes by adding a file which contains the input file format information. We can also get particular information like nodal coordinates or boundary conditions as an output. 127 This software meshes a 2D model domain using triangle or quadrilateral elements. It discretizes the model to stablest possible elements, that is, triangles and quadrilaterals which are as close as possible to equilateral triangles and squares. Then it checks if the generated elements is the best possible one, if not it changes the element to better one. The element size can be changed from one point inside a model subdomain to another, so the software adjusts the element size in any part of the model. For example if we have a simple triangle model with 3 vertex 1, 2 and 3; we can ask the model to discretize the model such as we have 2 m size element at vertex 1, 6 m at vertex 2 and 8 m at vertex 3. The software will increase size of elements gradually from 2 at vertex 1 to 6 at vertex 2 and same in the direction of vertexes 1 to 3 and vertexes 2 to 3. The software can also mesh a 3D model domain using hexahedral elements. The size of element in any direction of point can be changed, we can request different element size in x, y and z direction in a point. It means that in a zone element size can vary from one location to another. This software contains 34 files (more than 8000 lines of MATLAB code). The software has graphical outputs which illustrate the meshed model, zones (subdomains) and their properties, boundary conditions, split and slippery nodes, and output nodes. It also displays element stresses and nodal displacements. Mesh also generates input files for GeoFEST and TECTON, enabling me to change the geometry and parameters of large models easily. Mesh requires an input file which contains data about the geometry and material properties of the model, as well as the FEM input filenames and addresses. It is very simple to edit the model and generate new FEM input files, so Mesh is a very helpful tool for investigating model sensitivity to variations in model parameters and geometry. See Appendix D for input file information and formats. The program Mesh has 4 parts. The first part is the main directory, which contains the eight files that are common to both the 2D and 3D versions of the code. The second part is the ―2Dimension‖ directory, containing code which is used to mesh 2D models. This directory contains 12 files and meshes the model into triangles and/or rectangles with different sizes for each zone. The third part is the ―3Dimension‖ directory. This directory contains 13 files of code for meshing simple 3D models. The code divides volumes into quadrilaterals. The elements size can 128 grow or shrink from one side to the other of each modeled subdomain. This part of the code is still under development. The fourth directory is the data directory, which can take any name which is mentioned in the main file. This directory contains the input file which has the initial information about the model. Output files (i.e., input files for GeoFEST and Tecton) are also saved in this directory. The files and directories are listed and briefly described below. MAIN DIRECTORY FILES: Main: This function is for entering the data file name and its directory. It is also for determining the graphical outputs (48 lines) F_Mesh: This function calls F_data then transfers data to the 2D or 3D directory to mesh the model. (96 lines) F_Data: This function reads and checks data from the input file. (239 lines) F_Geofest: This function prepares an input file for GeoFest. (514 lines) F_Tecton: This function prepares an input file for Tecton. (1221 lines) F_Qube: Result_Graf: This function draws cubes. (39 lines) This function displays output from GeoFEST or Tecton. (1309 lines) Result_Graf3D: This function displays output from 3D FE models. (384 lines) 2 DIMENSION FILES: F_Angle: F_Diamond: This function computes the angle between 3 points. (44 lines) This is version 3.1 of function F_Diamond. This function computes the nodes and elements inside a diamond. (398 lines) 129 F_Divider: This function find if the element size at the ends of line is diferent or not. If they are different, the line is divided gradually from fine mesh to coarse mesh. If not, it is divided into constant segments - almost equal to the size of elements on that zone. (65 lines) F_Element: This function computes elements between two lines. The elements are triangular or rectangular shape. The outputs of this function are f_Ele (number of elements) and elem (matrix of elements) as global parameters: (169 lines) F_Mesh2D: F_Node: This function meshes a 2D model. (1042 lines) This function adds the number of nodes for each point of matrix of lines. (30 lines) F_Pos_Point: This function divides a line into segments (creates a discretized line), and adds new nodes to the Node matrix. (65 lines) F_RegNode: This function registers any position (x, y) as a node. (20 lines) F_Tri_RotDir: This function determines the direction of rotation from one line to another one. (22 lines) F_Triangle: This function computes the nodes and elements inside a triangle. (142 lines) F_Zone: This function divides partitions into triangles and rectangles, and prepares it for meshing. (765 lines) ZonePlot: This function plots the elements and/or nodes in any desired zone. (76 lines) 3 DIMENSION FILES: F_LineMaker: This function divides the line into constant number of segments equal to the size of elements in the zone. (116 lines) F_Divider_Like: This function divides the lines proportionally like the other lines. (47 lines) 130 F_Divider_Seg: This function divides the lines gradually from fine mesh to coarse mesh. If the segment size is the same at both ends, it divides the line into constant segments –almost equal to the dimensions of elements in that zone. (50 lines) F_Element_Cube: This function makes elements of the volume. The volume is a divided cubes and elements are quadrilaterals. f_Ele (number of elements) and elem (matrix of elements) are the global parameters. (45 lines) F_FindDivid: This function finds if any end of line connects to other zone. If it is connected to other zone, line is divided gradually from fine mesh to coarse mesh. If not, it is divided into constant segment -almost equal to the size of elements on that zone. (82 lines) F_LineMaker: This function divided line between two points to appropriate number of segments with appropriate size. (116 lines) F_Mesh3D: F_Node3D: This function meshes a 3D model. (641 lines) This function divides line to segment and creates a discreatized lines and adds new nodes to Node matrix. (29 lines) F_Pos_Pn3D: This function determines the position of 3D points of divided line. It uses vector V_P which shows proportionally how far is the inside points from first point. (58 lines) F_RegNode: This function registers any position (x, y) as a node. (20 lines) F_SurfMaker: This function divides surface surrounded between 4 lines and registers the nodes on the surface. (34 lines) F_Surface: This function meshes the surface between 4 lines if it has not been meshed yet. (137 lines) F_VolMaker: This function meshes the 3D model. (30 lines) 131 DATA FILE: F_Cascadia: This function computes and partitions zones of the Cascadia profile model, based on Lithoprobe data. (623 lines). 132 Appendix C Lagrange multiplier The Lagrange multiplier method is a technique for applying constraints in optimization problems. We can use this method to constrain the finite element solution. To apply this method first of all we need to identify the constraint equation which is a function of displacement like g(u)=0. This constraint equation imposes particular condition to the solution. The finite element global equation is: K u F f u K u F (B.1) Where K is stiffness matrix, u is the displacement vector and F is the force vector. So our objective is: min f (u ) subject to g (u ) 0 min L where L f (u ) g (u ) (B.2) (the Lagrange multiplier) and u are unknowns. G(u) is the constrain equation for instance if we want to force the node number 2 to move 2 cm then g(u) will be ( ) In general g(u) is a set of equation system. L u 0 min L L 0 Where f g u u g u 0 ( A) (B.3) ( B) f is the global equation matrix. So equation (A) is: u K u F g 0 u (B.4) For a linear constraint equation we will have: g (u) a1u1 a2 u 2 an u n C 0 (B.5) 133 So the equation (B) is: g (u ) g g g u1 u2 un C 0 u 1 u 2 u n (B.6) The new equation system is: g u1 k nn g u 2 g u1 u F g 1 1 u 2 u 2 F2 g u F u n n n C g 0 u n (B.7) The solution of this system will satisfy the condition of the constraint equation. 134 Appendix D Input file for mesh generating code (Mesh) We need to define the model for Mesh, using a Matlab file (.m file). The following command is inserted into the top of the input file. function [InOutFile,Title,dim, Points, Blocks, INF, BC, Velocity, Winkler, Split, Slippery,... CorVelo,Interface, Continent, CorBC, CorSplit, OutNode, V_OutElem, OutPar,... Es_Mod]=Filename (INF, Velocity, Winkler, Split, Slippery, CorVelo, Interface,… Continent, CorBC, CorSplit, OutNode, V_OutElem, Es_Mod) Input output file name InOutFile determine the GeoFEST input file name with the extension Gin, and the GeoFEST output file name with the extension Gou. InOutFile = 'PointReyesV21-20'; Title ―Title‖ is the title of the project, which describes the project and gives some information for future use. Title = ['Point Reyes profile, 3 strike slip fault. This model is for investigating the viscosity of LC and mantle'] Dimension of space dim = 3; Points The matrix ―points‖ contains the points that are used to define the different model zones (that is, model subregions bounded by a contrast in material properties or a fault). These points define the edges and corners of zones. The first columns of the matrix are: 135 Points = [ # of point, x of point , y of point, z of point] The fifth column which will be added by the program is node number of each point. Its structure is: Points = [ # of point, x of point ,y of point, z of point, # of node] One can add columns 6, 7 and 8 to indicate preferred size of elements at each point in the x, y, and z directions. These data are obtained from the ―F_DataIn‖ file. Points = [ 111, 112, 113, 114, 115, ]; 50000, 20000, 0, 0, 50000, 80000, 20000, 0, 0, 8000, 90000, 40000, 0, 0, 2000, 40000, 40000, 0, 0, 8000, 10000, 50000, 0, 0, 2000, 1000, 0; 1000, 0; 1000, 2000; 1000, 0; 1000, 0; Blocks First column: Number of the block Second column: Points which define the block, numbered as indicated on figure D1: [1,2,3,4,5,6,7,8] Figure D.1. This figure shows the shape of a zone with its vertex points. Third column: P-wave velocity in the block (m/s) 136 Fourth column: Density of the material of the block (kg/m^3) Fifth column: Viscosity of the block (Pa s, use zero for an elastic material) Sixth column: The power-law exponent for viscosity, use 1 if the block is a Maxwell linear viscoelastic material. Seventh column: The desired size of model elements which will be used to discretize the block. The numbers in the vector indicate element size respectively in x, y and z direction. Blocks= { 1, [318,319,329,328,118,119,129,128], 5091, 2549, 2, [3171,318,328,3271,1171,118,128,1271],5091, 2549, 3, [217,2171,2271,227,117,1171,1271,127],5091, 2549, 4, [2169,217,227,2269,1169,117,127,1269],5732, 2662, 5, [316,3169,3269,326,116,1169,1269,126],5732, 2662, }; 0, 0, [2000,2000,2000]; 0, 0, [2000,2000,2000]; 0, 0, [2000,2000,2000]; 0, 0, [2000,2000,2000]; 0, 0, [2000,2000,2000]; Displacement Boundary Condition First column: Boundary condition ID number. Second column: This parameter is the node number(s) for the boundary condition. For imposed displacement at a point, it is one node number. For a displacement boundary condition along a line, two node numbers are given (one node at each end of the line). For a quadrilateral surface, numbers for the four nodes at the corners are given (as in the example below). Third column: Vector with three values representing x, y, and z degrees of freedom for the boundary condition. 0 indicates that the boundary is fixed, and 1 indicates that it is free to move, in the relevant direction. Fourth column: The value of BC displacement (in meters) for the x and y directions. If the nodes are fixed only in one direction, this element only will have one member. BC = { 1, [318, 319, 119, 118], 2, [3171, 318, 118, 1171], }; [0,1,0], [0,0,0]; [0,1,0], [0,0,0]; 137 Velocity Boundary Condition First column: Boundary condition ID number Second column: This parameter is the node number(s) for the velocity boundary condition. For imposed velocity at a point, it is one node number. For a velocity boundary condition along a line, two node numbers are given (one node at each end of the line). For a quadrilateral surface, numbers for the four nodes at the corners are given (as in the example below). Third column: Vector with three values representing x, y, and z degrees of freedom for the boundary condition. 0 indicates that the boundary is fixed, and 1 indicates that it is free to move, in the relevant direction. Fourth column: Value of velocity for x and y directions (m/s). If have not fixed one of directions -x or y- so it does not matter what value you put for it. BC is applied to model in the first process but velocity is applied to the model every time step. Velocity = { 1, [319,119,129,329], 2, [311,111,121,321], }; [0,0,0], [0,-V/2,0]; [0,0,0], [0,V/2,0]; Split nodes First column: ID number of fault strand Second column: This parameter is the node number(s) for the split-node fault. For slip at a point, it is one node number. For slip along a line, two node numbers are given (one node at each end of the line). For slip along a quadrilateral surface, numbers for the four nodes at the corners are given (as in the example below). Third column: Slip amount (split) at each point Fourth column: B vector (Figure D.1). Fifth column: S vector (Figure D.1). 138 Sixth column: Repeat time (s). Seventh column: Start time (s). Figure D.2. This figure illustrates S and B vector at the fault (split) nodes. Split={1, [113,213,223,123], [1.2,1.2,1.5,1.5], [0,0,-1], [0,1,0], 200, 0; 2, [115,215], [09,1.2], [0,0,-1], [0,1,0], 200, 0; }; Output nodes First column: Second column: Output node, line, or surface ID number This parameter is the node number(s) for model output. For output at a point, it is one node number. For output at all nodes along a line, two node numbers are given (one node at each end of the line). For output at all nodes falling along within a quadrilateral surface, numbers for the four nodes at the corners are given (as in the example below). OutNode = { 1, [118,119,129,128]; 2, [1171,118,128,1271]; 3, [117,1171,1271,127]; 139 4, 5, }; [1169,117,127,1269]; [116,1169,1269,126]; Time step and output time information Number of scheduled earthquake cycle of output times. OutCyc = [11,12,77,78]; Number of scheduled times to print output in each cycle. All times are with respect to the time of the earthquake. OutTime = [0.01,0.1,1,2,3,4,5,10,15,20,40,70,100]; a: use present time step size until this time is reached (in the time units of the problem) b: implicit/explicit parameter; 0<alpha<1; usually set to 1.0. (dimensionless) c: time step size (are in the time units of the problem) a = [0.1;1;30;200]; b = [1;1;1;1]; c = [0.01;0.1;1;10]; The rest of parameter are adjusted by the code (please don’t change them) Eq_Prd = QCS; Period of earthquake cycle Number of scheduled times for each cycle to print output; this is the number of lines to follow in the block immediately below. TimePar = [a,b,c]; Timestep parameters: repeat this line for each time group e.g.: (a) (b) (c) 140 turnOn = 1; Turn-on time for the boundary velocities in the block immediately following enter value in units of physical (simulation) lime. Earthquake and backup parameters: [a,b,c] a = length(TimePar(:,1)); Number of time groups: number of different intervals of time step size. This number of lines will follow in the next input block. b = 1; Reform steps: number of time steps to take before reforming the stiffness matrix. The standard Hughes implicit algorithm calls for this to be 1. When using the iterative conjugate gradient solver, there is no advantage to setting time savings maybe realized by setting it to a modest number such as 5, but this is a heuristic shortcut that is without rigorous justification. The stability of the time-stepping ordinarily guaranteed by the implicit scheme may not be assured if this is set to a number greater than 1. c = 200; Backup steps: number of time steps to take before saving the simulation state to disk as backup or restart. d = length([Split{:,1}]); Number of fault strand. No_Prd is the number of earthquake perild we want to model and Eq_Prd is the period of earthquake cycle. qkPar = [a,b,c,d,No_Prd,Eq_Prd]; RNod = 9; Total number of nodes to report on in output file. Set zero if no nodal output is requested and set to -1 to output every node in the grid and other number if you need results for particular nodes in output file (It does not matter what number you input in this case. The code will calculate number of nodes). REle = -1; Total number of elements to report in output file. Set zero if no element output is requested and set to -1 to output every element in the grid. LastQuake = This parameter shows when the lat earthquake happened. strand = [Split{:,1};Split{:,6};Split{:,7}]'; 141 OutPar ={OutTime,TimePar,qkPar,RNod,REle,strand,turnOn,OutCyc,LastQuake}; All these parameters are computed above to form the output parameters. 142
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Viscoelastic finite-element models of the earthquake...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Viscoelastic finite-element models of the earthquake cycle along plate-boundary faults Vaghri, Ali 2011
pdf
Page Metadata
Item Metadata
Title | Viscoelastic finite-element models of the earthquake cycle along plate-boundary faults |
Creator |
Vaghri, Ali |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | In this thesis I describe three related projects on the mechanics of plate bounding faults and fault systems. All of the projects involve viscoelastic earthquake cycle modeling, using the finite element method, and investigating how surface deformation is affected by viscoelastic relaxation or viscous shear zone creep. In the first project I develop test models to study how contrasts in material properties across a fault can affect surface deformation. I model contrasts in elastic plate thickness and in viscosity profiles across a strike-slip fault. I find that an offset between the surface trace of a fault and its position at depth (where it creeps), due to a non-vertical dip, may make interseismic deformation look asymmetric. Rheological contrasts do not yield dramatic asymmetry for models with realistic viscosity profiles and the sense of asymmetry in viscoelastic models actually reverses during the interseismic interval. In the second project, I model a set of parallel strike-slip faults, which allow relative motion of the North American and Pacific Plates in northern California. Geologic information allows me to take the timing of large earthquakes on these faults into account in the modeling, and I compare modeled surface deformation with GPS data from the region. I find that to explain the GPS velocity field, I must incorporate the dip of the faults into my model. The required dip for the Green Valley Fault in particular is consistent with new, double-difference hypocenter data from the region. The final project involves modeling Cascadia Subduction Zone earthquake cycles. This subduction zone has been modeled in the past, but always using a kinematic device (“backslip”) for modeling slip on the subduction interface. I developed a detailed profile model of the Cascadia Subduction Zone, based on data from the Lithoprobe project. Since the Subduction Zone has a curved profile and its slip is driven by relative motion of the model boundaries, the models produced unrealistic uplift along the North America plate boundary. I tested several approaches for preserving the geometry of the subduction interface. I suggest that a new implementation of split nodes in the finite-element method would be required to make this correction work. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-09-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution 3.0 Unported |
DOI | 10.14288/1.0053367 |
URI | http://hdl.handle.net/2429/37136 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by/3.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2011_fall_vaghri_ali.pdf [ 4.84MB ]
- Metadata
- JSON: 24-1.0053367.json
- JSON-LD: 24-1.0053367-ld.json
- RDF/XML (Pretty): 24-1.0053367-rdf.xml
- RDF/JSON: 24-1.0053367-rdf.json
- Turtle: 24-1.0053367-turtle.txt
- N-Triples: 24-1.0053367-rdf-ntriples.txt
- Original Record: 24-1.0053367-source.json
- Full Text
- 24-1.0053367-fulltext.txt
- Citation
- 24-1.0053367.ris
Full Text
Cite
Citation Scheme:
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>
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-0053367/manifest