Seismic waves may exhibit significant dispersion and attenuation in reservoir rocks due to pore-scale fluid flow. Fluid flow at the microscopic scale is referred to as squirt flow and occurs in very compliant pores, such as grain contacts or microcracks, that are connected to other stiffer pores. We have performed 3D numerical simulations of squirt flow using a finite-element approach. Our 3D numerical models consist of a pore space embedded into a solid grain material. The pore space is represented by a flat cylinder (a compliant crack) whose edge is connected with a torus (a stiff pore). Grains are described as a linear isotropic elastic material, whereas the fluid phase is described by the quasistatic linearized compressible Navier-Stokes momentum equation. We obtain the frequency-dependent effective stiffness of a porous medium and calculate dispersion and attenuation due to fluid flow from a compliant crack to a stiff pore. We compare our numerical results against a published analytical solution for squirt flow and analyze the effects of its assumptions. Previous interpretation of the squirt flow phenomenon based mainly on analytical solutions is verified, and some new physical effects are identified. The numerical and analytical solutions agree only for the simplest model in which the edge of the crack is subjected to zero fluid pressure boundary condition while the stiff pore is absent. For the more realistic model that includes the stiff pore, significant discrepancies are observed. We identify two important aspects that need improvement in the analytical solution: the calculation of the frame stiffness moduli and the frequency dependence of attenuation and dispersion at intermediate frequencies.

Seismic wave propagation in porous rocks that are saturated with a liquid exhibits significant dispersion and attenuation due to fluid flow at the pore scale, so-called squirt flow. This phenomenon takes place in compliant flat pores such as microcracks and grain contacts that are connected to stiffer isometric pores. Accurate quantitative description is crucial for inverting rock and fluid properties from seismic attributes such as attenuation. Up to now, many analytical models for squirt flow were proposed based on simplified geometries of the pore space. These models were either not compared with a numerical solution or showed poor accuracy. We present a new analytical model for squirt flow which is validated against a three-dimensional numerical solution for a simple pore geometry that has been classically used to explain squirt flow; that is why we refer to it as classical geometry. The pore space is represented by a flat cylindrical (penny-shaped) pore whose curved edge is fully connected to a toroidal (stiff) pore. Compared with correct numerical solutions, our analytical model provides very accurate predictions for the attenuation and dispersion across the whole frequency range. This includes correct low- and high-frequency limits of the stiffness modulus, the characteristic frequency, and the shape of the dispersion and attenuation curves. In a companion paper (Part 2), we extend our analytical model to more complex pore geometries. We provide as supplementary material Matlab and symbolic Maple routines to reproduce our main results.

Summary Two-phase flow equations that couple solid deformation and fluid migration have opened new research trends in geodynamical simulations and modelling of subsurface engineering. Physical nonlinearity of fluid-rock systems and strong coupling between flow and deformation in such equations lead to interesting predictions such as spontaneous formation of focused fluid flow in ductile/plastic rocks. However, numerical implementation of two-phase flow equations and their application to realistic geological environments with complex geometries and multiple stratigraphic layers is challenging. This study documents an efficient pseudo-transient solver for two-phase flow equations and describes the numerical theory and physical rationale. We provide a simple explanation for all steps involved in the development of a pseudo-transient numerical scheme for various types of equations. Two different constitutive models are used in our formulations: a bilinear viscous model with decompaction weakening and a viscoplastic model that allows decompaction weakening at positive effective pressures. The resulting numerical models are used to study fluid leakage from high porosity reservoirs into less porous overlying rocks. The interplay between time-dependent rock deformation and the buoyancy of ascending fluids leads to the formation of localized channels. The role of material parameters, reservoir topology, geological heterogeneity and porosity is investigated. Our results show that material parameters control the propagation speed of channels while the geometry of the reservoir controls their locations. Geological layers present in the overburden do not stop the propagation of the localized channels but rather modify their width, permeability, and growth speed.

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 equations • Seismic attenuation due to squirt-flow is strongly affected by the roughness of the crack walls • The minimum and the hydraulic apertures significantly affect the energy dissipation process presented a comparison between numerical results and an analytical model for squirt flow. In general, accepted analytical models should reproduce the equations of Gassmann (1951) in the low-frequency limit (Chapman et al., 2002). The reason is that at the relaxed state for undrained boundary conditions (low-frequency limit), the time of a half period of a passing wave allows for fluid pressure to equilibrate through FPD. At the unrelaxed state (high-frequency limit), the fluid pressure has no time to equilibrate during a half period of a passing wave and the elastic properties of the saturated material are predicted by the formulation of Mavko and Jizba (1991), which assumes that no FPD occurs during the passage of the wave. At intermediate frequencies, FPD occurs inside the cracks during the passage of the wave and part of its energy is dissipated. Nevertheless, all analytical solutions assume smooth walls for the cracks despite the fact that crack walls in rocks have been observed to present complex profiles including wall roughnes...

Summary The efficient and accurate numerical modeling 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.

Abstract. 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.

Majority of the most powerful supercomputers on the world host hardware accelerators to sustain calculations at the petascale level and beyond. Graphical processing units (GPUs) are amongst widely employed hardware accelerators, initiating a revolution in high-performance computing (HPC) in the last decade. The three-dimensional calculations targeting billions of grid cells -technically impossible resolutions decades ago -became reality. This major breakthrough in HPC and supercomputing comes however with the cost of developing and reengineering scientific codes to efficiently utilize the available computing power. Increasing the low-level parallelism is the key. In Earth sciences, HPC and GPU-accelerated applications target in particular forward and inverse seismic modeling and geodynamics -fields where high spatial and temporal resolutions as well as large spatial domains are required. We here develop a multi-GPU implementation for applications in seismic modeling in porous media.Understanding seismic wave propagation in fluid-saturated porous media enables more accurate interpretation of seismic signals in Earth sciences. The two phase medium is represented by an elastic solid matrix (skeleton) saturated with a compressible viscous fluid. The dynamic response of such an isotropic two phase medium results in two longitudinal waves and one shear wave, as predicted by Frenkel (1944) (see also Pride and Garambois (2005)). The wave of the first kind (fast wave) is a true longitudinal wave where the solid matrix motion and the fluid velocity vector fields are in-phase. The wave of the second kind (slow wave) is a highly attenuated wave where the solid matrix motion and the fluid velocity vector fields are

Abstract. 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.