ABSTRACT We derive a novel seismic representation theorem for nonelastic dynamics in an arbitrarily shaped 3D volume enclosed by a surface S, building on Burridge et al. (1993) and extending Coppess et al. (2022). The novel result is the generalized seismic moment tensor with nine independent components, Mji(t)=∫S[λδijuknk+μ(uinj+ujni)−xj′ΔTi]dA, the antisymmetric part of which gives rise to the seismic torque, Nl(t)=−ϵlji∫Sxj′ΔTidA. The expressions are written in terms of surface integrals of perturbation tractions ΔTi (exerted on the source region by the surrounding solid), elastic displacements ui, lever arms xi′ extending from the source centroid to S, and elastic moduli λ and μ in an isotropic solid. Here, δij is the Kronecker delta, ϵlji the Levi-Civita symbol, and nj the outward unit normal on S. The expressions are valid for arbitrary rheology inside the source region. When combined with appropriate Green’s functions, the expressions provide a rigorous connection between nonelastic dynamics inside the source region and the displacement (at wavelengths large compared with the source dimensions) outside the source region. We apply the theorem to dynamic rupture simulations of basaltic caldera collapse earthquakes, a complex source process involving ring faulting, caldera block subsidence, and the mechanical coupling between the caldera block and the subcaldera magma chamber, yielding a time-dependent force, moment, and torque representing the dynamics. The results establish the correct interpretations of seismic representations in terms of the dynamics of collapse earthquakes. We empirically validate the theorem by demonstrating that the surface wavefield generated by the seismic representation converges to that of the dynamic rupture simulation, in the long wavelength limit. The theorem enables an efficient forward modeling approach when the long-period wavefield of exotic geophysical phenomena is simulated independently of the source dynamics, significantly reducing the complexity of modeling.

Physics-based simulations are critical for understanding natural hazards. The increasing complexity of numerical codes requires benchmark exercises to verify that different computational methods yield consistent results when solving the same governing equations. Here, we present an open-access web platform designed for the verification of earthquake dynamic rupture, seismic cycle, and tsunami simulations. The platform architecture utilizes a modular, serverless backend on Amazon Web Services (AWS) to provide scalable file processing and visualization. A lightweight static web application provides a secure interface for uploading and managing results, while the browser-based data visualization enables interactive analysis of time series and surface grid data. By using structured JavaScript Object Notation (JSON) text files to define benchmark structures, the system remains fully extensible, allowing the addition of new scenarios without modifying the underlying software logic. The platform hosts the "The Tsunami Problem Versions" (TTPV) 1 & 2, two benchmarks for 3D fully coupled earthquake dynamic rupture and tsunami generation, and provides a framework for earthquake cycle models. This community resource aims to build trust in numerical simulations and facilitate long-term collaborative code verification as modeling software continues to evolve.

Hurricane evolution is affected by turbulence in the hurricane boundary layer (HBL), which is typically measured using aircraft flights and towers. Through a case study of a landfalling hurricane, we show that seismoacoustic data can also be used for HBL turbulence analysis. We identified contributions of HBL turbulence in infrasound pressure and seismic displacement, validating our interpretation by combining large-eddy simulation, calibrated with meteorological data, with quasi-static elastic deformation modeling. The convection velocity of the turbulent pressure field is key to this pressure-displacement coupling. For atmospheric studies, continuous infrasound pressure serves as a proxy for 10-meter wind speed, and the inertial subrange of pressure spectra provides an estimate of the turbulent dissipation rate near the top of the surface layer at ~100 to 200 meters, complementing portable tower data at ~10 meters.

Relative plate motion in subduction zones transitions from frictional slip to viscous flow with increasing depth and temperature. The frictional-viscous transition can control the depth extent of megathrust earthquakes and episodic tremor and slip (ETS). Pore fluid pressure is a critical control on the transition, but models for its depth dependence are lacking. Here, we present a steady-state modeling framework to calculate the fluid pressure and shear stress profile along the subduction interface. We consider fluid production from dehydration reactions in the subducting oceanic lithosphere, calculated from thermodynamic equilibrium models. These fluids are channeled updip, though in some models we allow for fluid loss into the overriding plate. The fluid pressure is calculated from Darcy's law, with permeability depending on effective stress, temperature, and slip rate. We allow for both rate-state frictional sliding and thermally activated linear viscous flow, and solve for the partitioning of deformation between them. We apply the modeling framework to the Cascadia subduction zone. Our results show nearly uniform effective stress in the seismogenic zone, below which it decreases with depth and fluid pressure approaches lithostatic pressure. The frictional-viscous transition spans a wide range of depths, and mixed frictional-viscous deformation is predicted at the ETS source depth. A model with fluid leak-off into the overlying plate produces a more heterogeneous effective stress, with a local minimum near the major dehydration depth. Our results provide important insight into earthquake hazards and the mechanism of ETS in Cascadia, and the modeling framework is applicable to global subduction zones.

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.

The intensity of explosive volcanic eruptions is correlated with the amplitude of eruption tremor, a ubiquitously observed seismic signal during eruptions. Here we expand upon a recently introduced theoretical model that attributes eruption tremor to particle impacts and dynamic pressure changes in the turbulent flow above fragmentation (Gestrich et al., 2020). We replace their point source model with Rayleigh wave Green's functions with full Green's functions and account for depth variation of input fields using conduit flow models. The latter self-consistently capture covariation of input fields like particle velocity, particle volume fraction, and density. Body wave contributions become significant above 2-3 Hz, bringing the power spectral density (PSD) closer to observations. Conditions at the vent are not representative of flow throughout the tremor source region and using these values overestimates tremor amplitude. Particle size and its depth distribution alter the PSD and where dominant source contributions arise within the conduit. Solutions with decreasing mass eruption rate, representing a waning eruption, reveal a shift in the dominant tremor contribution from turbulence to particle impacts. Our work demonstrates the ability to integrate conduit flow modeling with volcano seismology studies of eruption tremor, providing an opportunity to link observations to eruptive processes.

Many geo-engineering applications, such as hydrocarbon extraction, geothermal energy production, and geologic carbon storage, involve the injection and/or extraction of fluids from the subsurface. A key hazard associated with these activities is induced seismicity. Fluid injection or extraction can perturb the stress state of the surrounding rock formation, potentially triggering the sliding of critically oriented faults and resulting in seismic events. Therefore, it is crucial to develop numerical models that can accurately capture reservoir behavior, including complex fluid migration, the mechanical response of the surrounding rock mass, and the physics of earthquake nucleation. In this work, we present a novel strategy, implemented in the GEOS simulation framework, for coupled poromechanical and quasi-dynamic earthquake simulations, incorporating fault frictional behavior described by a rate- and state-dependent friction law. The poromechanical equations are discretized using a low-order finite element method for the mechanical response, coupled with a finite volume method for fluid flow. Faults are modeled as lowerdimensional manifolds and discretized using surface elements at the interfaces between 3D matrix cells. Contact constraints are enforced using face-based, piecewise-constant Lagrange multipliers to represent fault tractions, and stabilization is achieved through face-based bubble functions that enrich the displacement space. Slip velocities, slip and the state evolution variable are considered to be piecewise-constant on fault elements. The quasi-static poromechanical equations are coupled with a quasi-dynamic earthquake model using a split-operator approach. At each timestep, the poromechanical equations are first solved under the assumption of fixed fault slip, after which the resulting fault tractions are used as input for the quasi-dynamic problem. The earthquake model is then solved locally for each fault element. We validate our approach against benchmark problems developed by the Sequences of Earthquakes and Aseismic Slip (SEAS) working group, supported by the Statewide California Earthquake Center (SCEC). Results demonstrate the robustness and accuracy of our method in capturing the complex interactions between fluid flow, fault mechanics, and earthquake nucleation processes. This work was partially performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Labo-ratory under Contract DE-AC52-07NA27344.

On 26 September 2022 two seismic events near the Danish island of Bornholm in the Baltic Sea were detected. The first event with a magnitude Mw 2.3 occurred at 00:03 UTC 40 km east-southeast of Bornholm. The determined location and the origin time of the event are consistent with data of the pressure decrease on one of the Nord Stream 2 pipelines. Another sequence of events occurred 17 hours later at 17:03 UTC around 60 km north-east of Bornholm with a maximum magnitude of Mw 2.7. It consists of three closely successive, but separable, single events. Using relative localisation methods and the gas pressure inside the pipeline recorded at the landing site in Germany, we can assign the epicentres of the three events to the locations of the leaks in the pipelines of Nord Stream 1 and 2. Based on comparable events in the region, which include both tectonic earthquakes and explosions, the explosive character of the investigated Nord Stream events can be verified. Infrasound signals associated with the destruction of the Nord Stream pipelines were recorded at two stations (I26DE in the Bavarian Forest and IKUDE near Kühlungsborn) in Germany. Particularly after the event sequence at 17:03 UTC, distinctive signals were registered whose characteristics indicate an explosive event with subsequent gas leakage at the surface. Our modelling of the sources shows that the measured seismic signals can sufficiently be explained by the instantaneous gas release. Synthetic seismograms for such a source and a subsurface model adapted for the study area show high consistency with the measured signals. Based on the released energy and the characteristics of the recorded waveforms, we conclude that the impulsive gas release from the burst gas pipes constitutes the dominant part of the signal source. The model places an upper limit of approximately 50 kg TNT equivalent on the yield of the chemical explosive component of the events, but we note that smaller yields may also be consistent with the data. We also carried out an analysis of the seismic signals of the event on the Balticconnector pipeline between Finland and Estonia on 8 October 2023 and found that again the instantaneous gas release can sufficiently explain the observed data. This supports a possible mechanical cause of the damage.

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 quasi‐dynamic benchmark problems BP6‐QD‐A/S/C that consider an 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 (‐A) and slip (‐S) law formulations for frictional state evolution, respectively, allowing us to better understand how various computational factors across codes affect the simulated evolution of pore pressure and aseismic slip. Comparisons of problems using the aging versus slip law, and a constant friction coefficient (‐C), illustrate how aseismic slip models 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 attained by ensuring sufficiently fine time‐stepping and consistent treatment of 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.

Observations of fluid-driven swarm seismicity expanding with the same diffusive space-time behavior as analytical solutions for aseismic slip have been interpreted as evidence that stress changes from aseismic slip trigger seismic slip. In some cases, aseismic slip is confirmed from crustal deformation measurements or sheared wellbore casing. Our work offers another interpretation of migrating seismicity that may be relevant when there is no independent evidence for aseismic slip. We present 2D earthquake sequence models that simulate constant-pressure fluid injection into a velocity-weakening rate-and-state fault, permeability enhancement with slip, and fluid transport. We also derive an analytical expression for a unilateral fluid-driven aseismic shear crack with constant friction and permeability enhancement. Our numerical results show that pressure diffusion and elastic stress transfer from seismic slip drive seismicity fronts that expand diffusively outward from the injector, and that analytical solutions for aseismic slip match this seismicity pattern, even when almost all slip is seismic. Aseismic slip at the slip front increases permeability, allowing fluid influx and pressurization, weakening the fault prior to slipping seismically. Adding roughness to the fault introduces heterogeneous fault normal stress, resulting in a broader distribution of event sizes, distributed seismicity behind the slip front, and swarm-like clusters of microseismicity. Even with this complexity, the nonplanar simulations produce a diffusively expanding slip front and mostly seismic slip. Our results suggest that aseismic slip solutions can be used to quantitatively interpret the space-time behavior of migrating swarms, even in cases with negligible measured aseismic slip.

The seismic hum in the -20 - 300 s period band is usually explained by the primary mechanism, where ocean infragravity waves exert pressure changes directly on the seafloor. However, there are some indications that atmospheric processes might also contribute to the seismic hum band. Hurricane landfall provides a unique opportunity to investigate the strong seismic ambient noise generated by turbulent pressure fluctuations from the atmosphere. We revisit Hurricane Isaac in 2012 which passed through the Transportable Array (TA) stations after its landfall on the Louisiana coast. Taking advantage of data recorded by co-located pressure sensors and seismometers, we propose the usage of the continuous wavelet transform to perform high-resolution timefrequency analysis, which reveals a clear signature of the hurricane eye and its surrounding eyewall in both pressure and seismic data. Our spectral analysis also captures diurnal cycles in atmospheric noise. Wavelet coherence analysis between pressure and vertical displacement shows a separation of two noise spectral bands dominated by the ocean and atmosphere, respectively, and we focus on the atmospheric noise with period 20 - 100 s, a subset of the seismic hum band. Observations consistently reveal high coherence between pressure and vertical displacement for the atmospheric noise, which suggests a local (-100 m - 1 km) scale coupling, in contrast to previous modeling that explored a non-local mechanism involving seismic wave generation from the entire hurricane. To test the local coupling hypothesis, we perform hurricane-scale numerical modeling to quantify the seismic contributions from the surface pressure fluctuations at different parts of the hurricane. We integrate surface wind re-analysis data and empirical relations derived from turbulence studies into the description of our input pressure source. Our results demonstrate the local quasi-static nature of the atmospheric noise from the hurricane, different from the previous model which consists of propagating seismic waves. Our modeling paves the way for further local-scale modeling of atmospheric noise based on realistic turbulent pressure fields. The combination of data and knowledge from both seismology and atmospheric sciences is essential to advancing the understanding and applications of atmosphere-generated seismic ambient noise.

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.

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 Kilauea 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 Kilauea caldera collapse of 2018.

Injection-induced seismicity and aseismic slip often involve the reactivation of long-dormant faults, which may have extremely low permeability prior to slip. In contrast, most previous models of fluid-driven aseismic slip have assumed linear pressure diffusion in a fault zone of constant permeability and porosity. Slip occurs within a frictional shear crack whose edge can either lag or lead pressure diffusion, depending on the dimensionless stress-injection parameter that quantifies the prestress and injection conditions. Here, we extend this foundational work by accounting for permeability enhancement and dilatancy, assumed to occur instantaneously upon the onset of slip. The fault zone ahead of the crack is assumed to be impermeable, so fluid flow and pressure diffusion are confined to the interior, slipped part of the crack. The confinement of flow increases the pressurization rate and reduction of fault strength, facilitating crack growth even for severely understressed faults. Suctions from dilatancy slow crack growth, preventing propagation beyond the hydraulic diffusion length. Our new two-dimensional and three-dimensional solutions can facilitate the interpretation of induced seismicity data sets. They are especially relevant for faults in initially low permeability formations, such as shale layers serving as caprock seals for geologic carbon storage, or for hydraulic stimulation of geothermal reservoirs.This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

We present an adjoint-based optimization method to invert for stress and frictional parameters used in earthquake modeling. The forward problem is linear elastodynamics with nonlinear rate-and-state frictional faults. The misfit functional quantifies the difference between simulated and measured particle displacements or velocities at receiver locations. The misfit may include windowing or filtering operators. We derive the corresponding adjoint problem, which is linear elasticity with linearized rate-and-state friction with time-dependent coefficients derived from the forward solution. The gradient of the misfit is efficiently computed by convolving forward and adjoint variables on the fault. The method thus extends the framework of full-waveform inversion to include frictional faults with rate-and-state friction. In addition, we present a space-time dual-consistent discretization of a dynamic rupture problem with a rough fault in antiplane shear, using high-order accurate summation-by-parts finite differences in combination with explicit Runge--Kutta time integration. The dual consistency of the discretization ensures that the discrete adjoint-based gradient is the exact gradient of the discrete misfit functional as well as a consistent approximation of the continuous gradient. Our theoretical results are corroborated by inversions with synthetic data. We anticipate that adjoint-based inversion of seismic and/or geodetic data will be a powerful tool for studying earthquake source processes; it can also be used to interpret laboratory friction experiments.

We introduce an energy stable, high-order-accurate finite difference approximation of the dynamic, pure bending Kirchhoff plate equations for complex geometries and spatially variable properties. We utilize the summation-by-parts (SBP) framework to discretize the biharmonic operator with variable coefficients, with attention given to free and clamped boundary conditions and corner conditions. Energy conservation is established by combining SBP boundary closures with weak enforcement of the boundary and interface conditions using a penalty (simultaneous approximation term, SAT) technique. Then we couple the plate equations to the shallow water equations to study flexural-gravity wave propagation, and prove that the semi-discrete system of equations is self-adjoint. We demonstrate the stability and accuracy properties of the method on curvilinear multiblock grids using the method of manufactured solutions. The method, which we provide in an open-source code, is then used to model ocean wave interactions with the Thwaites Glacier and Pine Island Ice Shelf in the Amundsen Sea off the coast of West Antarctica.

Numerical modeling of earthquake dynamics and derived insight for seismic hazard relies on credible, reproducible model results. The SEAS (Sequences of Earthquakes and Aseismic Slip) initiative has set out to facilitate community code comparisons, and verify and advance the next generation of physics-based earthquake models that reproduce all phases of the seismic cycle. With the goal of advancing SEAS models to robustly incorporate physical and geometrical complexities, here we present code comparison results from two new benchmark problems: BP1-FD considers full elastodynamic effects and BP3-QD considers dipping fault geometries. Eight modeling groups participated in each benchmark, allowing us to explore these physical ingredients across multiple codes and better understand associated numerical considerations. We find that numerical resolution and computational domain size are critical parameters to obtain matching results, with increasing domain-size requirements posing challenges for volume-based codes even in 2D settings. Codes for BP1-FD implemented different criteria for switching between quasi-static and dynamic solvers, which require tuning to obtain matching results. In BP3-QD, proper remote boundaries conditions consistent with specified rigid body translation are required to obtain matching surface displacements. With these numerical and mathematical issues resolved, we obtain good agreement among codes in long-term fault behavior, earthquake recurrence intervals, and rupture features of peak slip rates and stress drops for both benchmarks. Including full inertial effects generates events with larger slip rates and rupture speeds compared to the quasi-dynamic counterpart. For BP3-QD, both dip angle and sense of motion (thrust versus normal faulting) alter ground motion on the hanging and foot walls, and influence event patterns, with some sequences exhibiting similar-sized characteristic earthquakes, and others exhibiting several earthquakes of differing magnitudes. These findings underscore the importance of considering full dynamics and non-vertical dip angles in SEAS models, as both influence short and long-term earthquake behavior, and associated hazards.

Fault bends are known to act as barriers to rupture propagation in many earthquakes. A recent compilation of the surface rupture traces of paleo earthquakes quantifies the probability of rupture termination as a function of the bend angle. To understand the physical basis of these unique statistics, we carry out 2D quasi-dynamic earthquake sequence simulations on a fault with either a restraining or releasing (double-) bend. The fault is loaded by steady sliding from an adjacent creeping region, rather than with backslip, together with a stress relaxation method to avoid unphysical stress buildup due to the curvature of faults. This relaxation is a proxy for unmodeled off-fault secondary faulting. This ensures the existence of a long-term steady earthquake cycle and stress field, the latter of which is determined by the balance between the slip-induced stress changes and relaxation terms. We quantify the influence of the bend on rupture propagation by computing passing ratio (i.e., the fraction of ruptures that propagate through the bend). Our simulations approximately reproduce the Biasi-Wesnousky empirical law for a wide range of parameters (e.g., bend width, stress relaxation time, background stress). Also, we reproduce geologic observations that the long-term slip rate of a fault has local minima at restraining bends, without assuming spatially varying loading rates using the backslip approach. Additionally, we find that restraining and releasing bends have different earthquake cycles. For restraining bends, many ruptures are arrested before reaching the center of the bend, and the passing ratio decreases with increasing the bend angle. For releasing bends, ruptures stop after passing the center of the bend, and the passing ratio abruptly drops from near unity to zero at around 20∘. This difference can be qualitatively explained by an energy balance approach. Our model has a potential to understand the seismogenesis of nonplanar faults and we can readily extend our model into 3D specific fault systems for hazard assessment.

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.