Katherine R Coppess

and 2 more

Fragmentation plays a critical role in eruption explosivity by influencing the eruptive jet and plume dynamics that may initiate hazards such as pyroclastic flows. The mechanics and progression of fragmentation during an eruption are challenging to constrain observationally, limiting our understanding of this important process. In this work, we explore seismic radiation associated with unsteady fragmentation. Seismic force and moment tensor fluctuations from unsteady fragmentation arise from fluctuations in fragmentation depth and wall shear stress (e.g., from viscosity variations). We use unsteady conduit flow models to simulate perturbations to a steady-state eruption from injections of heterogeneous magma (specifically, variable magma viscosity due to crystal volume fraction variations). Changes in wall shear stress and pressure determine the seismic force and moment histories, which are used to calculate synthetic seismograms. We consider three heterogeneity profiles: Gaussian pulse, sinusoidal, and stochastic. Fragmentation of a high-crystallinity Gaussian pulse produces a distinct very-long-period (VLP) seismic signature and associated reduction in mass eruption rate, suggesting joint use of seismic, infrasound, and plume monitoring data to identify this process. Simulations of sinusoidal injections quantify the relation between the frequency or length scale of heterogeneities passing through fragmentation and spectral peaks in seismograms, with velocity seismogram amplitudes increasing with frequency. Stochastic composition variations produce stochastic seismic signals similar to observed eruption tremor, though computational limitations restrict our study to frequencies less than 0.25 Hz. We suggest that stochastic fragmentation fluctuations could be a plausible eruption tremor source.

Yuyun Yang

and 1 more

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.

Yuyun Yang

and 3 more

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.

Yuyun Yang

and 1 more

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.

Lauren S. Abrahams

and 4 more

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.

Junle Jiang

and 18 more

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.

Leighton M Watson

and 5 more

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.

Katherine R Coppess

and 2 more

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.

Lauren S Abrahams

and 2 more

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 quasi-periodic slow slip events on the Whillans Ice Plain in West Antarctica. We quantify the influence of layer thickness on sliding stability, specifically whether steady loading of the system produces steady sliding or sequences of stick-slip events. We do this using both linear stability analysis and nonlinear earthquake sequence simulations. We restrict our attention to the 2D antiplane shear problem, but anticipate that our findings generalize to the more complex 2D in-plane and 3D 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 (Η). We find that λ_c is proportional to sqrt(Η) for small Η and independent of Η for large Η. 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 (W_cr), with W_cr proportional to sqrt(H) 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.

Kali Allison

and 1 more

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.