Abstract Seismic ruptures that extend over kilometers along faults originate in shear zones only millimeters thick. We present a hybrid numerical framework coupling the spectral boundary integral method (SBIM) for the poroelastic bulk with a staggered‐grid finite difference method (FDM) for the fault gouge. This SBI‐FD approach enables cross‐scale simulations, resolving large‐scale slip dynamics while capturing pore pressure gradients and strain localization across narrow gouge layers. The pore pressure evolves within the gouge and interacts with the surrounding bulk in both the along‐fault and across‐fault directions. Using this framework, we investigate how bulk properties affect rupture dynamics, focusing on the influence of bulk poroelasticity on shear localization and pore pressure evolution within fault zones. We then examine whether the poroelastic bulk can be approximated by a reduced boundary‐layer approximation. Our simulations indicate that bulk poroelasticity can substantially influence the migration of shear localization, which may affect how geological records are interpreted. Assuming an undrained condition provides a good approximation of poroelastic bulk behavior for the rapid rupture studied in the paper. Additionally, we capture the migration of localized layers during fault slip and examine their significance for pseudotachylyte formation in natural fault zones.
Wang et al. (Sun,) studied this question.