Geophysical and geological studies provide evidence for cyclic changes in fault-zone pore fluid pressure that synchronize with or at least modulate seismic cycles. A hypothesized mechanism for this behavior is fault valving arising from temporal changes in fault zone permeability. In our study, we investigate the coupled dynamics of rate and state friction, along-fault fluid flow, and permeability evolution. Permeability decreases with time, and increases with slip. Linear stability analysis shows that steady slip with constant fluid flow along the fault zone is unstable to perturbations, even for velocity-strengthening friction with no state evolution, if the background flow is sufficiently high. We refer to this instability as the “fault valve instability.’ The propagation speed of the fluid pressure and slip pulse can be much higher than expected from linear pressure diffusion, and it scales with permeability enhancement. Two-dimensional simulations with spatially uniform properties show that the fault valve instability develops into slow slip events, in the form of aseismic slip pulses that propagate in the direction of fluid flow. We also perform earthquake sequence simulations on a megathrust fault, taking into account depth-dependent frictional and hydrological properties. The simulations produce quasi-periodic slow slip events from the fault valve instability below the seismogenic zone, in both velocity-weakening and velocity-strengthening regions, for a wide range of effective normal stresses. A separation of slow slip events from the seismogenic zone, which is observed in some subduction zones, is reproduced when assuming a fluid sink around the mantle wedge corner.

Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have rapidly progressed to address fundamental problems in fault mechanics and provide self-consistent, physics-based frameworks to interpret and predict geophysical observations across spatial and temporal scales. To advance SEAS simulations with rigor and reproducibility, we pursue community efforts to verify numerical codes in an expanding suite of benchmarks. Here we present code comparison results from a new set of benchmark problems BP6-QD-A/S/C that consider a single aseismic slip transient induced by changes in pore fluid pressure consistent with fluid injection and diffusion in fault models with different treatments of fault friction. Ten modeling groups participated in problems BP6-QD-A and BP6-QD-S considering rate-and-state fault models using the aging and slip law formulations for frictional state evolution, respectively, allowing us to explore these ingredients across multiple codes and better understand how various computational factors affect the simulated evolution of pore pressure and aseismic slip. Comparisons of problems using the aging versus slip law illustrate how models of aseismic slip can differ in the timing and amount of slip achieved with different treatments of fault friction given the same perturbations in pore fluid pressure. We achieve excellent quantitative agreement across participating codes, with further agreement being found by ensuring sufficiently fine time-stepping and consistent treatment of remote boundary conditions. Our benchmark efforts offer a community-based example to reveal sensitivities of numerical modeling results, which is essential for advancing multi-physics SEAS models to better understand and construct reliable predictive models of fault dynamics.

A document by Katherine Coppess. Click on the document to view its contents.

A document by Taiyi Wang. Click on the document to view its contents.

It is widely recognized that fluid injection can trigger fault slip. However, the processes by which the fluid-rock interactions facilitate or inhibit slip are poorly understood and some are neglected or oversimplified in most models of injection-induced slip. In this study, we perform a 2D antiplane shear investigation of aseismic slip that occurs in response to fluid injection into a permeable fault governed by rate-and-state friction. We account for pore dilatancy and permeability changes that accompany slip, and quantify how these processes affect pore pressure diffusion, which couples to aseismic slip. The fault response to injection has two phases. In the first phase, slip is negligible and pore pressure closely follows the standard linear diffusion model. Pressurization of the fault eventually triggers aseismic slip in the immediate vicinity of the injection site. In the second phase, the aseismic slip front expands outward and dilatancy causes pore pressure to depart from the linear diffusion model. Aseismic slip front overtakes pore pressure contours, with both subsequently advancing at constant rate along fault. We quantify how prestress, initial state variable, injection rate, and frictional properties affect the migration rate of the aseismic slip front, finding values ranging from less than 50 to 1000 m/day for typical parameters. Additionally, we compare to the case when porosity and permeability evolution are neglected. In this case, the aseismic slip front migration rate and total slip are much higher. Our modeling demonstrates that porosity and permeability evolution, especially dilatancy, fundamentally alters how faults respond to fluid injection.

Recent research in real-time tsunami early warning can be broadly classified into two approaches. The first involves the use of seismic and regional geodetic data to calculate the tsunami wavefield indirectly through the estimation of earthquake source parameters. The second directly reconstructs the tsunami wavefield using data assimilation of ocean-bottom pressure sensor data such as those from DONET and S-NET (Maeda et al. 2015, Gusman et al. 2016). Data assimilation interpolates between the numerical solution and the observations to make the forecast more consistent with real data. Currently, the most popular method for forecasting the waveform is optimal interpolation, which uses a Kalman filter (KF) like approach, but holds the Kalman gain matrix fixed to reduce the runtime. This approach, coupled with tsunami Green’s functions, is very efficient and generates useful predictions. Here, we demonstrate that more accurate and stable forecasts can be obtained using the ensemble KF (enKF), a more computationally efficient variant of KF, in which the gain matrix is updated according to the physical model and the evolution of the error covariance matrix. The ensemble representation is a form of dimensionality reduction, in that only a small ensemble is propagated, instead of the joint distribution including the full covariance matrix. This method also provides a means to obtain the probability distribution of the forecast at each grid point location. We use a scenario tsunami in the Cascadia subduction zone, generated from a 2D fully-coupled dynamic rupture simulation (Lotto et al., submitted 2018). Randomly perturbed tsunami wave height data is used in the assimilation process, as we propagate the wave using a 1D linear shallow water code on a staggered grid. Better waveform agreement is achieved even in the early stages of assimilation, with much less fluctuation compared to optimal interpolation. We also explore spatial and temporal aliasing effects, in terms of the relation between observation station spacing and wavelength, as well as between assimilation and forecast time intervals. Although enKF is computationally more expensive, we are working on a fast, parallelized GPU implementation, which will significantly reduce the runtime, taking us a step closer to reliable real-time tsunami early warning.

Fluids influence fault zone strength and the occurrence of earthquakes, slow slip events, and aseismic slip. We introduce an earthquake sequence model with fault zone fluid transport, accounting for elastic, viscous, and plastic porosity evolution, with permeability having a power-law dependence on porosity. Fluids, sourced at a constant rate below the seismogenic zone, ascend along the fault. While the modeling is done for a vertical strike-slip fault with 2D antiplane shear deformation, the general behavior and processes are anticipated to apply also to subduction zones. The model produces large earthquakes in the seismogenic zone, whose recurrence interval is controlled in part by compaction-driven pressurization and weakening. The model also produces a complex sequence of slow slip events (SSEs) beneath the seismogenic zone. The SSEs are initiated by compaction-driven pressurization and weakening and stalled by dilatant suctions. Modeled SSE sequences include long-term events lasting from a few months to years and very rapid short-term events lasting for only a few days; slip is ~1-10 cm. Despite ~1-10 MPa pore pressure changes, porosity and permeability changes are small and hence fluid flux is relatively constant except in the immediate vicinity of slip fronts. This contrasts with alternative fault valving models that feature much larger changes in permeability from the evolution of pore connectivity. Our model demonstrates the important role that compaction and dilatancy have on fluid pressure and fault slip, with possible relevance to slow slip events in subduction zones and elsewhere.

From interpreting data to scenario modeling of subduction events, numerical modeling has been crucial for studying tsunami generation by earthquakes. Seafloor instruments in the source region feature complex signals containing a superposition of seismic, ocean acoustic, and tsunami waves. Rigorous modeling is required to interpret these data and use them for tsunami early warning. However, previous studies utilize separate earthquake and tsunami models, with one-way coupling between them and approximations that might limit the applicability of the modeling technique. In this study, we compare four earthquake-tsunami modeling techniques, highlighting assumptions that affect the results, and discuss which techniques are appropriate for various applications. Most techniques couple a 3D Earth model with a 2D depth-averaged shallow water tsunami model. Assuming the ocean is incompressible and that tsunami propagation is negligible over the earthquake duration leads to technique (1), which equates earthquake seafloor uplift to initial tsunami sea surface height. For longer duration earthquakes, it is appropriate to follow technique (2), which uses time-dependent earthquake seafloor velocity as a time-dependent forcing in the tsunami mass balance. Neither technique captures ocean acoustic waves, motivating newer techniques that capture the seismic and ocean acoustic response as well as tsunamis. Saito et al. (2019) propose technique (3), which solves the 3D elastic and acoustic equations to model the earthquake rupture, seismic wavefield, and response of a compressible ocean without gravity. Then, sea surface height is used as a forcing term in a tsunami simulation. A superposition of the earthquake and tsunami solutions provides the complete wavefield, with one-way coupling. The complete wavefield is also captured in technique (4), which utilizes a fully-coupled solid Earth and ocean model with gravity (Lotto & Dunham, 2015). This technique, recently incorporated into the 3D code SeisSol, simultaneously solves earthquake rupture, seismic waves, and ocean response (including gravity). Furthermore, we show how technique (3) follows from (4) subject to well-justified approximations.

Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent, physics-based framework to connect, interpret, and predict diverse geophysical observations across spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare and verify numerical codes based on boundary-element, finite-element, and finite-difference methods, in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space, where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs for sufficiently large model domains and small grid spacings. However, discrepancies in rupture fronts of the initial event are influenced by the free surface and various computational factors. The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to numerical resolution and domain-size-dependent loading. Despite such variability, key properties of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer a community-based example to improve numerical simulations and reveal sensitivities of model observables, which are important for advancing SEAS models to better understand earthquake system dynamics.

Infrasound observations are increasingly used to constrain properties of volcanic eruptions. In order to better interpret infrasound observations, however, there is a need to better understand the relationship between eruption properties and sound generation. Here we perform two-dimensional computational aeroacoustic simulations where we solve the compressible Navier-Stokes equations for pure-air with a large-eddy simulation approximation. We simulate idealized impulsive volcanic eruptions where the exit velocity is specified and the eruption is pressure-balanced with the atmosphere. Our nonlinear simulation results are compared with the commonly-used analytical linear acoustics model of a compact monopole source radiating acoustic waves isotropically in a half space. The monopole source model matches the simulations for low exit velocities (<100 m/s or M ~ 0.3 where M is the Mach number); however, the two solutions diverge as the exit velocity increases with the simulations developing lower peak amplitude, more rapid onset, and anisotropic radiation with stronger infrasound signals recorded above the vent than on Earth’s surface. Our simulations show that interpreting ground-based infrasound observations with the monopole source model can result in an underestimation of the erupted volume for eruptions with sonic or supersonic exit velocities. We examine nonlinear effects and show that nonlinear effects during propagation are relatively minor for the parameters considered. Instead, the dominant nonlinear effect is advection by the complex flow structure that develops above the vent. This work demonstrates the need to consider anisotropic radiation patterns and jet dynamics when interpreting infrasound observations, particularly for eruptions with sonic or supersonic exit velocities.

Recent inversions of InSAR ground surface deformation in the Delaware Basin have revealed aseismic slip on conjugate, semi-optimally oriented normal faults. Slip, occurring over 3-5 years, extends approximately 1 km down-dip, over 10 km along strike, and reaches 23 cm. Several disposal wells are located nearby, suggesting that pressure changes from injection might have induced slip. We develop and calibrate 2D and pseudo-3D coupled pore pressure diffusion and rate-state models with velocity-strengthening friction. Pore pressure diffusion is limited to a high-permeability fault damage zone, and the net influx of fluid is adjusted to match the observed slip. A pressure increase of 1-2 MPa is required to initiate slip. Most aseismic slip occurs at approximately constant friction and ~5 MPa additional pressure increase is needed to reach ~20 cm slip. Our work suggests that models of the type developed here can be used to operationally manage injection-induced aseismic slip.

Explosive volcanic eruptions radiate seismic waves as a consequence of pressure and shear traction changes within the conduit/chamber system. Kinematic source inversions utilize these waves to determine equivalent seismic force and moment tensor sources, but relation to eruptive processes is often ambiguous and nonunique. In this work, we provide an alternative, forward modeling approach to calculate moment tensor and force equivalents of a model of eruptive conduit flow and chamber depressurization. We explain the equivalence of two seismic force descriptions, the first in terms of traction changes on the conduit/chamber walls, and the second in terms of changes in magma momentum, weight, and momentum transfer to the atmosphere. Eruption onset is marked by a downward seismic force, associated with loss of restraining shear tractions from fragmentation. This is followed by a much larger upward seismic force from upward drag of ascending magma and reduction of magma weight remaining in the conduit/chamber system. The static force is upward, arising from weight reduction. We calculate synthetic seismograms to examine the expression of eruptive processes at different receiver distances. Filtering these synthetics to the frequency band typically resolved by broadband seismometers produces waveforms similar to very long period (VLP) seismic events observed in strombolian and vulcanian eruptions. However, filtering heavily distorts waveforms, accentuating processes in early, unsteady parts of eruptions and eliminating information about longer time scale depressurization and weight changes that dominate unfiltered seismograms. The workflow we have introduced can be utilized to directly and quantitatively connect eruption models with seismic observations.