Biot's equations describe the physics of hydromechanically coupled systems establishing the widely recognized theory of poroelasticity. This theory has a broad range of applications in Earth and biological sciences as well as in engineering. The numerical solution of Biot's equations is challenging because wave propagation and fluid pressure diffusion processes occur simultaneously but feature very different characteristic time scales. Analogous to geophysical data acquisition, high resolution and three dimensional numerical experiments lately redefined state of the art. Tackling high spatial and temporal resolution requires a high‐performance computing approach. We developed a multi‐ graphical processing units (GPU) numerical application to resolve the anisotropic elastodynamic Biot's equations that relies on a conservative numerical scheme to simulate, in a few seconds, wave fields for spatial domains involving more than 1.5 billion grid cells. We present a comprehensive dimensional analysis reducing the number of material parameters needed for the numerical experiments from ten to four. Furthermore, the dimensional analysis emphasizes the key material parameters governing the physics of wave propagation in poroelastic media. We perform a dispersion analysis as function of dimensionless parameters leading to simple and transparent dispersion relations. We then benchmark our numerical solution against an analytical plane wave solution. Finally, we present several numerical modeling experiments, including a three‐dimensional simulation of fluid injection into a poroelastic medium. We provide the Matlab, symbolic Maple, and GPU CUDA C routines to reproduce the main presented results. The high efficiency of our numerical implementation makes it readily usable to investigate three‐dimensional and high‐resolution scenarios of practical applications.

SUMMARY: The efficient and accurate numerical modelling of Biot’s equations of poroelasticity requires the knowledge of the exact stability conditions for a given set of input parameters. Up to now, a numerical stability analysis of the discretized elastodynamic Biot’s equations has been performed only for a few numerical schemes. We perform the von Neumann stability analysis of the discretized Biot’s equations. We use an explicit scheme for the wave propagation and different implicit and explicit schemes for Darcy’s flux. We derive the exact stability conditions for all the considered schemes. The obtained stability conditions for the discretized Biot’s equations were verified numerically in one-, two- and three-dimensions. Additionally, we present von Neumann stability analysis of the discretized linear damped wave equation considering different implicit and explicit schemes. We provide both the Matlab and symbolic Maple routines for the full reproducibility of the presented results. The routines can be used to obtain exact stability conditions for any given set of input material and numerical parameters.

Shear strain localization refers to the phenomenon of accumulation of material deformation in narrow slip zones. Many materials exhibit strain localization under different spatial and temporal scales, particularly rocks, metals, soils, and concrete. In the Earth's crust, irreversible deformation can occur in brittle as well as in ductile regimes. Modeling of shear zones is essential in the geodynamic framework. Numerical modeling of strain localization remains challenging due to the non‐linearity and multi‐scale nature of the problem. We develop a numerical approach based on graphical processing units (GPU) to resolve the strain localization in two and three dimensions of a (visco)‐hypoelastic‐perfectly plastic medium. Our approach allows modeling both the compressible and incompressible visco‐elasto‐plastic flows. In contrast to symmetric shear bands frequently observed in the literature, we demonstrate that using sufficiently small strain or strain rate increments, a non‐symmetric strain localization pattern is resolved in two‐ and three‐dimensions, highlighting the importance of high spatial and temporal resolution. We show that elasto‐plastic and visco‐plastic models yield similar strain localization patterns for material properties relevant to applications in geodynamics. We achieve fast computations using three‐dimensional high‐resolution models involving more than 1.3 billion degrees of freedom. We propose a new physics‐based approach explaining spontaneous stress drops in a deforming medium.

The interplay between compaction‐driven fluid flow and plastic yielding within porous media is investigated through numerical modeling. We establish a framework for understanding the dynamics of fluid flow in deforming porous materials that corresponds to the equations describing solitary porosity wave propagation. A concise derivation of the coupled fluid flow and poro‐viscoelastoplastic matrix behavior is presented, revealing a connection to Biot's equations of poroelasticity and Gassmann's theory in the elastic limit. Our findings demonstrate that fluid overpressure resulting from channelized fluid flow initiates the formation of new shear zones. Through three‐dimensional simulations, we observe that the newly formed shear zones exhibit a parabolic shape. Furthermore, plasticity exerts a significant influence on both the velocity of fluid flow and the shape of fluid channels. Importantly, our study highlights the potential of spontaneous channeling of porous fluids to trigger seismic events by activating both new and pre‐existing faults.

The ramped temperature oxidation (RTO) test is a screening method used to assess the stability of a reservoir for air-injection Enhanced Oil Recovery (EOR) and to evaluate the oxidation behavior of oil samples. It provides valuable kinetic data for specific cases. The RTO test allows for the analysis of various characteristics, such as temperature evolution, peak temperatures, oxygen uptake, carbon dioxide generation, oxidation and combustion front velocity, recovered and burned hydrocarbons, and residual coke. The adaptation of RTO experiments to in situ combustion (ISC) modeling involves validation and history matching based on numerical simulation of RTO tests, using 3D digital models of experimental setup. The objective is to estimate the kinetic parameters for a customized reaction model that accurately represents ISC. Within this research, the RTO test was provided for bitumen samples related to the Samara oil region. A 3D digital model of the RTO test is constructed using CMG STARS, a thermal hydrodynamic simulator. The model is designed with multiple layers and appropriate heating regimes to account for uncertainties in the experimental setup and to validate the numerical model. The insulation of the setup affects radial heat transfers and helps to control the observed temperature levels. The modified traditional reaction model incorporates thermal cracking of Asphaltenes, low-temperature oxidation (LTO) of Asphaltenes and Maltenes, and high-temperature combustion of coke. Additionally, the model incorporates high-temperature combustion of light oil in the vapor phase, which is generated through Asphaltenes cracking and LTO reactions.

We propose a coupled hydro-mechanical-chemical model and its 1D numerical implementation. We demonstrate its application to the model filtration of a multicomponent fluid in deforming and reacting host rocks, considering changes in the densities, phase proportions, and the chemical compositions of the coexisting phases. The presented 1D numerical implementation is illustrated by the example of soapstone formation from serpentinite during the filtration of Н₂О−CО₂ fluid with a low CО₂ concentration coupled with the viscous deformation of the mineral matrix, considering the MgO−SiO₂−Н₂О−CО₂ system. The numerical results show the propagation of a porosity wave by means of a viscous (de)compaction mechanism accompanied by the formation of an elongated zone with higher filtration properties. After the formation of such a channel, the formation and propagation of the reaction fronts occur and are associated with the transformation of the mineral composition of the original rock. During H₂O−CO₂ fluid filtration, starting with 1 wt % dissolved CO₂, carbonization of hydrated serpentinite starts, specifically antigorite transforms to magnesite and talc.