Andrew M. Bradley

The Role of Thermal Pressurization and Dilatancy in Controlling the Rate of Fault Slip

Geophysical observations have shown that transient slow slip events, with average slip speeds v on the order of 10−8 to 10−7 m/s, occur in some subduction zones. These slip events occur on the same faults but at greater depth than large earthquakes (with slip speeds of order ∼ 1 m/s). We explore the hypothesis that whether slip is slow or fast depends on the competition between dilatancy, which decreases fault zone pore pressure p, and thermal pressurization, which increases p. Shear resistance to slip is assumed to follow an effective stress law τ = f (σ− p) ≡ f σ. We present two-dimensional quasi-dynamic simulations that include rate-state friction, dilatancy, and heat and pore fluid flow normal to the fault. We find that at lower background effective normal stress (σ), slow slip events occur spontaneously, whereas at higher σ, slip is inertially limited. At intermediate σ, dynamic events are followed by quiescent periods, and then long durations of repeating slow slip events. In these cases, accelerating slow events ultimately nucleate dynamic rupture. Zero-width shear zone approximations are adequate for slow slip events, but substantially overestimate the pore pressure and temperature changes during fast slip when dilatancy is included.

Time-dependent dike propagation from joint inversion of seismicity and deformation data

Dike intrusions both deform Earths surface and induce propagating earthquake swarms. We develop methods to utilize both deformation and seismicity from brittle, volcano-tectonic earthquakes to image time-dependent dike propagation. Dieterichs (1994) seismicity-rate theory is used to relate dike-induced stress changes to seismicity rate and is combined with elastic Greens functions relating dike opening to deformation. Different space-time patterns of seismicity develop if earthquakes occur at the same depth as the dike compared to above/below the dike. In the former, seismicity initiates near the dikes leading edges but shuts off as the dike tips pass and seismogenic volumes fall into stress shadows. In the latter, seismicity continues at a decaying rate after the tips pass. We focus on lateral propagation and develop a nonlinear inversion method that estimates dike length and pressure as a function of time. The method is applied to the 2007 Fathers Day intrusion in Kilauea Volcano. Seismicity is concentrated at ∼3 km depth, comparable to geodetic estimates of dike depth, and decays rapidly in time. With lateral propagation only and a vertical dike-tip line, it is difficult to fit both GPS data and the rapid down-rift jump in seismicity, suggesting significant vertical propagation. For the events to have occurred below the dike requires a very short aftershock decay time, hence unreasonably high background stressing rate. The rapid decay is better explained if the dike extends somewhat below the seismicity. We suggest that joint inversion is useful for studying the diking process and may allow for improved short-term eruption forecasts.

Software for Efficient Static Dislocation–Traction Calculations in Fault Simulators

Online Material: Software packages with documentation. Quasistatic or quasidynamic rate–state friction (QRSF) simulators are used to study the mechanics of faults (e.g., Shibazaki and Shimamoto, 2007). The displacement discontinuity method (DDM; Crouch and Starfield, 1983) meshes the fault or faults into N elements and constructs a matrix of Green’s functions (GFs) relating slip to stress. The simulator evolves strength and slip in time. Usually, the most expensive part of a simulation time step is the matrix–vector product (MVP) of the slip distribution with the DDM matrix; the straightforward implementation performs O ( N 2) operations. A simulator performs thousands to millions of MVPs with the same DDM matrix. My free and open‐source software (FOSS) package hmmvp speeds up simulations by an asymptotic factor of a little less than N faster than the straightforward implementation for both forming an approximation to and performing MVP with the GF matrix. In QRSF simulations, rupture tip length scales as f  r ∝ μ  ′ d c /( bσ ) for the aging evolution law and f  r multiplied by a factor that depends on slip speed and background values for the slip law. In this case, μ  ′= μ /(1− ν ), μ is the shear modulus, ν is Poisson’s ratio, b is the constant multiplying the state term in rate–state friction, σ is the effective normal stress, and d c is the characteristic slip distance for friction evolution (Rubin and Ampuero, 2009). Rupture tips must be well resolved in simulations. If f  r is nonuniform on the fault, discretization also can be nonuniform for efficiency. However, not all DDM operators are accurate on nonuniform meshes. My FOSS package dc3dm implements a method, IGA, that is as accurate on a nonuniform mesh as the standard method is on a uniform mesh. dc3dm uses hmmvp to efficiently approximate …

Designing vector-friendly compact BLAS and LAPACK kernels

Many applications, such as PDE based simulations and machine learning, apply blas/lapack routines to large groups of small matrices. While existing batched blas APIs provide meaningful speedup for this problem type, a non-canonical data layout enabling cross-matrix vectorization may provide further significant speedup. In this paper, we propose a new compact data layout that interleaves matrices in blocks according to the SIMD vector length. We combine this compact data layout with a new interface to blas/lapack routines that can be used within a hierarchical parallel application. Our layout provides up to 14X, 45X, and 27X speedup against OpenMP loops around optimized dgemm, dtrsm and dgetrf kernels, respectively, on the Intel Knights Landing architecture. We discuss the compact batched blas/lapack implementations in two libraries, KokkosKernels and Intel® Math Kernel Library. We demonstrate the APIs in a line solver for coupled PDEs. Finally, we present detailed performance analysis of our kernels.

Bounding the moment deficit rate on crustal faults using geodetic data: Methods

The geodetically-derived interseismic Moment Deficit Rate (MDR) provides a first-order constraint on earthquake potential, and can play an important role in seismic hazard assessment, but quantifying uncertainty in MDR is a challenging problem that has not been fully addressed. We establish criteria for reliable MDR estimators, evaluate existing methods for determining the probability density of MDR, and propose and evaluate new methods. Geodetic measurements moderately far from the fault provide tighter constraints on MDR than those nearby. Previously used methods can fail catastrophically under predictable circumstances. The bootstrap method works well with strong data constraints on MDR, but can be strongly biased when network geometry is poor. We propose two new methods: the Constrained Optimization Bounding Estimator (COBE) assumes uniform priors on slip rate (from geologic information) and MDR, and can be shown through synthetic tests to be a useful, albeit conservative estimator; the Constrained Optimization Bounding Linear Estimator (COBLE) is the corresponding linear estimator with Gaussian priors rather than point-wise bounds on slip rates. COBE matches COBLE with strong data constraints on MDR. We compare results from COBE and COBLE to previously published results for the interseismic MDR at Parkfield, on the San Andreas Fault, and find similar results; thus the apparent discrepancy between MDR and the total moment release (seismic and afterslip) in the 2004 Parkfield earthquake remains.

Constraining the Magmatic System at Mount St. Helens (2004–2008) Using Bayesian Inversion With Physics-Based Models Including Gas Escape and Crystallization

Physics-based models of volcanic eruptions track conduit processes as functions of depth and time. When used in inversions, these models permit integration of diverse geological and geophysical data sets to constrain important parameters of magmatic systems. We develop a 1-D steady state conduit model for effusive eruptions including equilibrium crystallization and gas transport through the conduit and compare with the quasi-steady dome growth phase of Mount St. Helens in 2005. Viscosity increase resulting from pressure-dependent crystallization leads to a natural transition from viscous flow to frictional sliding on the conduit margin. Erupted mass flux depends strongly on wall rock and magma permeabilities due to their impact on magma density. Including both lateral and vertical gas transport reveals competing effects that produce nonmonotonic behavior in the mass flux when increasing magma permeability. Using this physics-based model in a Bayesian inversion, we link data sets from Mount St. Helens such as extrusion flux and earthquake depths with petrological data to estimate unknown model parameters, including magma chamber pressure and water content, magma permeability constants, conduit radius, and friction along the conduit walls. Even with this relatively simple model and limited data, we obtain improved constraints on important model parameters. We find that the magma chamber had low (<5 wt %) total volatiles and that the magma permeability scale is well constrained at ∼10−11.4m2 to reproduce observed dome rock porosities. Compared with previous results, higher magma overpressure and lower wall friction are required to compensate for increased viscous resistance while keeping extrusion rate at the observed value.