Search

Fractures significantly impact the elastic and hydraulic properties of geological structures, influencing fields like geothermal energy, hydrocarbon exploration, nuclear waste disposal and |$\mathrm{ \mathrm{CO}}_2$| storage. Characterizing these formations is challenging due to the scale disparity between seismic wavelengths and fracture sizes. This study leverages decades of analytical and numerical advancements to evaluate the effective mechanical properties of fractured solids at the mesoscopic scale. A novel numerical method for modelling quasi-static Biot's poroelastic equations using graphics processing units (GPUs) is introduced for simulating hydromechanically coupled systems. Capable of handling up to 133 million voxel elements on a single GPU, this method offers unprecedented spatial resolution to model complex fracture networks. The GPU-accelerated solver, FastBiot_QS, delivers exceptional performance, achieving a computational speedup of approximately 520 times compared to central processing unit-based methods. The solver's accuracy is rigorously validated in 1-D and 3-D setups. Simulations reveal that fracture clustering and spatial distribution significantly affect seismic attenuation and modulus dispersion. Clusters of interconnected fractures lead to higher attenuation at higher frequencies, while sparsely distributed fractures result in higher attenuation at lower frequencies. Simulations with log-normal and uniform distributions present intermediate behaviours between densely clustered and sparsely distributed fractures. The study can improve interpretations of seismic data and hydraulic properties in fractured media. [ABSTRACT FROM AUTHOR]

The Accelerated Pseudo-Transient (APT) method is a matrix-free approach used to solve partial differential equations (PDEs), characterized by its reliance on local operations, which makes it highly suitable for parallelization. With the advent of the memory-wall phenomenon around 2005, where memory access speed overtook floating-point operations as the bottleneck in high-performance computing, the APT method has gained prominence as a powerful tool for tackling various PDEs in geosciences. Recent advancements have demonstrated the APT method's computational efficiency, particularly when applied to quasi-static nonlinear problems using Graphical Processing Units (GPUs). This manuscript presents a comprehensive analysis of the APT method, focusing on its application to quasi-static elastic, viscoelastic, and coupled hydro-mechanical problems, specifically those governed by quasi-static Biot's poroelastic equations, across 1D, 2D, and 3D domains. We systematically investigate the optimal numerical parameters required to achieve rapid convergence, offering valuable insights into the method's applicability and efficiency for a range of physical models. Our findings are validated against analytical solutions, underscoring the robustness and accuracy of the APT method in both homogeneous and heterogeneous media. We explore the influence of boundary conditions, non-linearities, and coupling on the optimal convergence parameters, highlighting the method's adaptability in addressing complex and realistic scenarios. To demonstrate the flexibility of the APT method, we apply it to the nonlinear mechanical problem of strain localization using a poro-elasto-viscoplastic rheological model, achieving extremely high resolutions – 10,0002 voxels in 2D and 5123 voxels in 3D – that, to our knowledge, have not been previously explored for such models. Our study contributes significantly to the field by providing a robust framework for the effective implementation of the APT method in solving challenging geophysical problems. Importantly, the results presented in this paper are fully reproducible, with Matlab, symbolic Maple scripts, and CUDA C codes made available in a permanent repository. [ABSTRACT FROM AUTHOR]

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. Plain Language Summary: In this study, we looked at how fluids move through porous rocks and how they interact with the solid frame of the rocks. The physics was explored in two‐ and three‐dimensions by leveraging the power of high‐performance computing (HPC) based on graphical processing units (GPU). We found that two key processes occur at the same time: fluid flow gets concentrated into channels due to the changing pressure and interaction with the solid material, and it also forms dike‐like structures as it pushes into newly formed shear zones. Importantly, our study highlights the potential of spontaneous channeling of porous fluids to trigger seismic events by activating both new and pre‐existing faults. This research underscores the complex relationship between fluid flow dynamics and geomechanical processes, offering insights into the mechanisms underlying earthquake initiation. Key Points: We present the numerical modeling of fully coupled fluid flow and poro‐viscoelastoplastic matrix flow with decompaction weakeningWe show that fluid overpressure at the tip of the fluid flow channel triggers the development of non‐symmetric shear bandsWe discover that plastic yielding accelerates the fluid flow and modifies the fluid flow pattern [ABSTRACT FROM AUTHOR]

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. Plain Language Summary: Strain localization is the accumulation of strain in narrow regions of rocks and other materials like metals, soils, and concrete, occurring at different scales. The strength of most geomaterials, particularly rocks, is strongly pressure‐dependent, with strength increasing with increasing pressure. We developed efficient numerical algorithms using High‐Performance Computing (HPC) and graphical processing units (GPUs) to model strain localization in 2D and 3D for applications in geodynamics and earthquake physics. Unlike previous models, our method reveals non‐symmetrical patterns by using very small strain increments, highlighting the need for high‐detail modeling. We found that elasto‐plastic and visco‐plastic models show similar strain patterns for relevant materials. Our method also achieves fast, detailed computations with over 1.3 billion variables and offers a new explanation for sudden stress drops in deforming materials. Key Points: We resolve material instability during deformation resulting in a non‐symmetric pattern of strain localizationWe demonstrate the similarity in patterns of strain localization between frictional and time‐dependent plasticity modelsWe achieve fast numerical simulations in high‐resolution model setups in three dimensions involving more than 500 million degrees of freedom [ABSTRACT FROM AUTHOR]

Understanding the underlying mechanisms of seismic attenuation and moduli dispersion in fluid-saturated cracked porous rocks is of great importance for the development of non-invasive methods to characterize the subsurface. Wave-induced fluid flow at the pore scale, so-called squirt flow, is responsible for seismic attenuation and moduli dispersion at sonic and ultra-sonic frequencies and may be relevant at seismic frequencies. The squirt flow associated attenuation is usually quantified using analytical models. However, numerical experiments suggest that the squirt flow related dissipation is sensitive to fine details of the pore geometry, which can only be modelled numerically. Most of the existing numerical studies explore this phenomenon using simplified models, and there is a lack of numerical studies that model the physics in realistic pore geometries with sufficient numerical resolution. As a result, the impact of wave-induced fluid flow on the effective static and time-dependent mechanical characteristics in realistic settings is still poorly understood. I address these issues by developing a numerical method to model the effective mechanical properties of a hydromechanically coupled system at the pore scale suitable for graphical processing units. A numerical evaluation of attenuation and modulus dispersion due to squirt flow in models based on 3-D microtomography images of cracked Carrara marble is presented. It is shown that the local hydraulic conductivity can be quantitatively estimated from the numerically evaluated effective properties. The accuracy of the numerical results is carefully analysed. This study improves the understanding of the underlying mechanisms of attenuation and moduli dispersion in fluid-saturated cracked rocks. The new method can be applied to model squirt flow for entire laboratory samples in the centimetre scale which was not possible a decade ago. [ABSTRACT FROM AUTHOR]

This research paper presents a comprehensive approach to investigating dry granular collapses in three dimensions, by combining analytical, numerical, and experimental methods. The experimental investigation utilised a novel apparatus to study granular collapses in the laboratory. It is demonstrated that a quasistatic understanding of granular collapses can accurately predict the final normalised run-out distances for dynamic granular collapses. Our approach involved establishing a correlation between the angle of repose and the initial aspect ratio of the granular column. It is also shown that the material point method (MPM) is suitable for modelling granular collapses in three dimensions. Our in-house solver was further validated using experimental evidence under an explicit formulation, resulting in good agreement between the numerical and experimental results. The findings demonstrate the effectiveness of our in-house solver for three-dimensional granular collapse modelling. [ABSTRACT FROM AUTHOR]

We propose an explicit GPU-based solver within the material point method (MPM) framework using graphics processing units (GPUs) to resolve elastoplastic problems under two- and three-dimensional configurations (i.e. granular collapses and slumping mechanics). Modern GPU architectures, including Ampere, Turing and Volta, provide a computational framework that is well suited to the locality of the material point method in view of high-performance computing. For intense and non-local computational aspects (i.e. the back-and-forth mapping between the nodes of the background mesh and the material points), we use straightforward atomic operations (the scattering paradigm). We select the generalized interpolation material point method (GIMPM) to resolve the cell-crossing error, which typically arises in the original MPM, because of the C0 continuity of the linear basis function. We validate our GPU-based in-house solver by comparing numerical results for granular collapses with the available experimental data sets. Good agreement is found between the numerical results and experimental results for the free surface and failure surface. We further evaluate the performance of our GPU-based implementation for the three-dimensional elastoplastic slumping mechanics problem. We report (i) a maximum 200-fold performance gain between a CPU- and a single-GPU-based implementation, provided that (ii) the hardware limit (i.e. the peak memory bandwidth) of the device is reached. Furthermore, our multi-GPU implementation can resolve models with nearly a billion material points. We finally showcase an application to slumping mechanics and demonstrate the importance of a three-dimensional configuration coupled with heterogeneous properties to resolve complex material behaviour. [ABSTRACT FROM AUTHOR]

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. Key Points: We present the dimensional analysis of Biot's equationsWe perform three dimensional numerical simulations of poroelastic wave propagationWe propose a multi‐graphical processing units implementation resolving over 1.5 billion grid cells in a few seconds with near ideal parallel efficiency [ABSTRACT FROM AUTHOR]

We present an efficient MATLAB-based implementation of the material point method (MPM) and its most recent variants. MPM has gained popularity over the last decade, especially for problems in solid mechanics in which large deformations are involved, such as cantilever beam problems, granular collapses and even large-scale snow avalanches. Although its numerical accuracy is lower than that of the widely accepted finite element method (FEM), MPM has proven useful for overcoming some of the limitations of FEM, such as excessive mesh distortions. We demonstrate that MATLAB is an efficient high-level language for MPM implementations that solve elasto-dynamic and elasto-plastic problems. We accelerate the MATLAB-based implementation of the MPM method by using the numerical techniques recently developed for FEM optimization in MATLAB. These techniques include vectorization, the use of native MATLAB functions and the maintenance of optimal RAM-to-cache communication, among others. We validate our in-house code with classical MPM benchmarks including (i) the elastic collapse of a column under its own weight; (ii) the elastic cantilever beam problem; and (iii) existing experimental and numerical results, i.e. granular collapses and slumping mechanics respectively. We report an improvement in performance by a factor of 28 for a vectorized code compared with a classical iterative version. The computational performance of the solver is at least 2.8 times greater than those of previously reported MPM implementations in Julia under a similar computational architecture. [ABSTRACT FROM AUTHOR]

Understanding the properties of cracked rocks is of great importance in scenarios involving CO2 geological sequestration, nuclear waste disposal, geothermal energy, and hydrocarbon exploration and production. Developing noninvasive detecting and monitoring methods for such geological formations is crucial. Many studies show that seismic waves exhibit strong dispersion and attenuation across a broad frequency range due to fluid flow at the pore scale known as squirt flow. Nevertheless, how and to what extent squirt flow affects seismic waves is still a matter of investigation. To fully understand its angle- and frequency-dependent behavior for specific geometries, appropriate numerical simulations are needed. We perform a three-dimensional numerical study of the fluid–solid deformation at the pore scale based on coupled Lamé–Navier and Navier–Stokes linear quasistatic equations. We show that seismic wave velocities exhibit strong azimuth-, angle- and frequency-dependent behavior due to squirt flow between interconnected cracks. Furthermore, the overall anisotropy of a medium mainly increases due to squirt flow, but in some specific planes the anisotropy can locally decrease. We analyze the Thomsen-type anisotropic parameters and adopt another scalar parameter which can be used to measure the anisotropy strength of a model with any elastic symmetry. This work significantly clarifies the impact of squirt flow on seismic wave anisotropy in three dimensions and can potentially be used to improve the geophysical monitoring and surveying of fluid-filled cracked porous zones in the subsurface. [ABSTRACT FROM AUTHOR]

We explore the impact of roughness in crack walls on the P wave modulus dispersion and attenuation caused by squirt flow. For that, we numerically simulate oscillatory relaxation tests on models having interconnected cracks with both simple and intricate aperture distributions. Their viscoelastic responses are compared with those of models containing planar cracks but having the same hydraulic aperture as the rough wall cracks. In the absence of contact areas between crack walls, we found that three apertures affect the P wave modulus dispersion and attenuation: the arithmetic mean, minimum aperture, and hydraulic aperture. We show that the arithmetic mean of the crack apertures controls the effective P wave modulus at the low‐ and high‐frequency limits, thus representing the mechanical aperture. The minimum aperture of the cracks tends to dominate the energy dissipation process and, consequently, the characteristic frequency. An increase in the confining pressure is emulated by uniformly reducing the crack apertures, which allows for the occurrence of contact areas. The contact area density and distribution play a dominant role in the stiffness of the model, and in this scenario, the arithmetic mean is not representative of the mechanical aperture. On the other hand, for a low percentage of minimum aperture or in the presence of contact areas, the hydraulic aperture tends to control the characteristic frequency. Analyzing the local energy dissipation, we can more specifically visualize that a different aperture controls the energy dissipation process at each frequency, which means that a frequency‐dependent hydraulic aperture might describe the squirt flow process in cracks with rough walls. Key Points: We solve the quasi‐static linearised Navier‐Stokes equations coupled to elasticity equationsSeismic attenuation due to squirt‐flow is strongly affected by the roughness of the crack wallsThe minimum and the hydraulic apertures significantly affect the energy dissipation process [ABSTRACT FROM AUTHOR]

Understanding the properties of cracked rocks is of great importance in scenarios involving CO2 geological sequestration, nuclear waste disposal, geothermal energy, and hydrocarbon exploration and production. Developing non-invasive detecting and monitoring methods for such geological formations is crucial. Many studies show that seismic waves exhibit strong dispersion and attenuation across a broad frequency range due to fluid flow at the pore scale known as squirt flow. Nevertheless, how and to what extent squirt flow affects seismic waves is still a matter of investigation. To fully understand its angle- and frequency-dependent behavior for specific geometries appropriate numerical simulations are needed. We perform a three-dimensional numerical study of the fluid-solid deformation at the pore scale based on coupled Lame-Navier and Navier-Stokes linear quasistatic equations. We show that seismic wave velocities exhibit strong azimuth-, angle- and frequency-dependent behavior due to squirt flow between interconnected cracks. We show that the overall anisotropy of a medium mainly increases due to squirt flow but in some specific planes the anisotropy can locally decrease. We analyze the Thomsen-type anisotropic parameters and adopt another scalar parameter which can be used to measure the anisotropy strength of a model with any elastic symmetry. This work significantly clarifies the impact of squirt flow on seismic wave anisotropy in three dimensions and can potentially be used to improve the geophysical monitoring and surveying of fluid-filled cracked porous zones in the subsurface. [ABSTRACT FROM AUTHOR]

Biot's theory describes the coupled solid-fluid interaction in a porous medium. In the earth sciences, poroelasticity is essential in many applications, for example, in hydrogeology, in seismic monitoring of CO2 reservoirs, etc. In the present work, we solve Biot's equations for particle velocity and stress in the time domain with a finite difference approach on a staggered grid. Depending on the medium's properties, besides the elastic compressional and shear waves, one may observe the propagating slow wave or the diffusive static slow mode. In such a case, when the diffusive slow mode is present, the time-stepping should be very small in order to reach a stable simulation of Biot's equations. We present a new treatment of this problem based on the so-called pseudo-transient method. The idea is that at each time step, another pseudo iteration causes the slow mode to attenuate quickly. As a result, very fine time-stepping is unnecessary and the time-stepping is controlled by a standard Courant stability condition for the fast P-wave. Furthermore, we accelerate a finite-difference code using an NVIDIA graphics card with the CUDA programming language. [ABSTRACT FROM AUTHOR]