Yehuda Ben-Zion

Elastodynamic Analysis for Slow Tectonic Loading with Spontaneous Rupture Episodes on Faults with Rate- and State-dependent Friction

We present an efficient and rigorous numerical procedure for calculating the elastodynamic response of a fault subjected to slow tectonic loading processes of long duration within which there are episodes of rapid earthquake failure. This is done for a general class of rate- and state-dependent friction laws with positive direct velocity effect. The algorithm allows us to treat accurately, within a single computational procedure, loading intervals of thousands of years and to calculate, for each earthquake episode, initially aseismic accelerating slip prior to dynamic rupture, the rupture propagation itself, rapid post seismic deformation which follows, and also ongoing creep slippage throughout the loading period in velocity-strengthening fault regions. The methodology is presented using the two-dimensional (2-D) antiplane spectral formulation and can be readily extended to the 2-D in-plane and 3-D spectral formulations and, with certain modifications, to the space-time boundary integral formulations as well as to their discretized development using finite difference or finite element methods. The methodology can be used to address a number of important issues, such as fault operation under low overall stress, interaction of dynamic rupture propagation with pore pressure development, patterns of rupture propagation in events nucleated naturally as a part of a sequence, the earthquake nucleation process, earthquake sequences on faults with heterogeneous frictional properties and/or normal stress, and others. The procedure is illustrated for a 2-D crustal strike-slip fault model with depth-variable properties. For lower values of the state-evolution distance of the friction law, small events appear. The nucleation phases of the small and large events are very similar, suggesting that the size of an event is determined by the conditions on the fault segments the event is propagating into rather than by the nucleation process itself. We demonstrate the importance of incorporating slow tectonic loading with elastodynamics by evaluating two simplified approaches, one with the slow tectonic loading but no wave effects and the other with all dynamic effects included but much higher loading rate.

Characterization of Fault Zones

There are currently three major competing views on the essential geometrical, mechanical, and mathematical nature of faults. The standard view is that faults are (possibly segmented and heterogeneous) Euclidean zones in a continuum solid. The continuum-Euclidean view is supported by seismic, gravity, and electromagnetic imaging studies; by successful modeling of observed seismic radiation, geodetic data, and changes in seismicity patterns; by detailed field studies of earthquake rupture zones and exhumed faults; and by recent high resolution hypocenter distributions along several faults. The second view focuses on granular aspects of fault structures and deformation fields. The granular view is supported by observations of rock particles in fault zone gouge; by studies of block rotations and the mosaic structure of the lithosphere (which includes the overall geometry of plate tectonics); by concentration of deformation signals along block boundaries; by correlation of seismicity patterns on scales several times larger than those compatible with a continuum framework; and by strongly heterogeneous wave propagation effects on the earth’s surface. The third view is that faults are fractal objects with rough surfaces and branching geometry. The fractal view is supported by some statistical analysis of regional hypocenter locations; by long-range correlation of various measurements in geophysical boreholes; by the fact that observed power-law statistics of earthquakes are compatible with an underlying scale-invariant geometrical structure; by geometrical analysis of fault traces at the earth’s surface; and by measurements of joint and fault surfaces topography.

Wrinkle‐like slip pulse on a fault between different materials

Pulses of slip velocity can propagate on a planar interface governed by a constant coefficient of friction, where the interface separates different elastic materials. Such pulses have been found in two-dimensional plane strain finite difference calculations of slip on a fault between elastic media with wave speeds differing by 20%. The self-sustaining propagation of the slip pulse arises from interaction between normal and tangential deformation that exists only with a material contrast. These calculations confirm the prediction of Weertman [1980] that a dislocation propagating steadily along a material interface has a tensile change of normal traction with the same pulse shape as slip velocity. The self-sustaining pulse is associated with a rapid transition from a head wave traveling along the interface with the S wave speed of the faster material, to an opposite polarity body wave traveling with the slower S speed. Slip occurs during the reversal of normal particle velocity. The pulse can propagate in a region with constant coefficient of friction and an initial stress state below the frictional criterion. Propagation occurs in only one direction, the direction of slip in the more compliant medium, with rupture velocity near the slower S wave speed. Displacement is larger in the softer medium, which is displaced away from the fault during the passage of the slip pulse. Motion is analogous to a propagating wrinkle in a carpet. The amplitude of slip remains approximately constant during propagation, but the pulse width decreases and the amplitudes of slip velocity and stress change increase. The tensile change of normal traction increases until absolute normal traction reaches zero. The pulse can be generated as a secondary effect of a drop of shear stress in an asperity. The pulse shape is unstable, and the initial slip pulse can change during propagation into a collection of sharper pulses. Such a pulse enables slip to occur with little loss of energy to friction, while at the same time increasing irregularity of stress and slip at the source.

Dynamic ruptures in recent models of earthquake faults

Abstract We discuss several problems of dynamic rupture relevant to mechanics of earthquake faults, material sciences, and physics of spatially extended dissipative systems. The problems include dynamic rupture along an interface separating different elastic solids, dynamic rupture on a planar surface governed by strongly velocity-weakening friction, and elastodynamic calculations of long deformation history on a smooth fault in an elastic continuum. These separate problems share a number of methodological and conceptual issues that form recurring themes in the paper. An important methodological issue for computational schemes is dependency of numerical results on the used grid size. This arises inevitably in computer simulations when the assumed constitutive laws do not include a length scale (e.g., of shear or extensional displacement) over which material properties evolve. Such simulations do not have a stable underlying solution, to which they may converge with sufficient grid refinement. However, they may provide rough approximations—lacking at present a rigorous foundation—to the behavior of systems containing elements of discreteness (associated with abrupt fluctuations) at scales relevant to observations of interest. Related important conceptual issues are connections between, or when appropriate separation of, small scale phenomena (e.g., nucleation of rupture, processes at rupture front) and large scale features of the response (e.g., overall space–time dimensions of rupture, statistics of many events). Additional recurring conceptual topics are crack vs. pulse modes of dynamic rupture, the stress under which earthquake faults slip, and the origin of spatio-temporal complexities of earthquakes. These seemingly different issues probably have one or more common origins. Dynamic rupture on an interface between different solids, strongly velocity-weakening friction on a homogeneous fault, and strong fault zone heterogeneities can all produce narrow self-healing slip pulses with low dynamic stress (and low associated frictional heat) during the active part of slip. Strong fault heterogeneities probably play the dominant role in producing the observed earthquake complexities. Improved understanding of the discussed problems will require establishing connections between discrete and continuum descriptions of mechanical failure processes, generalization of current models to realistic three-dimensional dynamic models, and high-resolution laboratory and in-situ observations over broad scales of space and time. These challenging problems provide by their subject matter and involved great difficulties important targets for multi-disciplinary research by engineers, earth scientists, and physicists.

Distributed damage, faulting, and friction

We present a formulation for mechanical modeling of geological processes in the seismogenic crust using damage rheology. The seismogenic layer is treated as an elastic medium where distributed damage, modifying the elastic stiffness, evolves as a function of the deformation history. The model damage rheology is based on thermodynamic principles and fundamental observations of rock deformation. The theoretical analysis leads to a kinetic equation for damage evolution having two principal coefficients. The first is a criterion for the transition between strength degradation and recovering (healing), and is related to friction. The second is a rate coefficient of damage evolution which can have different values or functional forms for positive (degradation) and negative (healing) evolution. We constrain these coefficients by fitting model predictions to laboratory data, including coefficient of friction in sawcut setting, intact strength in fracture experiments, first yielding in faulting experiments under three-dimensional strain, onset and evolution of acoustic emission, and dynamic instability. The model damage rheology accounts for many realistic features of three-dimensional deformation fields associated with an earthquake cycle. These include aseismic deformation, gradual strength degradation, development of process zones and branching faults around high-damage areas, strain localization, brittle failure, and state dependent friction. Some properties of the model damage rheology (e.g., cyclic stick-slip behavior with possible accompanying creep) are illustrated with simplified analytical results. The developments of the paper provide an internally consistent framework for simulating long histories of crustal deformation, and studying the coupled evolution of regional earthquakes and faults. This is done in a follow up work.

Slip patterns and earthquake populations along different classes of faults in elastic solids

Numerical simulations of slip instabilities on a vertical strike-slip fault in an elastic half-space are performed for various models belonging to two different categories. The first category consists of inherently discrete cellular fault models. Such are used to represent fault systems made of segments (modeled by numerical cells) that can fail independently of one another. Their quasi-independence is assumed to provide an approximate representation of strong fault heterogeneity, due to geometric or material property disorder, that can arrest ruptures at segment boundaries. The second category consists of models having a well-defined continuum limit. These involve a fault governed by rate- and state-dependent friction and are used to evaluate what types of property heterogeneity could lead to the quasi-independent behavior of neighboring fault segments assumed in the first category. The cases examined include models of a cellular fault subjected to various complex spatial distributions of static to kinetic strength drops, and models incorporating rate- and state-dependent friction subjected to various spatial distributions of effective stress (normal stress minus pore pressure). The results indicate that gradual effective stress variations do not provide a sufficient mechanism for the generation of observed seismic response. Strong and abrupt fault heterogeneity, as envisioned in the inherently discrete category, is required for the generation of complex slip patterns and a wide spectrum of event sizes. Strong fault heterogeneity also facilitates the generation of rough rupture fronts capable of radiating high-frequency seismic waves. The large earthquakes in both categories of models occur on a quasi-periodic basis; the degree of periodicity increases with event size and decreases with model complexity. However, in all discrete segmented cases the models generate nonrepeating sequences of earthquakes, and the nature of the large (quasi-periodic) events is highly variable. The results indicate that expectations for regular sequences of earthquakes and/or simple repetitive precursory slip patterns are unrealistic. The frequency-size (FS) statistics of the small failure episodes simulated by the cellular fault models are approximately self-similar with b ≈ 1.2 and bA ≈ 1, where b and bA are b values based on magnitude and rupture area, respectively. For failure episodes larger than a critical size, however, the simulated statistics are strongly enhanced with respect to self-similar distributions defined by the small events. This is due to the fact that the stress concentrated at the edge of a rupture expanding in an elastic solid grows with the rupture size. When the fault properties (e.g., geometric irregularities) are characterized by a narrow range of size scales, the scaling of stress concentrations with the size of the failure zone creates a critical rupture area terminating the self-similar earthquake statistics. In such systems, events reaching the critical size become (on the average) unstoppable, and they continue to grow to a size limited by a characteristic model dimension. When, however, the system is characterized by a broad spectrum of size scales, the above phenomena are suppressed and the range of (apparent) self-similar FS statistics is broad and characterized by average b and bA values of about 1. The simulations indicate that power law extrapolations of low-magnitude seismicity will often underestimate the rate of occurrence of moderate and large earthquakes. The models establish connections between features of FS statistics of earthquakes (range of self-similar regimes, local maxima) and structural properties of faults (dominant size scales of heterogeneities, dimensions of coherent brittle zones). The results suggest that observed FS statistics can be used to obtain information on crustal thickness and fault zone structure.

Earthquake failure sequences along a cellular fault zone in a three‐dimensional elastic solid containing asperity and nonasperity regions

Numerical simulations of earthquake failure sequences along a discrete cellular fault zone are performed for a three-dimensional (3-D) model representing approximately the central San Andreas fault. The model consists of an upper crust overlying a lower crust and mantle region, together defining an elastic half-space with a vertical half-plane fault. The fault contains a region where slip is calculated on a uniform grid of cells governed by a static/kinetic friction law and regions where slip is prescribed so as to represent tectonic loading, aseismic fault creep, and adjacent great earthquakes. The computational region models a 70-km-long and 17.5-km-deep section of the San Andreas fault to the NW of the great 1857 rupture zone. Different distributions of stress drops on failing computational cells are used to model asperity (“Parkfield asperity”) and nonasperity fault regions. The model is “inherently discrete” and corresponds to a situation in which a characteristic size of geometric disorder within the fault (i.e., cell size, here a few hundreds of meters) is much larger than the “nucleation size” (of the order of tens of centimeters to tens of meters) based on slip weakening or state evolution slip distances. The computational grid is loaded by a constant plate motion imposed at the lower crust, upper mantle, and creeping fault regions and by a “staircase” slip history imposed at the 1857 and 1906 rupture zones. Stress transfer along and outside the fault due to the imposed loadings and failure episodes along the computational grid is calculated using 3-D elastic dislocation theory. The resulting displacement field in the computational region is compatible with geodetic and seismological observations only when the asperity and nonasperity regions are characterized by significantly different average stress drops. The frequency-magnitude statistics of the simulated failure episodes are approximately self-similar for small events, with b ≈ 1.2 (the b value of statistics based on rupture area bA is about 1) but are strongly enhanced with respect to self-similarity for events larger than a critical size. This is interpreted as a direct manifestation of our 3-D elastic stress transfer calculations; beyond certain rupture area and potency (seismic moment divided by rigidity) release values, the event is usually unstoppable, and it continues to grow to a size limited by a characteristic model dimension. This effect is not accounted for by cellular automata and block-spring models in which the adopted simplified stress transfer laws fail to scale properly with increasing rupture size. The simulations suggest that local maxima in observed frequency-magnitude statistics correspond to dimensions of coherent brittle zones, such as the width of the seismogenic layer or the length of a fault segment bounded by barriers. The analysis indicates that a single cell size, representing approximately a single scale of geometric disorder, cannot induce self-similarity in a 3-D elastic model over a broad range of magnitudes. A representation of geometric disorder covering a range of scales may thus be required to generate a wide domain of self-similar Gutenberg-Richter statistics. Our simulations show a great diversity in the mode of failure of the Parkfield asperity; the earthquakes themselves define an irregular sequence of events. The modeling, like many other discrete fault models, suggests that expectations for periodic Parkfield earthquakes and/or simple precursory patterns repeating from one event to the other are unrealistic.

Dynamic simulations of slip on a smooth fault in an elastic solid

We report on numerical simulations of slip evolution along a two dimensional (slip varies only with depth) vertical strike-slip fault in an elastic half-space, using a framework incorporating full inertial elastodynamics. The model is a follow-up on earlier quasi-static and quasi-dynamic simulations of deformations along smooth fault systems in elastic continua. The fault is driven below a crustal depth of 24 km by a constant plate velocity of 35-mm/yr. Deformation at each fault location in the crustal zone is the sum of slip rate contributions from rate- and state-dependent friction and power law creep, where both processes have temperature-dependent (and hence depth-varying) coefficients and both take place locally under the same stress. The simulations employ two versions of rate- and state-dependent friction: a “slip” law, which requires nonzero velocity for state evolution, and an ”ageing/slowness“ law, which incorporates state evolution and restrengthening in stationary contact. The assumed constitutive laws and distribution of frictional parameters are compatible with laboratory experiments. The elastodynamic calculations are based on spectral representations of variables and a new algorithm providing a unified computational framework for calculations of long deformational histories containing short periods of rapid instabilities. The simulations show dynamic rupture propagation and wave phenomena not accounted for in the previous quasi-static and quasi-dynamic works. However, the results are qualitatively similar to those obtained by corresponding quasi-static and quasi-dynamic calculations. Slip histories along a smooth fault, simulated here with full elastodynamics for various constitutive laws and model parameters, consist mostly of quasi-periodic large events. This finding indicates that inertial dynamics does not provide a generic mechanism for generating spatio-temporal complexities of slip. On the other hand, calculations done for cases representing, approximately, strongly disordered systems do show rich slip histories with a range of event sizes. This result is compatible with our previous conclusions that the origin of observed broad distributions of earthquake sizes is strong fault zone heterogeneities. The fully dynamic calculations illustrate the evolution of nucleation phases of instabilities associated with accelerating and expanding creep. Final slip values of model earthquakes in full elastodynamic calculations are larger than those of corresponding quasi-static and quasi-dynamic events. The dynamic overshoot in simulations with the slip version of friction is larger than in those employing the ageing/slowness law.

Properties of seismic fault zone waves and their utility for imaging low‐velocity structures

A two-dimensional solution for the scalar wave equation in a model of two vertical layers between two quarter spaces is used to study properties of seismic waves in a laterally heterogeneous low-velocity structure. The waves, referred to as seismic fault zone waves, include head waves, internal fault zone reflections, and trapped waves. The analysis aims to clarify the dependency of the waves on media velocities, media attenuation coefficients, layer widths, and source-receiver geometry. Additional calculations with extreme low-velocity layers provide examples that may be relevant for volcanic and geothermal domains. The interference patterns controlling seismic fault zone waves change with the number of internal reflections in the low-velocity structure. This number increases with propagation distance along the structure, decreases with fault zone width, and increases (for given length scales) with the velocity contrast. The relative lateral position of the source within the low-velocity layer modifies die length scales associated with internal reflections and influences the resulting interference pattern. Low values of Q affect considerably the dominant period and overall duration of the waves. Thus there are significant tradeoffs between propagation distance along the structure, fault zone width, velocity contrast, source location within the fault zone, and Q. The lateral and depth receiver coordinates determine the particular point where the interference pattern is sampled and observed motion is a strong function of both coordinates. The zone connecting sources generating fault zone waves and observation points with appreciable wave amplitude can be over an order of magnitude larger than the fault zone width. Calculations for cases with layer P wave velocity of ∼200 m s−1, modeling a vertical dike or crack with fluid and gas, show conspicuous persisting oscillations. The results resemble aspects of seismic data in volcanic domains. If these waves exist in observed records, their explicit recognition and modeling will help to separate source and structural effects and aid in the interpretation of volcano-seismology signals. Although the tradeoffs in parameters governing seismic fault zone waves are significant, each variable has its own signature, and the parameters may be constrained by additional geophysical data. Simultaneous modeling of many waveforms with an appropriate solution can resolve the various parameters and provide a high-resolution structural image.

Earthquake cycle, fault zones, and seismicity patterns in a rheologically layered lithosphere

We study the coupled evolution of earthquakes and faults in a model consisting of a seismogenic upper crust governed by damage rheology over a viscoelastic substrate. The damage rheology has two types of functional coefficients: (1) a “generalized internal friction” separating states associated with material degradation and healing and (2) damage rate coefficients for positive (degradation) and negative (healing) changes. The evolving damage modifies the effective elastic properties of material in the upper crust as a function of the ongoing deformation. This simulates the creation and healing of fault systems in the upper seismogenic zone. In addition to the vertically averaged thin sheet approximation we introduce a Green function for three-dimensional elastic half-space for the instantaneous component of deformation. The formulation accounts in an internally consistent manner for evolving deformation fields, evolving fault structures, aseismic energy release, and spatiotemporal seismicity patterns. These developments allow us to simulate long histories of crustal deformation and to study the simultaneous evolution of regional earthquakes and faults for various model realizations. To focus on basic features of a large strike-slip fault system, we first consider a simplified geometry of the seismogenic crust by prescribing initial conditions consisting of a narrow damage zone in an otherwise damage-free plate. For this configuration, the model generates an earthquake cycle with distinct interseismic, preseismic, coseismic, and postseismic periods. Model evolution during each period is controlled by a subset of physical properties, which may be constrained by geophysical, geodetic, rock mechanics, and seismological data. In the more generic case with a random initial damage distribution, the model generates large crustal faults and subsidiary branches with complex geometries. The simulated statistics depend on the space-time window of the observational domain. The results indicate that long healing timescale, τh, describing systems with relatively long memory, leads to the development of geometrically regular fault systems and the characteristic frequency-size earthquake distribution. Conversely, short τh (relatively short memory) leads to the development of a network of disordered fault systems and the Gutenberg-Richter earthquake statistics. For intermediate values of τh the results exhibit alternating overall switching of response from periods of intense seismic activity and the characteristic earthquake distribution to periods of low seismic activity and Gutenberg-Richter statistics.