Jeremy E. Kozdon

Earthquake Ruptures with Strongly Rate-Weakening Friction and Off-Fault Plasticity, Part 1: Planar Faults

Abstract Observations demonstrate that faults are fractal surfaces with deviations from planarity at all scales. We study dynamic rupture propagation on self-similar faults having root mean square (rms) height fluctuations of order 10 -3 to 10 -2 times the profile length. Our 2D plane strain models feature strongly rate-weakening fault friction and off-fault Drucker–Prager viscoplasticity. The latter bounds otherwise unreasonably large stress concentrations in the vicinity of bends. Our choice of a cohesionless yield function prevents tensile stress states and thus fault opening. A consequence of strongly rate-weakening friction is the existence of a critical background stress level above which self-sustaining rupture propagation, in the form of self-healing slip pulses, first becomes possible. Around this level, at which natural faults are expected to operate, ruptures become extremely sensitive to fault roughness and exhibit substantial fluctuations in rupture velocity. Except for shallow inclinations of the maximum compressive stress to the fault (less than about 20°), the fluctuations are anticorrelated with the local fault slope. These accelerations and decelerations of the rupture, together with naturally emerging slip heterogeneity, excite waves of all wavelengths and result in ground acceleration spectra that are flat at high frequency, consistent with observed strong motion records.

Simulation of Dynamic Earthquake Ruptures in Complex Geometries Using High-Order Finite Difference Methods

We develop a stable and high-order accurate finite difference method for problems in earthquake rupture dynamics in complex geometries with multiple faults. The bulk material is an isotropic elastic solid cut by pre-existing fault interfaces that accommodate relative motion of the material on the two sides. The fields across the interfaces are related through friction laws which depend on the sliding velocity, tractions acting on the interface, and state variables which evolve according to ordinary differential equations involving local fields.The method is based on summation-by-parts finite difference operators with irregular geometries handled through coordinate transforms and multi-block meshes. Boundary conditions as well as block interface conditions (whether frictional or otherwise) are enforced weakly through the simultaneous approximation term method, resulting in a provably stable discretization.The theoretical accuracy and stability results are confirmed with the method of manufactured solutions. The practical benefits of the new methodology are illustrated in a simulation of a subduction zone megathrust earthquake, a challenging application problem involving complex free-surface topography, nonplanar faults, and varying material properties.

Interaction of Waves with Frictional Interfaces Using Summation-by-Parts Difference Operators: Weak Enforcement of Nonlinear Boundary Conditions

We present a high-order difference method for problems in elastodynamics involving the interaction of waves with highly nonlinear frictional interfaces. We restrict our attention to two-dimensional antiplane problems involving deformation in only one direction. Jump conditions that relate tractions on the interface, or fault, to the relative sliding velocity across it are of a form closely related to those used in earthquake rupture models and other frictional sliding problems. By using summation-by-parts (SBP) finite difference operators and weak enforcement of boundary and interface conditions, a strictly stable method is developed. Furthermore, it is shown that unless the nonlinear interface conditions are formulated in terms of characteristic variables, as opposed to the physical variables in terms of which they are more naturally stated, the semi-discretized system of equations can become extremely stiff, preventing efficient solution using explicit time integrators.The use of SBP operators also provides a rigorously defined energy balance for the discretized problem that, as the mesh is refined, approaches the exact energy balance in the continuous problem. This enables one to investigate earthquake energetics, for example the efficiency with which elastic strain energy released during rupture is converted to radiated energy carried by seismic waves, rather than dissipated by frictional sliding of the fault. These theoretical results are confirmed by several numerical tests in both one and two dimensions demonstrating the computational efficiency, the high-order convergence rate of the method, the benefits of using strictly stable numerical methods for long time integration, and the accuracy of the energy balance.

Choosing weight functions in iterative methods for simple roots

Weight functions with a parameter are introduced into an iteration process to increase the order of the convergence and enhance the behavior of the iteration process. The parameter can be chosen to restrict extraneous fixed points to the imaginary axis and provide the best basin of attraction. The process is demonstrated on several examples.

Stable coupling of nonconforming, high-order finite difference methods

A methodology for handling block-to-block coupling of nonconforming, multiblock summation-by-parts finite difference methods is proposed. The coupling is based on the construction of projection operators that move a finite difference grid solution along an interface to a space of piecewise defined functions; we specifically consider discontinuous, piecewise polynomial functions. The constructed projection operators are compatible with the underlying summation-by-parts energy norm. Using the linear wave equation in two dimensions as a model problem, energy stability of the coupled numerical method is proven for the case of curved, nonconforming block-to-block interfaces. To further demonstrate the power of the coupling procedure, we show how it allows for the development of a provably energy stable coupling between curvilinear finite difference methods and a curved-triangle discontinuous Galerkin method. The theoretical results are verified through numerical simulations on curved meshes as well as eigenval...

Multidimensional upstream weighting for multiphase transport on general grids

The governing equations for multiphase flow in porous media have a mixed character, with both nearly elliptic and nearly hyperbolic variables. The flux for each phase can be decomposed into two parts: (1) a geometry- and rock-dependent term that resembles a single-phase flux; and (2) a mobility term representing fluid properties and rock–fluid interactions. The first term is commonly discretized by two- or multipoint flux approximations (TPFA and MPFA, respectively). The mobility is usually treated with single-point upstream weighting (SPU), also known as dimensional or donor cell upstream weighting. It is well known that when simulating processes with adverse mobility ratios, SPU suffers from grid orientation effects. An important example of this, which will be considered in this work, is the displacement of a heavy oil by water. For these adverse mobility ratio flows, the governing equations are unstable at the modeling scale, rendering a challenging numerical problem. These challenges must be addressed in order to avoid systematic biasing of simulation results. In this work, we present a framework for multidimensional upstream weighting for multiphase flow with buoyancy on general two-dimensional grids. The methodology is based on a dual grid, and the resulting transport methods are provably monotone. The multidimensional transport methods are coupled with MPFA methods to solve the pressure equation. Both explicit and fully implicit approaches are considered for time integration of the transport equations. The results show considerable reduction of grid orientation effects compared to SPU, and the explicit multidimensional approach allows larger time steps. For the implicit method, the total number of non-linear iterations is also reduced when multidimensional upstream weighting is used.

EARTHQUAKE RUPTURES ON ROUGH FAULTS

Natural fault surfaces exhibit roughness at all scales, with root-mean-square height fluctuations of order 10− 3 to 10− 2 times the profile length. We study earthquake rupture propagation on such faults, using strongly rate-weakening fault friction and off-fault plasticity. Inelastic deformation bounds stresses to reasonable values and prevents fault opening. Stress perturbations induced by slip on rough faults cause irregular rupture propagation and the production of incoherent high-frequency ground motion.

A Suite of Exercises for Verifying Dynamic Earthquake Rupture Codes

We describe a set of benchmark exercises that are designed to test if computer codes that simulate dynamic earthquake rupture are working as intended. These types of computer codes are often used to understand how earthquakes operate, and they produce simulation results that include earthquake size, amounts of fault slip, and the patterns of ground shaking and crustal deformation. The benchmark exercises examine a range of features that scientists incorporate in their dynamic earthquake rupture simulations. These include implementations of simple or complex fault geometry, off‐fault rock response to an earthquake, stress conditions, and a variety of formulations for fault friction. Many of the benchmarks were designed to investigate scientific problems at the forefronts of earthquake physics and strong ground motions research. The exercises are freely available on our website for use by the scientific community.

Monotone Multi-dimensional Upstream Weighting on General Grids

The governing equations for multi-phase flow in porous media often have a mixed elliptic and (nearly) hyperbolic character. The total flux for each phase consists of two parts: a geometry and rock dependent term that resembles a single-phase flux and a mobility term representing fluid properties and rock-fluid interactions. The geometric term is commonly discretized by two or multi point flux approximations (TPFA and MPFA, respectively). The mobility is usually treated with single point upstream weighting (SPU), also known as dimensional or donor cell upstream weighting. It is well known that when simulating processes with adverse mobility ratios, e.g. gas injection, SPU yields grid orientation effects. For these physical processes, the governing equations are unstable on the scale at which they are discretized, rendering a challenging numerical problem. These challenges must be addressed in order to avoid systematic biasing of simulation results and to improve the overall performance prediction of enhanced oil recovery processes. In this work, we present a framework for multi-dimensional upstream weighting for multi-phase flow with gravity on general two-dimensional grids. The methodology is based on a dual grid, and the resulting transport methods are provably monotone. The multi-dimensional transport methods are coupled with MPFA methods to solve the pressure equation. Both explicit and fully implicit approaches are considered for treatment of the transport equations. The results show considerable reduction of grid orientation effects compared to SPU, and the explicit multi-dimensional approach allows larger time steps. For the implicit method, the total number of non-linear iterations is also reduced when multi-dimensional upstream weighting is used.

Boundary conditions and stability of a perfectly matched layer for the elastic wave equation in first order form

In computations, it is now common to surround artificial boundaries of a computational domain with a perfectly matched layer (PML) of finite thickness in order to prevent artificially reflected waves from contaminating a numerical simulation. Unfortunately, the PML does not give us an indication about appropriate boundary conditions needed to close the edges of the PML, or how those boundary conditions should be enforced in a numerical setting. Terminating the PML with an inappropriate boundary condition or an unstable numerical boundary procedure can lead to exponential growth in the PML which will eventually destroy the accuracy of a numerical simulation everywhere. In this paper, we analyze the stability and the well-posedness of boundary conditions terminating the PML for the elastic wave equation in first order form. First, we consider a vertical modal PML truncating a two space dimensional computational domain in the horizontal direction. We freeze all coefficients and consider a left half-plane problem with linear boundary conditions terminating the PML. The normal mode analysis is used to study the stability and well-posedness of the resulting initial boundary value problem (IBVP). The result is that any linear well-posed boundary condition yielding an energy estimate for the elastic wave equation, without the PML, will also lead to a well-posed IBVP for the PML. Second, we extend the analysis to the PML corner region where both a horizontal and vertical PML are simultaneously active. The challenge lies in constructing accurate and stable numerical approximations for the PML and the boundary conditions. Third, we develop a high order accurate finite difference approximation of the PML subject to the boundary conditions. To enable accurate and stable numerical boundary treatments for the PML we construct continuous energy estimates in the Laplace space for a one space dimensional problem and two space dimensional PML corner problem. We use summation-by-parts finite difference operators to approximate the spatial derivatives and impose boundary conditions weakly using penalties. In order to ensure numerical stability of the discrete PML, it is necessary to extend the numerical boundary procedure to the auxiliary differential equations. This is crucial for deriving discrete energy estimates analogous to the continuous energy estimates. Numerical experiments are presented corroborating the theoretical results. Moreover, in order to ensure longtime numerical stability, the boundary condition closing the PML, or its corresponding discrete implementation, must be dissipative. Furthermore, the numerical experiments demonstrate the stable and robust treatment of PML corners.