Search

Geophysical and geological studies provide evidence for cyclic changes in fault‐zone pore fluid pressure that synchronize with or at least modulate slip events. A hypothesized explanation is fault valving arising from temporal changes in fault zone permeability. In our study, we investigate how the coupled dynamics of rate and state friction, along‐fault fluid flow, and permeability evolution can produce slow slip events. 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, which scales with permeability enhancement, can be much higher than expected from linear pressure diffusion. 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. Plain Language Summary: Slow slip events are observed in subduction zones worldwide. Their mechanism is not well understood, but geophysical and geological research suggests a relation with recurring changes in fluid pressure within the fault zone. Here we explore the fault valve mechanism for slow slip events using mathematical and computational models that couple fluid flow through fault zones with frictional slip on faults. The fault valve mechanism (arising from cyclic changes in the permeability or resistance to fluid flow) produces pulses of high fluid pressure, accompanied by slow slip, that advance along the fault in the direction of fluid flow. We quantify the conditions under which this occurs as well as observable properties like the propagation speed and rate of occurrence of slow slip events. We also perform simulations of subduction zone slow slip events using fault zone and frictional properties that vary with depth in a realistic manner. The simulations show that the fault valve mechanism can produce slow slip events with approximately the observed rate of occurrence, while also highlighting some discrepancies with observations that must be addressed in future work. Key Points: We analyze the dynamics of fault slip with fault‐zone fluid flow and fault‐parallel permeability enhancement with slip and sealing with timeFault‐valve instability produces unidirectional aseismic slip and pore pressure pulses even with velocity‐strengthening frictionSubduction zone earthquake cycle simulations show that the fault‐valve instability can produce slow slip events below the seismogenic zone [ABSTRACT FROM AUTHOR]

All instrumented basaltic caldera collapses have generated Mw > 5 very long period earthquakes. However, previous studies of source dynamics have been limited to lumped models treating the caldera block as rigid, leaving open questions related to how ruptures initiate and propagate around the ring fault, and the seismic expressions of those dynamics. We present the first 3D numerical model capturing the nucleation and propagation of ring fault rupture, the mechanical coupling to the underlying viscoelastic magma, and the associated seismic wavefield. We demonstrate that seismic radiation, neglected in previous models, acts as a damping mechanism reducing coseismic slip by up to half, with effects most pronounced for large magma chamber volume/ring fault radius or highly compliant crust/compressible magma. Viscosity of basaltic magma has negligible effect on collapse dynamics. In contrast, viscosity of silicic magma significantly reduces ring fault slip. We use the model to simulate the 2018 Kı̄lauea caldera collapse. Three stages of collapse, characterized by ring fault rupture initiation and propagation, deceleration of the downward‐moving caldera block and magma column, and post‐collapse resonant oscillations, in addition to chamber pressurization, are identified in simulated and observed (unfiltered) near‐field seismograms. A detailed comparison of simulated and observed displacement waveforms corresponding to collapse earthquakes with hypocenters at various azimuths of the ring fault reveals a complex nucleation phase for earthquakes initiated on the northwest. Our numerical simulation framework will enhance future efforts to reconcile seismic and geodetic observations of caldera collapse with conceptual models of ring fault and magma chamber dynamics. Plain Language Summary: Caldera collapse manifests as the rapid subsidence of a kilometer‐scale block of crust circumscribed by a near‐circular fault on top of a volcano. The subsidence of the caldera block is caused by the eruption‐induced withdrawal of magma and reduction in pressure in the underlying magma chamber. All scientifically instrumented caldera collapses at volcanoes with low‐viscosity magma are accompanied by earthquakes of magnitude 5 and above. How do magma viscosity and the seismic wave radiation influence the amount of slip per earthquake on the fault? What can we learn about the dynamics of these earthquakes from seismic records? We address these questions by performing computer simulations of caldera collapse earthquakes and compare the results to the seismic records from the Kı̄lauea caldera collapse of 2018. Key Points: Seismic wave radiation through ring fault and magma, as well as high magma viscosity, reduce fault slip by up to half during collapseRupture propagation, downward momentum transfer via magma pressure waves, and chamber pressurization are identified in unfiltered seismogramsComparison between simulated and observed near‐field seismograms from Kı̄lauea 2018 reveals complex nucleation phase on the ring fault [ABSTRACT FROM AUTHOR]

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. Plain Language Summary: Water in the crust plays an important role in controlling the strength of fault zones and frictional sliding, which manifest as earthquakes and slow slip events that do not produce ground shaking. In this study, we perform computer modeling of earthquake sequences that are coupled to the evolution of fluid pressure and rock properties. In particular, compaction or dilation of the water‐filled pore space in rock drives changes in fluid pressure and influences the fault's frictional resistance to slip. The model quantifies the effects of compaction and dilation on both large earthquakes and slow slip events, providing specific predictions regarding slow slip event properties, pressure changes, and changes in fluid flow that might be testable with geophysical, geologic, and geochemical data. Key Points: Fluid pressure changes from pore compaction and dilatancy influence slow slip events in our modelModeled slow slip events span long‐term events lasting for months to years to rapid short‐term events for a few daysEarthquake recurrence interval is controlled in part by compaction‐driven pressurization and weakening [ABSTRACT FROM AUTHOR]

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 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 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 (ultra long period time scale depressurization and weight changes that dominate unfiltered seismograms. Our workflow can be utilized to directly and quantitatively connect eruption models with seismic observations. Plain Language Summary: Volcanic eruptions radiate seismic waves that can be recorded by seismometers placed on and around a volcano. Analysis of seismic data enables one to study eruptions, in particular the processes occurring in the magma‐filled conduit and chamber that feeds the eruption. One process of particular interest is fragmentation, in which magma containing a mixture of liquid melt and gas bubbles breaks apart in the conduit and erupts explosively from the vent. We perform computer simulations of explosive eruptions and then use the output of those simulations to predict seismic radiation. We examine the seismograms produced by this workflow to identify features that are diagnostic of process, such as fragmentation, that occur at different times in the eruption. These predictions will guide interpretation of seismic data from real eruptions. Key Points: We provide expressions for the seismic force in eruptions that can be evaluated using outputs from unsteady conduit flow modelsWe generate and analyze seismograms from vulcanian eruption models to identify the seismic expression of fragmentation and other processesOur modeling suggests that eruptive mass could be inferred by extending seismic inversions to periods of minutes to tens of minutes [ABSTRACT FROM AUTHOR]

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. Plain Language Summary: Earthquakes and fault zone processes occur over time scales ranging from milliseconds to millennia and longer. Computational models are increasingly used to simulate sequences of earthquakes and aseismic slip (SEAS). These simulations can be connected to diverse geophysical observations, offering insights into earthquake system dynamics. To improve these simulations, we pursue community efforts to design benchmarks for 3D SEAS problems. We involve earthquake researchers around the globe to compare simulation results using different numerical codes. We identify major factors that contribute to the discrepancies among simulations. For example, the spatial dimension and resolution of the computational model can affect how earthquakes start and grow, as well as how frequently they recur. Code comparisons are more challenging when we consider the Earth's surface in the simulations. Fortunately, we find that several key characteristics of earthquakes are accurately reproduced in simulations, such as the duration, total movement, maximum speed, and stress change on the fault, even when model resolutions are not ideal. These exercises are important for promoting a new generation of advanced models for earthquakes. Understanding the sensitivity of simulation outputs will help test models against real‐world observations. Our community efforts can serve as a useful example to other geoscience communities. Key Points: We pursue community efforts to develop code verification benchmarks for three‐dimensional earthquake rupture and crustal faulting problemsWe assess the agreement and discrepancies of seismic and aseismic fault behavior among simulations based on different numerical methodsOur comparisons lend confidence to numerical codes and reveal sensitivities of model observables to major computational and physical factors [ABSTRACT FROM AUTHOR]

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. Plain Language Summary: Volcanic eruptions are noisy phenomena. During an eruption material is thrown into the atmosphere, pushing air out of the way and generating low frequency sound waves termed infrasound. We use infrasound observations to learn about the properties of volcanic eruptions. However, our understanding of the complex processes that generate sound during a volcanic eruption is limited. In order to address this, we perform simulations of volcanic eruptions and the associated infrasound signal. We compare our simulation results to an analytical model that is commonly used to interpret volcano infrasound observation. We show that for low exit velocities (up to 100 m/s or M ≈ 0.3 where M is the Mach number) the analytical model does a good job in explaining the infrasound observations and the radiation pattern. However, for higher exit velocities the analytical model overpredicts the peak amplitude of the infrasound signal, underpredicts the erupted volume, and does not account for the directionality of the radiation pattern. This work quantifies some of the complexities that should be considered when interpreting infrasound observations and is a step towards developing more sophisticated source models for volcanic eruptions. Key Points: Aeroacoustic simulations from the start‐up of a pressure‐balanced volcanic jet model fluid flow and infrasound radiationCompact monopole source model underpredicts the erupted volume for eruptions with sonic or supersonic exit velocitiesInfrasound radiation pattern depends on jet dynamics and is highly anisotropic for eruptions with sonic or supersonic exit velocities [ABSTRACT FROM AUTHOR]

It is widely recognized that fluid injection can trigger aseismic fault slip. However, the processes by which the fluid‐rock interactions facilitate or inhibit slip are poorly understood and some are oversimplified in most models of injection‐induced slip. In this study, we perform a 2D anti‐plane 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 porosity and permeability changes that accompany slip, including dilatancy, and quantify how these processes affect pore pressure diffusion, which couples to aseismic slip. 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 eventually triggers aseismic slip close to the injection site. In the second phase, aseismic slip front expands outward and dilatancy causes pore pressure to depart from the linear diffusion model. We quantify how prestress, injection rate, permeability and other fluid transport properties affect the slip front migration rate, finding rates ranging from 10 to 1,000 m/day for typical parameters. The migration rate is strongly influenced by the fault's closeness to failure and injection rate. The total slip on the fault, on the other hand, is primarily determined by the injected volume, with minimal sensitivity to injection rate. Additionally, we show that when dilatancy is neglected, slip front migration rate and total slip can be several times higher. Our modeling demonstrates that porosity and permeability evolution, especially dilatancy, fundamentally alters how faults respond to fluid injection. Plain Language Summary: The underground injection of fluids during wastewater disposal, geothermal operations, and other energy‐production activities has been linked to the occurrence of earthquakes. In addition to earthquakes, fluid injection can also trigger aseismic slip on faults, that is, frictional sliding that occurs so slowly that seismic waves and ground shaking are not produced. Here, we perform computer modeling of fluid injection and aseismic slip, exploring how the injection rate and fluid transport properties influence the aseismic slip response. We speculate that additional complexity in frictional properties and other conditions would cause aseismic slip to be accompanied by numerous, small earthquakes (microseismicity), as is often observed during injection. We quantify the rate at which aseismic slip migrates outward from the injection site and compare predicted migration rates to observed microseismicity patterns. Our model also predicts fluid pressure changes, slip, rock deformation surrounding the fault, and fluid flow paths that might be measurable and used to validate the modeling. Key Points: Modeling of constant rate fluid injection into a fault predicts steadily propagating aseismic slip frontMigration rate of aseismic slip front increases with injection rate and ranges from 10 to 1,000 m/day for typical parametersDilatancy and permeability enhancement alter system response as compared to linear pore pressure diffusion [ABSTRACT FROM AUTHOR]

Localized frictional sliding on faults in the continental crust transitions at depth to distributed deformation in viscous shear zones. This brittle‐ductile transition (BDT), and/or the transition from velocity‐weakening (VW) to velocity‐strengthening (VS) friction, are controlled by the lithospheric thermal structure and composition. Here, we investigate these transitions, and their effect on the depth extent of earthquakes, using 2D antiplane shear simulations of a strike‐slip fault with rate‐and‐state friction. The off‐fault material is viscoelastic, with temperature‐dependent dislocation creep. We solve the heat equation for temperature, accounting for frictional and viscous shear heating that creates a thermal anomaly relative to the ambient geotherm which reduces viscosity and facilitates viscous flow. We explore several geotherms and effective normal stress distributions (by changing pore pressure), quantifying the thermal anomaly, seismic and aseismic slip, and the transition from frictional sliding to viscous flow. The thermal anomaly can reach several hundred degrees below the seismogenic zone in models with hydrostatic pressure but is smaller for higher pressure (and these high‐pressure models are most consistent with San Andreas Fault heat flow constraints). Shear heating raises the BDT, sometimes to where it limits rupture depth rather than the frictional VW‐to‐VS transition. Our thermomechanical modeling framework can be used to evaluate lithospheric rheology and thermal models through predictions of earthquake ruptures, postseismic and interseismic crustal deformation, heat flow, and the geological structures that reflect the complex deformation beneath faults. Key Points: Earthquake modeling predicts viscous shear zones a few kilometers wide below faultsShear heating creates a thermal anomaly of tens to hundreds of degrees in the lower crustShear heating raises the brittle‐ductile transition and can limit rupture depth [ABSTRACT FROM AUTHOR]

We quantify sliding stability and rupture styles for a horizontal interface between an elastic layer and stiffer elastic half‐space with a free surface on top and rate‐and‐state friction on the interface. This geometry includes shallowly dipping subduction zones, landslides, and ice streams. Specific motivation comes from quasiperiodic slow slip events on the Whillans Ice Plain in West Antarctica. We quantify the influence of layer thickness on sliding stability, specifically whether a steadily loaded system produces steady sliding or stick‐slip sequences. We do this using both linear stability analysis and nonlinear earthquake sequence simulations. We restrict our attention to the 2‐D antiplane shear problem but anticipate that our findings generalize to more complex 2‐D in‐plane and 3‐D problems. Steady sliding with velocity‐weakening rate‐and‐state friction is linearly unstable to Fourier mode perturbations having wavelengths greater than a critical wavelength (λc). We quantify the dependence of λc on the rate‐and‐state friction parameters, elastic properties, loading, and the layer thickness (H). Confirming previous studies, we find that λc ∝ H1/2 for small H and is independent of H for large H. The linear stability analysis provides insight into nonlinear earthquake sequence dynamics of a nominally velocity‐strengthening interface containing a velocity‐weakening region of width W. Sequence simulations reveal a transition from steady sliding at small W to stick‐slip events when W exceeds a critical width (Wcr), with Wcr ∝ H1/2 for small H. Overall, this study demonstrates that the reduced stiffness of thin layers promotes instability, with implications for sliding dynamics in thin layer geometries. Key Points: We quantified sliding stability of an elastic layer over a stiffer half‐space with rate‐and‐state friction at the interfaceReduced stiffness of thin layers with a free surface on top promotes instabilityThe linear stability analysis explains conditions for the onset of slow slip sequences [ABSTRACT FROM AUTHOR]

Very long period (VLP) seismic events (with dominant periods of 15 to 40 s), observed from 2007 to 2018 at the summit of Kīlauea Volcano, Hawai'i, arise from resonant oscillations in the shallow magma plumbing system. Utilizing an oscillation model developed in the companion paper (Liang et al., 2020), we perform Bayesian inversions on seismic data from four representative VLP events separately for the parameters of the shallow conduit‐reservoir system, exploring both sphere and crack reservoir geometries. Both sphere and crack geometries are preferentially located ∼1–2 km beneath the northeast edge of Halema'uma'u crater and produce similar fits to the data. Considering a reasonable range for reservoir storativity, magma density, and density contrast between the top and bottom of the conduit, we favor a spherical reservoir with a radius of 0.8 to 1.2 km and a short conduit of less than a few hundred meters. For this geometry, buoyancy from density stratification in the conduit provides the dominant restoring force for the VLP oscillation. Viscosity is constrained within an order of magnitude for each event (e.g., approximately 2 to 23 Pa s for one event versus 27 to 513 Pa s for another). Changes in VLP period T and quality factor Q can be explained by changes in viscosity, density stratification, and/or conduit/reservoir geometry. In particular, observed fluctuations in Q over short time intervals (e.g., hours) with minimal changes in T apparently require rapid changes of magma viscosity by over an order of magnitude, assuming geometry remains unchanged, possibly reflecting changes in volatile content, bubble concentration, or conduit flow regime. Key Points: Inversions of VLP seismic data favor a short conduit and spherical reservoir over a single crackBuoyancy instead of reservoir stiffness is likely to dominate restoring force of VLP oscillationsChanges in VLP period and decay rate suggest time‐varying viscosity and density stratification [ABSTRACT FROM AUTHOR]

Very long period (VLP, 2–100 s) seismic signals at basaltic volcanoes like Kilauea, Hawai'i, and Mount Erebus, Antarctica are likely from resonant oscillations of magma within the shallow plumbing system. The system consists of conduits connected to cracks (dikes and sills) or reservoirs of other shapes. A quantitative understanding of wave propagation and resonance in a coupled conduit‐crack system is required to interpret observations. In this work, we idealize the system as an axisymmetric conduit coupled to a tabular crack, accounting for fluid inertia, compressibility, and viscosity, buoyancy, and crack wall elasticity. We perform time domain simulations and eigenmode analyses of the governing equations, linearized about a rest state. The fundamental mode or conduit‐reservoir mode reflects the balance of conduit magma inertia with buoyancy (and, for small cracks, crack wall elasticity). Magma oscillates in an effectively incompressible manner within the conduit, deflating and inflating the crack, which couples to the surrounding solid to produce observable surface displacements. For sufficiently low viscosity magmas, viscous effects are confined to boundary layers. Shorter period modes are primarily reverberating crack waves with negligible coupling to the conduit. Finally, we introduce an approximate reduced model for the conduit‐reservoir mode, which can also handle more general reservoir geometries (e.g., spherical chambers). The reduced model connects the observable VLP period and quality factor to two uniquely constrained parameters: the inviscid oscillation period T0 and the viscous diffusion time τvis across the conduit radius. Our models can be extended to study the seismic response of more complex magmatic systems. Key Points: Eigenmodes of coupled conduit‐crack system explain VLP seismic events and crack wave resonancesConduit‐reservoir mode with period of tens of seconds features magma oscillating in conduit driven by buoyancy and reservoir stiffnessReduced model connects VLP period and quality factor to geometry and magma properties and highlights parameter trade‐offs [ABSTRACT FROM AUTHOR]