Simulating Subterranean Fluid Injection through Iteration on the VirtualQuake Model

Abstract: This work extends the VirtualQuake earthquake simulation framework to incorporate the effects of fluid injection on fault stability and induced seismicity. Reworking VirtualQuake into a system using stress point sources, instead of rectangular segments, the new model offers increased geometric flexibility, greater stability, and the re-addition of cross-fault interactions. This approach is paired with the addition of fluid injection modeling, through the distribution of inflationary stress sources, according to invasion percolation, simulating both the stress effect of the injection on nearby faults, and deformation from the injection itself. The model captures both immediate and long-term impacts of injection cycles, including hydraulic fracturing processes and post-injection pressure dissipation. Results show that while single injections produce limited stress changes, repeated injections generate persistent high-pressure regions that progressively destabilize nearby faults, increasing the likelihood of seismic events. This model evolution offers a tool for the evaluation and characterization of long term risks from commercial injection.

• The paper introduces a novel point source tiling approach that replaces rectangular patches to accurately simulate complex fault geometries.
• It integrates invasion percolation and advanced pressure dynamics to capture the effects of subsurface fluid injection on seismic hazards.
• Convergence studies confirm that the new method achieves computational efficiency while preserving key physical insights in multi-cycle injection simulations.

Simulating Subterranean Fluid Injection through Iteration on the VirtualQuake Model

Introduction and Motivation

This paper systematically extends earthquake simulation frameworks to include the effects of anthropogenic subsurface fluid injection, with particular emphasis on hydraulic fracturing, carbon sequestration, geothermal energy extraction, and underground hydrogen storage. All these technologies alter subsurface stress fields, with induced seismicity becoming an increasingly critical hazard for both hydrocarbon production and emergent green technologies. Existing approaches such as VirtualQuake capture physics-based stress interactions and rupture cascades, but their original implementations rely on discretized rectangular fault patches, introducing numerical instabilities and hampering flexible modeling of realistic, intersecting, or highly irregular faults systems.

Technical Innovations: Point Source Tiling and Fractal Fault Representation

The primary technical advancement achieved is the replacement of rectangular fault patch solutions with a tiled collection of elastic half-space point sources. This substantially improves geometric flexibility: faults can assume arbitrary or even fractal geometries, intersections are naturally accommodated, and the singularities at patch boundaries are eliminated. Modular node-based object-oriented code architecture enables arbitrary node implementations, fostering extensibility and composability.

Convergence studies demonstrate that a relatively sparse tiling of point sources attains accurate stress fields (Figures 1 and 2), with diminishing returns beyond modest increases in source density. This allows computational efficiency while preserving physical accuracy.

Figure 1: The convergence of a tiled point source system with increasing source density.

Figure 2: Percent difference in magnitude between a tiled fault with 4 and 16 sources, respectively.

To further represent fault surface complexity, the method supports sampling from two-dimensional fractional Brownian motion (FBM) fields to generate fault surface perturbations, parameterized by the Hurst exponent. This allows statistically accurate modeling of the observed fractal nature of fault topography, significant for its influence on slip behavior and permeability variations (Figures 3 and 4).

Figure 3: Behavioral differences of Fractional Brownian Motion with varying Hurst Parameter.

Figure 4: The result of a two-dimensional Brownian transformation on a sample fault geometry.

Physics-Based Earthquake Simulation Mechanism

The adapted modeling formalism follows a robust event-driven earthquake cycle approach. For each fault node, projected shear and normal stresses are computed using Green's functions, stored as matrices, and accumulated across time linearly according to predefined slip rates. Rupture nucleation is controlled by the Coulomb Failure Function (CFF), incorporating static friction and hydrostatic pore pressure. Upon CFF transiting zero, an earthquake cascade is triggered, with dynamic slip evolution determined by self-stiffness and seismic stress drop parameters. This approach enables accurate modeling of stress transfer and emergent earthquake moment magnitude, with explicit handling of stress shadowing and aftershock phenomena.

Incorporating Subsurface Fluid Injection: Invasion Percolation and Pressure Dynamics

A central contribution is the explicit simulation of fluid injection effects. Fluid front evolution is modeled via invasion percolation in a 3D lattice, capturing the cluster bursts and fractal growth characteristic of hydraulically driven fracture propagation. Local bond strengths can be biased directionally, reflecting increased permeability along previously activated fault planes.

Pressure evolution of injected fluid follows both Darcy and non-Darcy (Forchheimer) regimes depending on injection velocity, and each invaded site is treated as an inflationary point source. This allows fluid-induced stress perturbations to be integrated directly, modifying the CFF of nearby fault sections and thus their seismic hazard. The pressure distribution over time and space is computed with consideration of complex boundary and material heterogeneity, and the algorithm flexibly supports time-dependent injection schedules (Figure 5).

Figure 5: Pressure dynamics at the beginning of a series of fracking cycles.

Fracture Propagation, Multi-Stage Injection, and Long-Tailed Pressure Relaxation

The simulation accurately captures multi-cycle fracking operations, including distinct breakdown, sustained injection, and shut-in periods. The temporal evolution includes pore pressure dissipation governed by fracture volume, fluid and rock properties, and proppant retention. Notably, the model demonstrates that a single fracking operation leaves elevated pressures for extended periods (often over a year in low-permeability shale), with pressure relaxation governed by slow diffusive flow into the surrounding rock (Figures 6 and 7).

Figure 6: An example injection simulation, with associated spatial pressure drop through the fracture.

Figure 7: The beginning of pressure equalization over the course of a year, post injection.

This long-lived elevation of fluid pressure drives cumulative increases in local seismic hazard, as subsequent injections occur long before conditions can return to baseline. Simulation of 10 sequential injection cycles clearly shows monotonically increasing CFF on fault segments, with spatial patterning reflecting both fluid migration and geometry of proximal fault elements (Figure 8).

Figure 8: The cumulative effects on CFF from a sample fluid injection. Each voxel square is mapped to an inflationary source.

Implications and Future Directions

The extended VirtualQuake approach enables direct mapping from injection schedules and site properties to predicted stress evolutions and seismicity. The framework’s explicit physics and geometric flexibility position it as a platform for:

• Investigating operational changes (e.g., extended pauses between injections) and their potential for seismic risk mitigation.
• Generating synthetic datasets for data-driven risk assessment tools and seismicity forecasts.
• Parameterizing models with in-situ measurement to improve regional hazard assessments.

The implementation can be readily extended with additional physical modules (e.g., poroelasticity, time-dependent permeability, coupled chemical effects), and the modular node architecture streamlines incorporation of multi-physics interactions. Bridging this mechanistic modeling with machine learning surrogates is a natural progression, leveraging the simulation outputs as training corpora for rapid risk estimation.

Conclusion

This work presents a comprehensive, extensible simulation framework that accurately models subsurface fluid injection and its seismic implications by generalizing stress representations to point source tilings and incorporating invasion percolation for fluid migration. The approach enables physically rigorous and computationally efficient assessment of injection-induced seismicity, supporting both operational decision-making and theoretical advances in earthquake science.