Finite Difference Based Wave Simulation in Fractured Porous Rocks

Biot’s theory provides a framework for computing seismic wavefields in fluid saturated porous media. Here we implement a velocity-stress staggered grid 2D finite difference algorithm to model the wave-propagation in poroelastic media. The Biot’s equation of motion are formulated using a finite difference algorithm with fourth order accuracy in space and second order accuracy in time. Seismic wave propagation in reservoir rocks is also strongly affected by fractures and faults. We next derive the equivalent media model for fractured porous rocks using the linear slip model and perform numerical simulations in the presence of fractured interfaces. As predicted by Biot’s theory a slow compressional wave is observed in the particle velocity snapshots. In the layered model, at the boundary, the slow P-wave converts to a P-wave that travels faster than the slow P-wave. We finally conclude by commenting on the major details of our results.

💡 Research Summary

This paper presents a comprehensive numerical framework for simulating seismic wave propagation in fluid‑saturated porous rocks, explicitly accounting for the presence of fractures. The authors begin by recalling Biot’s theory, which predicts three coupled wave modes in a poroelastic medium: a fast compressional wave (P‑fast), a slow compressional wave (P‑slow), and a shear wave (S). To solve Biot’s equations of motion, they develop a two‑dimensional staggered‑grid finite‑difference scheme that places particle velocities and stresses on interleaved grids. Spatial derivatives are approximated with fourth‑order central differences, while temporal integration uses a second‑order Leap‑frog (or explicit central‑difference) method. This combination yields high accuracy while preserving the explicit nature of the algorithm, which is essential for large‑scale seismic simulations. Stability is governed by the Courant‑Friedrichs‑Lewy (CFL) condition, and the time step is adaptively chosen to satisfy the most restrictive wave speed in the model. To suppress artificial reflections, a perfectly matched layer (PML) is implemented at all domain boundaries, reducing reflected energy to below 10⁻⁴ of the incident amplitude.

The second major contribution is the incorporation of fractures through the linear slip model. In this approach, a fracture is idealized as a thin interface that permits relative tangential displacement (slip) proportional to the shear traction across it. The proportionality constant, the slip coefficient β, encapsulates the combined effects of fracture aperture, surface roughness, and infill material. By homogenizing a periodic array of such interfaces, the authors derive effective elastic stiffness tensors C* and effective permeability tensors κ* that replace the original Biot parameters in the governing equations. This “equivalent medium” retains the same mathematical structure as the original poroelastic model, allowing the same finite‑difference solver to be used without any additional code complexity. The model parameters (fracture density, orientation, and slip coefficient) can be varied to explore a wide range of realistic fracture networks.

Two numerical experiments are presented. The first is a benchmark in a homogeneous, unfractured poroelastic block. An impulsive source generates both fast and slow compressional waves as well as a shear wave. Snapshots of particle velocity clearly show the slow P‑wave traveling at a markedly lower phase velocity and exhibiting strong attenuation, in full agreement with analytical predictions from Biot’s theory. The second experiment introduces a layered medium in which one layer contains a set of parallel fractures modeled by the linear slip approach. When the slow P‑wave impinges on the fractured interface, part of its energy is converted into a fast P‑wave that propagates ahead of the incident slow wave. The conversion efficiency depends on the contrast in effective stiffness and permeability across the interface, as well as on the slip coefficient. Additionally, shear wave amplitudes are amplified at the interface due to the concentration of tangential strain, illustrating the role of fractures as shear‑wave guides.

The authors discuss the numerical performance of their scheme. The fourth‑order spatial discretization reduces numerical dispersion, allowing accurate phase‑velocity representation even with relatively coarse grids (10–15 points per wavelength for the fast P‑wave). The second‑order temporal scheme provides sufficient stability for the long‑time simulations required in reservoir‑scale studies. Computational cost scales linearly with the number of grid points, and the explicit nature of the algorithm makes it well‑suited for parallel implementation on modern GPU clusters.

In the concluding section, the paper emphasizes that the combination of a high‑order staggered‑grid finite‑difference solver with a linear‑slip‑based equivalent medium model offers a powerful and flexible tool for investigating wave phenomena in fractured porous rocks. The successful reproduction of the slow compressional wave and its conversion to a fast wave at fractured boundaries validates both the numerical implementation and the underlying physical model. The authors suggest several avenues for future work: extending the framework to three dimensions, incorporating non‑linear slip behavior to capture shear‑failure processes, and calibrating the slip coefficient against laboratory measurements or field observations. Ultimately, this methodology could improve seismic imaging, reservoir characterization, and monitoring of hydraulic fracturing operations by providing more realistic predictions of wavefield behavior in complex subsurface environments.