3D wave propagation and earthquake dynamic rupture simulations in complex poroelastic media

Numerical simulations of earthquakes and seismic wave propagation require accurate material models of the solid Earth. In contrast to purely elastic rheology, poroelasticity accounts for pore fluid pressure and fluid flow in porous media. Poroelastic effects can alter both the seismic wave field and the dynamic rupture characteristics of earthquakes. For example, the presence of fluids may affect cascading multi fault ruptures, potentially leading to larger-than-expected earthquakes. However, incorporating poroelastic coupling into the elastodynamic wave equations increases the computational complexity of numerical simulations compared to elastic or viscoelastic material models, as the underlying partial differential equations become stiff. In this study, we use a Discontinuous Galerkin solver with Arbitrary High-Order DERivative time stepping (ADER-DG) of the poroelastic wave equations implemented in the open-source software SeisSol to simulate 3D complex seismic wave propagation and 3D dynamic rupture in poroelastic media. We verify our approach for double couple point sources using independent methods, including a semi-analytical solution and a finite difference scheme, and a homogeneous full-space and a poroelastic layer-over-half-space model, respectively. In a realistic carbon capture and storage (CCS) reservoir scenario at the Sleipner site in the Utsira Formation, Norway, we model 3D wave propagation through poroelastic sandstone layers separated by impermeable shale. Our results show a sudden change in the pressure field across material interfaces, which manifests as a discontinuity when viewed at the length scale of the dominant wavelengths of S- or fast P-waves. Accurately resolving the resulting steep pressure gradient dramatically increases the computational demands, requiring high-resolution modeling. We show that the Gassmann elastic equivalent model yields almost identical results to the fully poroelastic model when focusing solely on solid particle velocities. We extend this approach using suitable numerical fluxes to 3D dynamic rupture simulations in complex fault systems, presenting the first 3D scenarios that combine poroelastic media with geometrically complex, multi fault rupture dynamics and tetrahedral meshes. Our findings reveal that in contrast to modeling wave propagation only, poroelastic materials significantly alter rupture characteristics compared to using elastic equivalent media since the elastic equivalent fails to capture the evolution of pore pressure. Particularly in fault branching scenarios, the Biot coefficient plays a key role in either promoting or inhibiting fault activation. In some cases, ruptures are diverted to secondary faults, while in others, poroelastic effects induce rupture arrest. In a fault zone dynamic rupture model, we find poroelasticity aiding pulse-like rupture. A healing front is induced by the reduced pore pressure due to reflected waves from the boundaries of the poroelastic damage zone. Our results highlight that poroelastic effects are important for realistic simulations of seismic waves and earthquake rupture dynamics. In particular, our poroelastic simulations may offer new insights into the complexity of multi-fault rupture dynamics, fault-to-fault interaction, and seismic wave propagation in realistic models of the Earth's subsurface.