https://orcid.org/0000-0002-6529-9256
SUMMARY
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.
1 INTRODUCTION
Fractures play a crucial role in earth sciences and civil engineering due to their profound impact on the elastic and hydraulic properties of geological structures. Characterizing fractured formations, particularly in terms of their hydraulic behaviour, is crucial for many practical applications, including geothermal energy production, hydrocarbon exploration, nuclear waste disposal and |$\mathrm{CO}_2$| storage. Direct correlation of measured seismic attributes to individual fractures proves challenging due to the disparity in scale between seismic wavelengths, which are typically much larger, and the sizes of the fractures. Nonetheless, it is feasible to gather information about fractured rocks within the context of an effective medium theory.
Over the past several decades, a variety of analytical and numerical approaches have been created to assess the effective mechanical properties of cracked solids. For the purposes of this study, these properties are classified into three distinct groups: (i) elasto-static, (ii) time- or frequency-dependent at the pore scale, and (iii) time- or frequency-dependent at the mesoscopic scale, utilizing Biot’s poroelastic equations. Time- or frequency-dependent properties are typically associated with (equivalent) viscoelastic behaviour. The focus of this study is specifically on the third group.
Elasto-static properties generally refer to the assessment of effective (static) elastic properties using different numerical approaches (Garboczi & Day 1995; Garboczi 1998; Saenger et al. 2004, 2016; Grechka & Kachanov 2006; Andrä et al. 2013; Kabel et al. 2016; Saxena et al. 2019) and others. Numerous textbooks, including Kachanov & Sevostianov (2018), offer detailed explanations of effective medium theories that are used to analytically determine the effective elastic properties of cracked solids.
At the pore scale, the effective viscoelastic properties of cracked solids can be determined through direct modelling that integrates coupled elasticity and Navier–Stokes equations, as demonstrated in studies by Quintal et al. (2016), Alkhimenkov et al. (2020a) and Lissa et al. (2020). Research efforts that aim to characterize cracked solids have examined both simplistic cracked configurations (Alkhimenkov et al. 2020b) and more complex, realistic cracked frameworks (Alkhimenkov 2024a).
At the mesoscopic scale, the effective viscoelastic properties of fractured solids are examined within the framework of poroelasticity. Studies utilizing Biot’s poroelastic equations to explore seismic attenuation and stiffness moduli dispersion include works by Masson & Pride (2007), Rubino et al. (2009), Wenzlau et al. (2010) and Quintal et al. (2011). Subsequent numerical investigations into the properties of fractured rocks, employing poroelastic models of interconnected fractures, have been conducted by Rubino et al. (2013), Quintal et al. (2014), Hunziker et al. (2018), Lissa et al. (2019, 2021) and Caspari et al. (2019). Different aspects of wave anelasticity have been studied by Ba et al. (2016) for squirt flow, by Ba et al. (2017) for patchy saturation and fabric heterogeneity, and by Zhang et al. (2021) for infinituple-porosity media.
At present, the performance of numerical computations is primarily limited by memory access speed rather than by the speed of floating-point operations. The advent of multicore technology has substantially advanced our modelling capabilities. Ever since the identification of the “memory wall” phenomenon in the early 2000s (Wulf & McKee 1995), there has been a steady rise in the disparity between computation speeds and memory access speeds. This shift has moved the focus from compute-bound to memory-bound algorithms, marking a significant change in computational strategy. The memory bandwidth of graphics processing units (GPUs) can surpass that of central processing units (CPUs) by several magnitudes. Applications that are specifically designed for GPUs can handle 3-D problems with resolutions that involve billions of grid cells (Alkhimenkov et al. 2021; Räss et al. 2022). In the realm of geosciences, GPU-accelerated algorithms have effectively been used to simulate compaction-driven fluid flow (Omlin et al. 2017, 2018; Räss et al. 2019), thermomechanical deformation of ice (Räss et al. 2020), and wave propagation utilizing elastodynamic Biot’s poroelastic equations (Alkhimenkov et al. 2021). The use of matrix-free pseudo-transient iterative methods for solving quasi-static equations has enabled modelling at exceptionally high-spatial resolutions (Räss et al. 2022). To the best of my knowledge, no studies have yet modelled the quasi-static Biot’s poroelastic equations on GPUs to characterize stochastic fracture networks in 3-D.
In this work, a novel numerical approach to model the effective viscoelastic properties of a hydromechanically coupled system at the mesoscopic scale is introduced, utilizing the quasi-static Biot’s poroelastic equations on GPUs. At this scale, the effective mechanical properties of a fractured porous system are effectively captured through equivalent viscoelastic properties. The implementation of this model is conducted in CUDA C (Compute Unified Device Architecture), optimized for Nvidia GPU devices. The design of this solver is founded on three principal concepts: concise numerical implementation, high-numerical resolution, and exceptional computational efficiency. The code is both concise and optimized to exploit the parallel capabilities of GPUs, achieving significant performance enhancements. It supports high-spatial resolutions with up to 133 million voxel elements using a single GPU, allowing for the detailed representation of complex geometries of the fractured formations. The computational efficiency of the solver ensures that even large domains involving millions of voxel elements can be processed in seconds or minutes. This method is applied to investigate attenuation and modulus dispersion caused by fluid pressure diffusion within a 3-D stochastic fracture networks at the mesoscopic scale, showcasing its practical application in realistic scenarios.
This article introduces several significant contributions to the field:
A novel numerical approach has been introduced, employing quasi-static Biot’s poroelastic equations for the evaluation of viscoelastic properties, specifically optimized for GPUs.
The new method has been rigorously validated against analytical solutions and 3-D finite-element (FE) methods, confirming its precision and effectiveness.
An extensive numerical analysis of attenuation and modulus dispersion has been conducted, examining the impacts of fluid pressure diffusion within 3-D stochastic fracture networks.
Simulations reveal that fracture clustering significantly affect seismic attenuation and modulus dispersion.
The application potential of this methodology for modelling viscoelastic properties of large reservoir models at the meter scale has been demonstrated, underscoring its practical utility.
2 MATHEMATICAL FORMULATION: QUASI-STATIC BIOT’S POROELASTIC EQUATIONS
In this study, we utilize the quasi-static Biot’s equations to model coupled porous-media flow and deformation, assuming a poroelastic behaviour of the continuum. A comprehensive list of symbols used in the analysis is provided in Table 1.
List of symbols.
| Symbol | Meaning |
|---|---|
| |$\sigma ^s_{ij}, \sigma ^f_{ij}$| | Solid and fluid stress tensor |
| |$\bar{\sigma }_{ij}$| | |$=(1-\phi )\sigma ^s_{ij} + \phi \sigma ^f_{ij}$|, total stress tensor |
| |$p_\mathrm{s}, p_\mathrm{f}$| | Solid and fluid pressure |
| |$\bar{p}$| | |$=(1-\phi )p_\mathrm{s} + \phi p_\mathrm{f}$|, total pressure |
| |$\tau ^s_{ij}, \tau ^f_{ij}$| | Solid and fluid deviatoric stress tensor |
| |$\bar{\tau }_{ij}$| | |$=(1-\phi )\tau ^s_{ij} + \phi \tau ^f_{ij}$|, total deviatoric stress tensor |
| |$v^s, v^f$| | Solid and fluid velocity |
| |$q^D_i$| | |$=\phi (v^f - v^s)$|, Darcy’s flux |
| |$\rho ^s$| | Solid density |
| |$K_\mathrm{s},K_\mathrm{f}$| | Elastic solid and fluid bulk modulus |
| |$G_\mathrm{s}, G_\mathrm{d} = G_\mathrm{u}$| | Elastic solid, drained and undrained shear modulus |
| |$K_{d},K_{u}$| | Elastic drained and undrained bulk modulus |
| |$\eta _\mathrm{f}$| | Fluid shear viscosity |
| k | Permeability |
| |$\phi$| | Porosity |
| |$\alpha$| | Biot–Willis coefficient |
| B | Skempton coefficient |
| |$\delta _{ij}$| | Kronecker delta |
| |$i,j,k=\overline{1..3}$| | Indexes |
| Symbol | Meaning |
|---|---|
| |$\sigma ^s_{ij}, \sigma ^f_{ij}$| | Solid and fluid stress tensor |
| |$\bar{\sigma }_{ij}$| | |$=(1-\phi )\sigma ^s_{ij} + \phi \sigma ^f_{ij}$|, total stress tensor |
| |$p_\mathrm{s}, p_\mathrm{f}$| | Solid and fluid pressure |
| |$\bar{p}$| | |$=(1-\phi )p_\mathrm{s} + \phi p_\mathrm{f}$|, total pressure |
| |$\tau ^s_{ij}, \tau ^f_{ij}$| | Solid and fluid deviatoric stress tensor |
| |$\bar{\tau }_{ij}$| | |$=(1-\phi )\tau ^s_{ij} + \phi \tau ^f_{ij}$|, total deviatoric stress tensor |
| |$v^s, v^f$| | Solid and fluid velocity |
| |$q^D_i$| | |$=\phi (v^f - v^s)$|, Darcy’s flux |
| |$\rho ^s$| | Solid density |
| |$K_\mathrm{s},K_\mathrm{f}$| | Elastic solid and fluid bulk modulus |
| |$G_\mathrm{s}, G_\mathrm{d} = G_\mathrm{u}$| | Elastic solid, drained and undrained shear modulus |
| |$K_{d},K_{u}$| | Elastic drained and undrained bulk modulus |
| |$\eta _\mathrm{f}$| | Fluid shear viscosity |
| k | Permeability |
| |$\phi$| | Porosity |
| |$\alpha$| | Biot–Willis coefficient |
| B | Skempton coefficient |
| |$\delta _{ij}$| | Kronecker delta |
| |$i,j,k=\overline{1..3}$| | Indexes |
List of symbols.
| Symbol | Meaning |
|---|---|
| |$\sigma ^s_{ij}, \sigma ^f_{ij}$| | Solid and fluid stress tensor |
| |$\bar{\sigma }_{ij}$| | |$=(1-\phi )\sigma ^s_{ij} + \phi \sigma ^f_{ij}$|, total stress tensor |
| |$p_\mathrm{s}, p_\mathrm{f}$| | Solid and fluid pressure |
| |$\bar{p}$| | |$=(1-\phi )p_\mathrm{s} + \phi p_\mathrm{f}$|, total pressure |
| |$\tau ^s_{ij}, \tau ^f_{ij}$| | Solid and fluid deviatoric stress tensor |
| |$\bar{\tau }_{ij}$| | |$=(1-\phi )\tau ^s_{ij} + \phi \tau ^f_{ij}$|, total deviatoric stress tensor |
| |$v^s, v^f$| | Solid and fluid velocity |
| |$q^D_i$| | |$=\phi (v^f - v^s)$|, Darcy’s flux |
| |$\rho ^s$| | Solid density |
| |$K_\mathrm{s},K_\mathrm{f}$| | Elastic solid and fluid bulk modulus |
| |$G_\mathrm{s}, G_\mathrm{d} = G_\mathrm{u}$| | Elastic solid, drained and undrained shear modulus |
| |$K_{d},K_{u}$| | Elastic drained and undrained bulk modulus |
| |$\eta _\mathrm{f}$| | Fluid shear viscosity |
| k | Permeability |
| |$\phi$| | Porosity |
| |$\alpha$| | Biot–Willis coefficient |
| B | Skempton coefficient |
| |$\delta _{ij}$| | Kronecker delta |
| |$i,j,k=\overline{1..3}$| | Indexes |
| Symbol | Meaning |
|---|---|
| |$\sigma ^s_{ij}, \sigma ^f_{ij}$| | Solid and fluid stress tensor |
| |$\bar{\sigma }_{ij}$| | |$=(1-\phi )\sigma ^s_{ij} + \phi \sigma ^f_{ij}$|, total stress tensor |
| |$p_\mathrm{s}, p_\mathrm{f}$| | Solid and fluid pressure |
| |$\bar{p}$| | |$=(1-\phi )p_\mathrm{s} + \phi p_\mathrm{f}$|, total pressure |
| |$\tau ^s_{ij}, \tau ^f_{ij}$| | Solid and fluid deviatoric stress tensor |
| |$\bar{\tau }_{ij}$| | |$=(1-\phi )\tau ^s_{ij} + \phi \tau ^f_{ij}$|, total deviatoric stress tensor |
| |$v^s, v^f$| | Solid and fluid velocity |
| |$q^D_i$| | |$=\phi (v^f - v^s)$|, Darcy’s flux |
| |$\rho ^s$| | Solid density |
| |$K_\mathrm{s},K_\mathrm{f}$| | Elastic solid and fluid bulk modulus |
| |$G_\mathrm{s}, G_\mathrm{d} = G_\mathrm{u}$| | Elastic solid, drained and undrained shear modulus |
| |$K_{d},K_{u}$| | Elastic drained and undrained bulk modulus |
| |$\eta _\mathrm{f}$| | Fluid shear viscosity |
| k | Permeability |
| |$\phi$| | Porosity |
| |$\alpha$| | Biot–Willis coefficient |
| B | Skempton coefficient |
| |$\delta _{ij}$| | Kronecker delta |
| |$i,j,k=\overline{1..3}$| | Indexes |
The stress–strain relation assuming small strains and the conservation of mass for porous media saturated with a fluid can be written as (Biot 1962; Wang 2017):
where the total pressure |$\bar{p}$| is defined as positive in compression (|$\bar{p}= -\bar{\sigma }_{kk} /3 \equiv -\mathrm{tr}(\bar{\sigma }_{ij}) /3$|), |$\mathrm{tr}$| denotes the trace of a tensor. The deviatoric part of the total stress tensor is
The material constants in the expression (1) are the the Biot–Willis coefficient, |$\alpha$|,
and the Skempton coefficient, B,
Other useful parameters include the undrained bulk modulus, |$K_\mathrm{ u}$|,
and the fluid storage coefficient, M,
Expression (5) is also known as the Gassmann equation (Gassmann 1951; Alkhimenkov 2023). The quasi-static Biot poroelastic equations represent the inertialess formulation of the full elasto-dynamic Biot equations. The conservation of linear momentum (the equilibrium equation) and Darcy’s law are
3 NUMERICAL SOLUTION STRATEGY
3.1 Discretization
Consider a volume V within 3-D Euclidean space |$E^3$|, bounded by a regular surface |$\partial V$|. The volume is defined as a cuboid with dimensions |$V = (-L_x/2,L_x/2)\times (-L_y/2,L_y/2)\times (-L_z/2,L_z/2)$|, |$L_x=L_y=L_z=L^{*}=1$|. Here, |$L^{*}$| represents the non-dimensional length of the computational domain. The spatial and temporal discretization of this volume is achieved using a regular staggered numerical grid, as outlined by Virieux (1986), which facilitates a conservative formulation comparable to a type of finite-volume method described by Dormy & Tarantola (1995).
To solve the quasi-static Biot equations numerically, as outlined in eqs (1), (2) and (7), it is necessary to invert the coefficient matrices specified in (1). This process yields the following expressions:
and
Note that the coefficient matrices in eq. (8) are symmetric by analogy with the system of eq. (1).
3.2 Accelerated pseudo-transient method
Rather than employing a direct solver, we opt for the iterative, matrix-free pseudo-transient method (Frankel 1950; Räss et al. 2022). This method is particularly suited for solving quasi-static equations with high resolution, such as |$N = 511^3 = 133\,000\,000$| voxels, which far exceed the computational capacity of modern direct solvers like MUMPS and PARDISO. These direct solvers are typically restricted to finite element models with up to approximately 5 million elements.
The accelerated pseudo-transient method iteratively solves a modified partial differential equation (PDE) by introducing inertial and relaxation terms, which vanish as the solution converges to the original PDE. Rather than solving inertialess equations directly, the method solves equations with inertia by attenuating dynamic fields, achieving the desired quasi-static solution as these fields diminish to a precise level. Essentially, the quasi-static solution emerges as an attractor within the hyperbolic framework (equations featuring both inertia and attenuation). This technique is also referred to as the relaxation method (Frankel 1950). This matrix-free method relies solely on local operations, eliminating the need for global memory access and enabling efficient parallelization on modern computing architectures, particularly GPUs.
To solve the quasi-static Biot poroelastic equations, the process involves two primary steps. Initially, inertial terms are incorporated into eq. (10); subsequently, inertial terms are responsible for wave propagation in pseudo physical time and space (i.e. hyperbolic) and viscous terms (treated as a Maxwell rheology integrated into eqs 8 and 9) are the physical quantities. Following these modifications, the quasi-static Biot’s poroelastic eqs (8)–(10) are reformulated to include pseudo-time |$\widetilde{t}$|. The pseudo-transient version of the equation for the total pressure, |$\bar{p}$|, is (see expression 8):
where |$\hat{\bar{p}}$| is the total pressure at the previous physical time step, which remains constant across pseudo-time iterations, |$\Delta t$| is the physical time step, and |$\widetilde{K}_1$| is a numerical parameter. The pseudo-transient version of the equation for the fluid pressure, |$p_\mathrm{f}$|, is (see expression 8):
where |$\hat{p}_\mathrm{ f}$| is the fluid pressure field at the previous physical time step. For the deviatoric stress the corresponding equation is
where |$\widetilde{G}$| is a numerical parameter. The conservation of linear momentum equations and the conservation of mass for the fluid equations become the following:
We have omitted the added mass coefficients |$\rho _\mathrm{f}$| in the expression (14) because they do not influence the outcome. This is due to the fact that the inertial terms |${\partial v^s_i}/{\partial \widetilde{t}}$| diminish to a specific level of precision, such as |$10^{-8}$|.
In eqs (14)–(13), the quantities |$\partial \widetilde{t}$|, |$\widetilde{\rho }_t$| and |$\widetilde{\rho }_a$| are numerical parameters. More precisely, the following combinations of the numerical parameters are needed |$\widetilde{K }_1 \Delta \widetilde{t}$|, |$\widetilde{G } \Delta \widetilde{t}$|, |$\Delta \widetilde{t} /\widetilde{\rho }_t$| and |$\Delta \widetilde{t} /\widetilde{\rho }_a$|. The optimal values for these parameters can be established through numerical methods or, in certain instances, through analytical approaches (Räss et al. 2022). Such an analysis of Biot’s poroelastc equations has been performed by Alkhimenkov & Podladchikov (2024). The following expression is employed for the first numerical combination:
where |$\widetilde{V}^{\mathrm{ LF}}_\mathrm{ p} \, \Delta \widetilde{t} = \widetilde{C} \Delta x$|, |$\widetilde{C} \le 1$| and the Strouhal number, |${\mathrm{St}}$|, is expressed as
The second numerical combination is
where |$r = {K}_\mathrm{ u}/ {G}_\mathrm{ u}$|. Calculation of |$\widetilde{K }_1 \Delta \widetilde{t}$| is also straightforward: |$\widetilde{K }_1 \, \Delta \widetilde{t} = r \, \widetilde{G } \, \Delta \widetilde{t}$|. The last combination, |$\Delta \widetilde{t} /\widetilde{\rho }_a$|, is
As a result, the system of eqs (11)–(14) can be solved. In 1-D, the optimal value of |$\mathrm{St}$| is (Alkhimenkov & Podladchikov 2024)
In 2-D and 3-D, the values are |${\mathrm{St}}_\text{opt} = 2\, \pi \, \sqrt{2}$| and |${\mathrm{St}}_\text{opt} = 2\, \pi \, \sqrt{3}$|, respectively (Alkhimenkov & Podladchikov 2024).
3.3 Non-dimensionalization
Non-dimensionalization is an essential technique for the analysis of partial differential equations. This method effectively highlights crucial characteristic scales, streamlines the equations and minimizes the number of input parameters needed for calculations. Dimensional analysis for the elastodynamic Biot’s equations has previously been conducted by Alkhimenkov et al. (2021). Building on this foundation, we adopt a similar approach for the non-dimensionalization of the quasi-static Biot’s equations. We define characteristic scales, specifically |$L_x^{*}$| for length and |$T^{*}$| for time. Subsequent substitutions are as follows:
where the superscript |$\,\, \breve{} \,\,$| denotes the dimensionless quantities.
The resulting dimensionless quasi-static Biot’s equations are:
where |$K_\mathrm{ G}= K_\mathrm{ d}/G_\mathrm{ d}$| and the coefficient |$I_2$| is
Now we set |$I_2 = 1$|, |$L_x=1$|, |$K_\mathrm{ d}=1$|, then |$T^{*} = \frac{\eta _\mathrm{ f} }{k }$|. Since |$T^{*} = 1/( \omega ^{*}) = 1/(2\, \pi \, f)$|, where f is the dimensionless frequency, the behaviour of the system is fully controlled by the fluid mobility |$D=\frac{k }{\eta _\mathrm{f} }$| for a given frequency |$f_i$|. Or, in other words, one can vary frequency |$f_i$| for a specific |$D=\frac{k }{\eta _\mathrm{f} }$| and obtain the same result. In the present analysis, we fix |$I_2$| and plot solutions as a function of a dimensionless frequency |$f=f_i$| (e.g. assuming that |$T^{*} =1$|, |$K_\mathrm{d}=1$| and |$L_x=1$|, which gives |$I_2 = \frac{\eta _\mathrm{f} }{k }$|). Below, we will deal only with the dimensionless equations, and, for simplicity, we will omit the overscore tilde symbol |$\, \breve{ } \,$| in the rest of the paper.
3.4 Boundary conditions
Direct tests are conducted to evaluate the dissipation and effective moduli (components of the stiffness tensor) as a function of the dimensionless frequency |$f_i$| (Alkhimenkov et al. 2020a). The continuum’s effective poroelastic behaviour is modelled using a viscoelastic complex-valued equivalent modulus at each discrete frequency |$f_i$|. This method yields an accurate representation of the quasi-static Biot’s poroelastic equations. This study focuses on analysing the dispersion and attenuation of the |$H(f_i)$| component, which governs the primary, or P wave, velocity dispersion and attenuation propagating in the |$z$|-direction.
The derived effective viscoelastic tensor is expressed in Voigt notation (|$C_{ijkl} = C_{mn}$|, where |$i,j, k, l = \overline{1,3}$| and |$m,n = \overline{1,6}$|). To compute the modulus |$C_{33}(f_i) \equiv H(f_i)$|, a velocity boundary condition |$v=1$| is applied perpendicularly at the top external wall in the |$z$|-direction, with a symmetric condition |$v=-1$| at the bottom external wall. All other domain boundaries are treated with standard reflecting boundary conditions. Stress and strain derivatives for a specific frequency |$f_i$| are spatially averaged over the numerical domain |$V$|. The effective modulus, |$H(f_i)$|, is then computed as:
where |$\langle \cdot \rangle$| indicates volume averaging across the numerical domain |$V$|,
and
where |$x \in V$| is a spatial coordinate. The inverse quality factor |$\frac{1}{Q(f_i)}$|, serving as a measure of dissipation (attenuation) for each component of the stiffness tensor, is defined by:
where |$H$| represents the complex-valued P-wave modulus, with |$\text{Im}$| and |$\text{Re}$| denoting the imaginary and real parts, respectively. Notably, our GPU-based numerical solver operates in the time domain using real quantities, necessitating an alternative expression for |$\frac{1}{Q(f_i)}$|. In anisotropic media, different components of the stiffness tensor govern the velocity and attenuation of wave propagation in varying directions. To quantify the attenuation associated with the |$C_{33}(f_i) \equiv H(f_i)$| component, we adopt the following expression for the inverse quality factor (modified after Aki & Richards (2002))
where |$\langle \sigma _{zz}^{*}(\omega _i) \rangle$| corresponds to the stress calculated at the second physical time step. The initial calculation for a given |$\Delta t_i$| results in |$\langle \sigma _{zz}(\omega _i) \rangle$|, while the subsequent calculation at the same |$\Delta t_i$| produces |$\langle \sigma _{zz}^{*}(\omega _i) \rangle$|. For further analytical verification, we apply the approximate Kramers–Kronig relations to estimate |$\frac{1}{Q}$| as follows:
This equation is routinely utilized throughout the manuscript to validate the numerical solver’s computations.
3.5 Implementation using GPUs
The numerical method outlined herein is structured around local computations, making it inherently suitable for parallelization on GPUs. The algorithm involves two nested for loops: the outer loop iterates over the discrete physical time steps |$t_i$|, while the inner loop addresses the pseudo-time |$\widetilde{t}_i$|, necessary to achieve the quasi-static equilibrium at each time step |$t_i$|. For determining the dispersion curve (effective P-wave modulus) as a function of dimensionless frequency f, only a single iteration over the physical time steps |$t_i$| is required. In contrast, the calculation of dissipation necessitates iterating over two physical time steps |$t_i$|. This approach allows for the quantification of dissipation through the difference in stress amplitudes observed between the consecutive physical time steps.
The resulting code is named FastBiot_QS. Initial code development occurred on a laptop equipped with a 13th Generation Intel Core i9-13900HX CPU (64GB RAM) and an NVIDIA GeForce RTX 4090 (16GB) laptop GPU. Subsequent large-scale 3-D simulations were executed on a node comparable to an NVIDIA DGX-1, featuring four NVIDIA Ampere A100 (80GB) GPUs and an AMD EPYC 7742 (512GB RAM) Server Processor. These simulations utilized a single GPU setup, employing either an RTX 4090 or an A100 GPU.
4 VALIDATION
4.1 Material parameters and non-dimensional framework
All simulations in this study are performed in a non-dimensional framework, a standard approach in high-performance computing to ensure robustness and accuracy of the solver. This framework enables universal applicability of the results, allowing them to be scaled to match different rock types by adjusting the scaling factors.
For instance:
To represent a porous sandstone with bulk and shear moduli |$K_{ {\mathrm{ d}}} = 20 \, \text{GPa}$| and |$G_\mathrm{d} = 20 \, \text{GPa}$|, the dimensionless |$K_\mathrm{d}$| and |$G_\mathrm{d}$| are scaled by |$20 \times 10^9$|.
Similarly, for a stiffer limestone with |$K_{ {\mathrm{ d}}} = 40 \, \text{GPa}$| and |$G = 40 \, \text{GPa}$|, the dimensionless |$K_\mathrm{d}$| and |$G_\mathrm{d}$| are scaled by |$40 \times 10^9$|.
Other parameters are scaled according to the chosen dimensionless timescale |$T^{*}$| and length scale |$L_x^{*} = 1$|, following the non-dimensional analysis (see Section 3.3 for details). This approach ensures that the solver’s results can be seamlessly transferred to dimensional quantities, making them applicable to a wide range of geological materials.
4.2 Validation of one-dimensional code implementation
4.2.1 One-dimensional layered media
The validation of the 1-D code implementation, FastBiot_QS, was rigorously tested against an analytical solution for layered poroelastic media (White et al. 1975; Carcione & Picotti 2006; Quintal et al. 2009). This comparison serves as a robust test of the coupled physics encapsulated within Biot’s poroelastic equations. The material parameters used in the simulation are given in Table 2. The parameters were non-dimensionalized by normalizing with |$G_\mathrm{d}$|, selecting a corresponding dimensionless timescale |$T^{*}$|, and assuming |$L_x^{*}=1$| (see Section 3.3 for details). Fig. 1 illustrates the comparison of the numerical results from the 1-D version of FastBiot_QS with the analytical solutions across a range of dimensionless frequencies f. The agreement between the numerical and analytical results substantiates the accuracy of the 1-D numerical implementation (FastBiot_QS).
Numerical results calculated by using the present GPU solver FastBiot_QS compared against the exact analytical solutions for layered media: (a) the real part of the |$H(f_i)$| component (the P-wave modulus) and (b) corresponding dimensionless attenuation |$1/Q(f_i)$|.
Properties of the model used in the numerical simulations: 1-D layered media
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | |$37\, \textrm {GPa}$| | |$37\, \textrm {GPa}$| |
| |$K_\mathrm{d}$| | |$4.8\, \textrm {GPa}$| | |$4.8\, \textrm {GPa}$| |
| |$K_\mathrm{f}$| | |$0.012 \, \textrm {GPa}$| | |$2.25\, \textrm {GPa}$| |
| |$G_\mathrm{d}$| | |$5.7\, \textrm {GPa}$| | |$5.7\, \textrm {GPa}$| |
| |$\eta _\mathrm{f}$| | |$0.00015\, \textrm {Pa} \cdot \textrm {s}$| | |$0.003\, \textrm {Pa} \cdot \textrm {s}$| |
| k | |$1 \times 10^{-12}\, \textrm {m}^2$| | |$1 \times 10^{-12}\, \textrm {m}^2$| |
| |$\phi$| | 0.3 | 0.3 |
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | |$37\, \textrm {GPa}$| | |$37\, \textrm {GPa}$| |
| |$K_\mathrm{d}$| | |$4.8\, \textrm {GPa}$| | |$4.8\, \textrm {GPa}$| |
| |$K_\mathrm{f}$| | |$0.012 \, \textrm {GPa}$| | |$2.25\, \textrm {GPa}$| |
| |$G_\mathrm{d}$| | |$5.7\, \textrm {GPa}$| | |$5.7\, \textrm {GPa}$| |
| |$\eta _\mathrm{f}$| | |$0.00015\, \textrm {Pa} \cdot \textrm {s}$| | |$0.003\, \textrm {Pa} \cdot \textrm {s}$| |
| k | |$1 \times 10^{-12}\, \textrm {m}^2$| | |$1 \times 10^{-12}\, \textrm {m}^2$| |
| |$\phi$| | 0.3 | 0.3 |
Properties of the model used in the numerical simulations: 1-D layered media
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | |$37\, \textrm {GPa}$| | |$37\, \textrm {GPa}$| |
| |$K_\mathrm{d}$| | |$4.8\, \textrm {GPa}$| | |$4.8\, \textrm {GPa}$| |
| |$K_\mathrm{f}$| | |$0.012 \, \textrm {GPa}$| | |$2.25\, \textrm {GPa}$| |
| |$G_\mathrm{d}$| | |$5.7\, \textrm {GPa}$| | |$5.7\, \textrm {GPa}$| |
| |$\eta _\mathrm{f}$| | |$0.00015\, \textrm {Pa} \cdot \textrm {s}$| | |$0.003\, \textrm {Pa} \cdot \textrm {s}$| |
| k | |$1 \times 10^{-12}\, \textrm {m}^2$| | |$1 \times 10^{-12}\, \textrm {m}^2$| |
| |$\phi$| | 0.3 | 0.3 |
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | |$37\, \textrm {GPa}$| | |$37\, \textrm {GPa}$| |
| |$K_\mathrm{d}$| | |$4.8\, \textrm {GPa}$| | |$4.8\, \textrm {GPa}$| |
| |$K_\mathrm{f}$| | |$0.012 \, \textrm {GPa}$| | |$2.25\, \textrm {GPa}$| |
| |$G_\mathrm{d}$| | |$5.7\, \textrm {GPa}$| | |$5.7\, \textrm {GPa}$| |
| |$\eta _\mathrm{f}$| | |$0.00015\, \textrm {Pa} \cdot \textrm {s}$| | |$0.003\, \textrm {Pa} \cdot \textrm {s}$| |
| k | |$1 \times 10^{-12}\, \textrm {m}^2$| | |$1 \times 10^{-12}\, \textrm {m}^2$| |
| |$\phi$| | 0.3 | 0.3 |
4.3 Validation of three-dimensional code implementation
4.3.1 Three-dimensional layered media
The 3-D GPU-accelerated version of FastBiot_QS was validated through a comparison against an analytical solution for 3-D layered poroelastic media, as referenced in White et al. (1975), Carcione & Picotti (2006), Quintal et al. (2009) and depicted in Fig. 2(a). The material parameters used in the simulation are given in Table 3. Following a similar approach to earlier validations, the Kramers–Kronig relations were employed to assess the implementation. This validation process confirms the accuracy of the 3-D numerical code, as demonstrated in Fig. 3.
(a) Sketch illustrating the geometry of the 3-D poroelastic layers. (b) Sketch illustrating the geometry of the (mesoscopic) pore space (intersecting fractures). (c) Sketch illustrating the geometry of the (mesoscopic) pore space (intersecting fractures) which is inclined by 45 deg.
Numerical results calculated by using the present 3-D GPU solver FastBiot_QS compared against the exact analytical solutions for 3-D layered media: (a) the real part of the |$H(f_i)$| component (the P-wave modulus) and (b) corresponding dimensionless attenuation |$1/Q(f_i)$|.
Properties of the model used in the numerical simulations: layered media.
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 1 |
| |$K_\mathrm{f}$| | |$10^{-5}$| | 0.2 |
| |$G_\mathrm{d}$| | 1 | 1 |
| |$\eta _\mathrm{f}$| | 1 | 1000 |
| k | 1 | 1 |
| |$\phi$| | 0.2 | 0.2 |
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 1 |
| |$K_\mathrm{f}$| | |$10^{-5}$| | 0.2 |
| |$G_\mathrm{d}$| | 1 | 1 |
| |$\eta _\mathrm{f}$| | 1 | 1000 |
| k | 1 | 1 |
| |$\phi$| | 0.2 | 0.2 |
Properties of the model used in the numerical simulations: layered media.
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 1 |
| |$K_\mathrm{f}$| | |$10^{-5}$| | 0.2 |
| |$G_\mathrm{d}$| | 1 | 1 |
| |$\eta _\mathrm{f}$| | 1 | 1000 |
| k | 1 | 1 |
| |$\phi$| | 0.2 | 0.2 |
| Parameter | Layer 1 | Layer 2 |
|---|---|---|
| |$K_\mathrm{ g}$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 1 |
| |$K_\mathrm{f}$| | |$10^{-5}$| | 0.2 |
| |$G_\mathrm{d}$| | 1 | 1 |
| |$\eta _\mathrm{f}$| | 1 | 1000 |
| k | 1 | 1 |
| |$\phi$| | 0.2 | 0.2 |
4.4 Validation of FastBiot_QS against a three-dimensional finite element solver
The 3-D GPU-accelerated code, FastBiot_QS, was rigorously validated against an FE solver for poroelastic equations in the displacement–pressure (|$u-p$|) formulation (Quintal et al. 2011). In this FE solver, poroelastic equations were integrated into the COMSOL Multiphysics partial differential equation module using a weak form representation. The entire spatial domain is discretized using an unstructured mesh composed of tetrahedral elements. The system of linear equations is solved using a direct PARDISO solver (Schenk & Gärtner 2004). A direct relaxation test is conducted to determine the |$C_{33}$| component (the P-wave modulus) of the stiffness matrix in Voigt notation. The methodology involves applying a displacement boundary condition of |$u = 10^{-6}\, \exp (i\omega t)$| to the top external wall in a perpendicular direction, while setting normal displacements to zero on other model walls. Initial displacement conditions are zeroed out. Resultant stress and strain calculations on the FE mesh are spatially averaged for each frequency. The complex-valued |$C_{ij}(\omega )$| component is then computed for each frequency. A thorough description of the 3-D boundary conditions can be found in several key references, such as Alkhimenkov et al. (2020a, b).
Two distinct geometries were utilized for this comparative analysis. The first geometry comprises three orthogonal fractures embedded within a poroelastic background material, as illustrated in Fig. 2(b). The second geometry mirrors the first but is inclined along the x-axis by 45 deg, a setup designed to validate the 3-D code’s performance in scenarios where the subdomains do not align with the mesh (Fig. 2c). The material parameters used in the simulation are given in Table 4. The mesh of the (mesoscopic) pore space (voxel-type elements) used in the GPU solver FastBiot_QS is presented in Fig. 4. The resolution of the model is |$N=319^3\approx 32\,000\,000$| grid cells. Figs 5 and 6 present the results of the 3-D GPU solver FastBiot_QS against those of the FEM solver as a function of the dimensionless frequency f. Despite the complex 3-D geometry, the overall agreement is excellent.
Sketch illustrating the mesh of the (mesoscopic) pore space used in the GPU solver FastBiot_QS. Panel (a) shows the full model. Panels (b)–(f) show zoomed in portions of the model. The voxel elements are clearly visible in panel (f).
Numerical results calculated by using the present GPU solver FastBiot_QS and the finite element method (FEM) for the geometry of the (mesoscopic) pore space presented in Fig. 2(b). Panel (a) shows the real part of the |$H(f_i)$| component (the P-wave modulus). Panel (b) shows the error magnitudes of H, expressed as percentages.. Panel (c) shows the corresponding dimensionless attenuation |$1/Q(f_i)$|.
Numerical results calculated by using the present GPU solver FastBiot_QS and the FEM for the geometry of the (mesoscopic) pore space presented in Fig. 2(c). Panel (a) shows the real part of the |$H(f_i)$| component (the P-wave modulus). Panel (b) shows the error magnitudes of H, expressed as percentages.. Panel (c) shows the corresponding dimensionless attenuation |$1/Q(f_i)$|.
Properties of the model used in the numerical simulations: stochastic fracture networks
| Parameter | Frame | Fractures |
|---|---|---|
| |$K_g$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 0.001 |
| |$K_\mathrm{f}$| | 0.1 | 0.1 |
| |$G_\mathrm{d}$| | 1.2083 | 0.1913 |
| |$\eta _\mathrm{f}$| | 0.001 | 0.001 |
| k | |$10^{-6}$| | 0.1 |
| |$\phi$| | 0.2 | 0.5 |
| Parameter | Frame | Fractures |
|---|---|---|
| |$K_g$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 0.001 |
| |$K_\mathrm{f}$| | 0.1 | 0.1 |
| |$G_\mathrm{d}$| | 1.2083 | 0.1913 |
| |$\eta _\mathrm{f}$| | 0.001 | 0.001 |
| k | |$10^{-6}$| | 0.1 |
| |$\phi$| | 0.2 | 0.5 |
Properties of the model used in the numerical simulations: stochastic fracture networks
| Parameter | Frame | Fractures |
|---|---|---|
| |$K_g$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 0.001 |
| |$K_\mathrm{f}$| | 0.1 | 0.1 |
| |$G_\mathrm{d}$| | 1.2083 | 0.1913 |
| |$\eta _\mathrm{f}$| | 0.001 | 0.001 |
| k | |$10^{-6}$| | 0.1 |
| |$\phi$| | 0.2 | 0.5 |
| Parameter | Frame | Fractures |
|---|---|---|
| |$K_g$| | 2 | 2 |
| |$K_\mathrm{d}$| | 1 | 0.001 |
| |$K_\mathrm{f}$| | 0.1 | 0.1 |
| |$G_\mathrm{d}$| | 1.2083 | 0.1913 |
| |$\eta _\mathrm{f}$| | 0.001 | 0.001 |
| k | |$10^{-6}$| | 0.1 |
| |$\phi$| | 0.2 | 0.5 |
It is important to note that both solvers may not converge to the exact solution due to the discrete nature of their methodologies. For instance, in the FEM solver, the use of tetrahedral elements can introduce geometric inaccuracies. Similarly, errors may also occur in the GPU solver FastBiot_QS due to its regular voxel grid structure, despite its high density. In conclusion, both solvers deliver satisfactory results, corroborating the efficiency of the newly developed 3-D GPU-based solver FastBiot_QS.
Although the geometries depicted in Figs 2(b) and (c) might appear simplistic, they exhibit complex physical behaviours. These geometries facilitate flow from fractures to the background matrix, a phenomenon clearly evident in Fig. 5 at low frequencies. This interaction is characterized by a peak in attenuation and an increase in stiffness at these lower frequencies. At higher frequencies, flow transitions from fracture-to-fracture, manifesting as a second peak in the attenuation curve and a further increase in stiffness as frequency increases.
4.5 Validation of FastBiot_QS against Kramers–Kronig relations
Fig. 7 presents the results of the 3-D GPU solver FastBiot_QS and the Kramers–Kronig relations calculated for the dimensionless attenuation |$1/Q(f_i)$| (green curve). In the 3-D GPU solver FastBiot_QS, dispersion |$H(f_i)$| and |$1/Q(f_i)$| are calculated independently and, in theory, must satisfy the causality principle, that is the Kramers–Kronig relations (Christensen 1982; Mikhaltsevitch et al. 2016). The Kramers–Kronig relations connect the real and imaginary parts of a complex function, ensuring that they are consistent with the principle of causality. In the context of quasi-static poroelasticity, these relations link the moduli dispersion (|$H(f_i)$|) and attenuation (|$1/Q(f_i)$|). As shown in Fig. 7, the dimensionless attenuation |$1/Q(f_i)$| calculated via the 3-D GPU solver FastBiot_QS and the Kramers–Kronig relations are in excellent agreement. It serves as a critical verification step, as moduli and attenuation are computed independently in the time domain in our approach, unlike frequency-domain solvers where these parameters are inherently linked. This step is particularly valuable for ensuring the consistency and reliability of our methodology.
Numerical results calculated by using the present GPU solver FastBiot_QS for the two geometries of the (mesoscopic) pore space presented in Figs 2(b) and (c). Panels (a) and (c) show the real part of the |$H(f_i)$| component (the P-wave modulus). Panels (b) and (d) show the corresponding dimensionless attenuation |$1/Q(f_i)$|. The dashed green curve (marker - cross) corresponds to the approximate Kramers–Kronig relations.
5 RESULTS
In this section, we present the essential information about fracture (crack) density and the results of five different simulations conducted to evaluate the seismic attenuation and modulus dispersion in fractured poroelastic media. Each simulation corresponds to a different configuration of fractures within the media. They include: (1) 54 interconnected fractures with regular spacing, (2) 54 non-connected fractures with regular spacing, (3) 54 fractures with a log-normal distribution of their positions and random orientations, (4) 54 fractures grouped into four clusters with respect to their positions and random orientations and (5) 54 fractures with a uniform distribution of their positions and random orientations.
5.1 Analysis of fracture configurations
5.1.1 Fracture (crack) density definitions
The elastic properties of cracked materials can be described as functions of crack density parameters. In this study, we explore fractured media but the theory developed for cracks can be also applied to fractures. In the case of isotropic random fracture orientations, the scalar crack density, introduced by Bristow (1960), is defined as:
where |$r_k$| are the radii of the cracks, k is the number of cracks and |$V$| represents the volume. For non-random orientations of fractures (cracks), this has been extended to the second-rank crack density tensor |$a$|, defined by:
where |$\mathbf {\mathit{ n}}$| is a unit normal to the |$k$|-th fracture (crack). Additionally, in the context of elastic properties, the fourth-rank crack density tensor |$b$| is used:
This tensor has a relatively minor role for traction-free fracture. The contribution of a fracture to these parameters is proportional to the cube of its size. In all simulations, we consider a volume |$V$| of 1 (normalized volume) containing 54 fractures, with a fracture radius |$r$| given by |$\frac{56}{511} \approx 0.1096$|.
5.1.2 Fracture density analysis
Fracture density is a key parameter in the 3-D simulations, as it fundamentally controls the stiffness and fluid pressure diffusion (FPD) behaviour within the system. Specifically, the density of fractures influences the medium’s connectivity and flow characteristics:
For low-fracture densities with uniformly distributed fractures, the likelihood of fracture intersections is minimal, resulting in negligible fracture-to-fracture flow.
For high-fracture densities, the presence of numerous intersecting fractures facilitates significant fracture-to-fracture flow, enhancing the overall connectivity of the system.
In this study, fracture density is kept constant across all 3-D simulations to provide a consistent basis for analysing the effects of other parameters.
5.1.2.1 Random orientations of fractures
In the case of random orientations of fractures, we can use the scalar fracture density definition |$\rho$|. Since all fractures have the same radius, the summation can be simplified:
Thus, the scalar crack density, |$\rho$|, is approximately 0.071 06.
5.1.2.2 Fractures oriented parallel to XY- and XZ-planes
There are 27 fractures in the |$\mathit{ XY}$| plane and 27 in the |$\mathit{ XZ}$| plane. The normal vectors for the cracks are as follows: for fractures in the |$XY$| plane, |$\mathbf {\mathit{ n}} = (0, 0, 1)$|; for fractures in the |$\mathit{ XZ}$| plane, |$\mathbf {\mathit{ n}} = (0, 1, 0)$|. For 27 fractures in the |$XY$| plane and 27 in the |$\mathit{ XZ}$| plane:
. Thus, the second-rank fracture (crack) density tensor |$a$| is:
The fourth-rank fracture (crack) density tensor |$b$| is:
and
Thus, the tensor |$b$| is:
5.1.3 Interconnected fractures
The first simulation involves a 3-D fracture network where all fractures are interconnected (Figs 8a–c). The resulting attenuation curve reveals two distinct peaks (Fig. 11). The first peak at lower frequencies is attributed to FPD between fractures and the surrounding matrix, indicative of fracture-to-background flow. This mechanism occurs because fluid pressure equilibration between fractures and the surrounding matrix is a slower process, dominating the lower frequency range. The second peak at higher frequencies is associated with fracture-to-fracture flow, reflecting fluid pressure equilibration within the interconnected fracture network. Interconnected fractures facilitate faster pressure equilibration, resulting in attenuation at higher frequencies.
3-D crack network of interconnected cracks. (a) XZ-slice of the 3-D model. (b) YZ-slice of the 3-D model. (c) Full 3-D model consisting of 54 interconnected cracks.
Sketch illustrating the numerical resolution used in the simulations. The model with non-connected fractures is shown.
3-D crack networks used in the simulations. (a) Log-Normal distribution of fracture. (b) Four clusters of fractures. (c) Uniform distribution of fractures.
Numerical results calculated by using the present 3-D GPU solver FastBiot_QS for five different fracture network distributions: (a) the real part of the |$H(f_i)$| component (the P-wave modulus) and (b) corresponding dimensionless attenuation |$1/Q(f_i)$|.
5.1.4 Non-connected fractures
The second simulation considers a 3-D fracture network with non-connected fractures (Fig. 8 d–f). The mesh of the (mesoscopic) pore space (voxel-type elements) used in the GPU solver FastBiot_QS is presented in Fig. 9. The resolution of the model is |$N=511^3\approx 133\,000\,000$| grid cells. The attenuation curve shows a single prominent peak at lower frequencies (Fig. 11). The absence of a second peak indicates that the primary attenuation mechanism is fracture-to-background flow, as there are no interconnected pathways for fracture-to-fracture fluid flow. The attenuation in this scenario is significantly higher compared to interconnected cracks because the isolated fractures cannot facilitate efficient pressure diffusion within the network. This highlights the impact of fracture connectivity on the attenuation behaviour and suggests that non-connected fractures lead to increased energy dissipation at lower frequencies.
5.1.5 Log-normal distribution of fractures
The third simulation involves a 3-D fracture network with a log-normal distribution of fracture positions within the domain V and uniform orientations of fractures (Fig. 10a). The results indicate that this distribution results in complex attenuation behaviour, exhibiting characteristics of both interconnected and non-connected crack networks (Fig. 11b). Due to the uniform orientations of fractures the overall dispersion and attenuation are smaller compared to the previous two scenarios where 27 fractures are parallel to XY-plane. However, at higher frequencies (|$f> 10^2$|), the attenuation is higher compared to the simulation of non-connected cracks. This indicates that some degree of connectivity is present.
5.1.6 Four clusters of cracks
The fourth simulation examines a 3-D crack network consisting of four distinct clusters of cracks (Fig. 10b). Each cluster represents a localized region of higher fracture density. Within each cluster, fracture-to-fracture flow dominates at higher frequencies, leading to stronger attenuation peaks compared to log-normal distribution (Fig. 11b). At low frequencies, fracture-to-background flow is the primary mechanism and the attenuation is quite small reflecting the clustered nature of the fracture network.
5.1.7 Uniform distribution of cracks
The fourth simulation examines a 3-D crack network consisting of fractures with random positions (Fig. 10c). The overall behaviour of the attenuation curve is similar to log-normal distribution of fractures but differ form the simulation where four clusters of fractures were considered.
5.2 Fluid pressure field
5.2.1 Interconnected versus non-connected fracture networks
Fig. 12 shows snapshots of the fluid pressure field for interconnected (a), (b) and non-connected (c), (d) fracture networks at two different frequencies, |$f_i=10$| and |$f_i=1000$|. At the low frequency |$f_i=10$|, fracture-to-background flow is evident, as shown in Figs 12(a) and (c). Note that in the interconnected fracture network, both horizontal and vertical fractures exhibit high-fluid pressure, while in the non-connected fracture network, only horizontal fractures exhibit high-fluid pressure. In both fracture networks, strong fluid pressure gradients near the walls of the fractures indicate fracture-to-background flow. At the higher frequency |$f_i=1000$|, in the interconnected fracture network, strong-fluid pressure gradients appear in the vertical fractures, indicating that this frequency is close to the second attenuation peak (fracture-to-fracture flow). Conversely, in the non-connected fracture network, no significant fluid pressure diffusion occurs.
Numerical results calculated using the present 3-D GPU solver FastBiot_QS for interconnected (a), (b) and non-connected (c), (d) fracture network distributions at two different frequencies |$f_i=10$| and |$f_i=1000$|: the fluid pressure |$p_\mathrm{f}$| field is shown as a 3-D view projected on the fracture surface (left) and as an XZ-plane view (right).
5.2.2 Four clusters of fractures
Fig. 13 shows snapshots of the fluid pressure field for the model of four clusters of fractures at two different frequencies, |$f_i=10$| and |$f_i=1000$|. At the low frequency |$f_i=10$|, the fluid pressure is high in the fractures, and fracture-to-background flow is taking place (Fig. 13c). At the higher frequency |$f_i=1000$|, fracture-to-fracture flow occurs. The fluid pressure is high in the nearly horizontal fractures, and fluid pressure gradients are present in the vertical fractures (Figs 13d and e).
Numerical results calculated using the present 3-D GPU solver FastBiot_QS for the model of four clusters of fractures at two different frequencies |$f_i=10$| and |$f_i=1000$|: the fluid pressure |$p_\mathrm{f}$| field is shown as a 3-D view projected on the fracture surface (a) and (d) and as an XZ-plane view (b, c) and (e, f).
6 DISCUSSION
6.1 Earlier two-dimensional studies of double-porous media
One of the first modelling studies to explore the geophysical relevance of double-porous media is Masson & Pride (2007), who presented a novel approach for modelling quasi-static poroelastic equations and investigated different geometries of double-porous materials. Subsequently, a frequency-domain finite element approach was developed by Rubino et al. (2009). Since then, many aspects of double-porous media within the framework of poroelasticity have been extensively studied.
An early study specifically examining fractures in poroelastic media is Rubino et al. (2013), which identified two distinct attenuation peaks: the first resulting from FPD between fractures and the background, and the second from FPD within interconnected fractures. These two peaks, along with their physical interpretations, were further investigated in 2-D by Rubino et al. (2014).
A stochastic fracture network in 2-D was rigorously analysed by Hunziker et al. (2018), where the two attenuation peaks were examined in detail. Additionally, Hunziker et al. (2018) conducted a 3-D study in COMSOL, comparing the results of a single fracture set in 3-D with those of a single fracture set in 2-D. This type of study was facilitated by the adaptive mesh refinement in 2-D, presented comprehensively by Favino et al. (2020).
These studies provide a solid foundation for extending the analysis to 3-D stochastic fracture networks, as demonstrated in this work.
6.2 Simulation cases and insights
This study investigates several 3-D fracture network configurations to demonstrate the capabilities of the proposed solver and validate its application to diverse geometries. While these simulations provide valuable initial insights, the scope of this work does not include exhaustive sensitivity analyses or comprehensive explorations of all configurations, which are planned for future research. These cases collectively underscore the solver’s robustness and versatility, establishing a foundation for future studies aimed at gaining broader and more generalizable insights through detailed sensitivity analyses and parameter variations.
6.2.1 Log-normal distribution of fracture positions
The log-normal distribution was selected as a representative statistical model commonly used in fracture network studies. Although real geological settings exhibit diverse fracture distributions, this case serves as a proof-of-concept for the solver’s ability to handle statistically distributed fractures.
6.2.2 Fracture clusters
The simulation of four clusters of fractures illustrates the solver’s ability to model complex, clustered configurations. This case highlights FPD within connected fractures inside clusters, resulting in two attenuation peaks. While broader conclusions about fracture clusters cannot be drawn from a single case study, this simulation provides a foundation for future work involving multiple configurations and parameter variations.
6.2.3 Random fracture positions
The simulation of randomly distributed fractures demonstrates the solver’s flexibility in handling stochastic fracture arrangements. It highlights the solver’s applicability to modelling irregular and non-deterministic fracture networks, which are common in geological formations.
6.3 Fracture clustering and stochastic spatial distribution of fractures
The results of these simulations underscore the significant impact of fracture clustering and distribution on seismic attenuation and modulus dispersion in 3-D fractured media. This observation suggests that the relationship between fracture clusters and fluid pressure diffusion warrants further investigation. The presence of clusters of fractures facilitates rapid fluid pressure equilibration within the fracture network, leading to stronger attenuation at higher frequencies. Conversely, non-connected fractures primarily contribute to fracture-to-background flow, resulting in higher attenuation at lower frequencies.
To further illustrate the impact of fracture clustering on the attenuation |$1/Q$|, let us plot the normalized |$1/Q$| by the first attenuation peak on a linear scale on the y-axis. Such a plot can provide insights into fracture connectivity. Two end-members, interconnected and non-connected fractures, provide the low and high values of the attenuation curve |$1/Q$| (Fig. 14, black and red curves). The other three simulations fall between these limits. The simulation with four clusters of cracks corresponds to the highest connectivity within the clusters, which is indicated by a prominent second attenuation peak (Fig. 14, pink squares). The other two distributions, log-normal and uniform, show a similar degree of connectivity, with the uniform distribution having slightly more interconnected fractures (Fig. 14, green curve).
Numerical results calculated using the present 3-D GPU solver FastBiot_QS for five different fracture network distributions: dimensionless attenuation |$1/Q(f_i)$| normalized by each maximum at low frequencies (note the linear scale on the y-axis).
The simulation involving clusters of cracks is particularly insightful for understanding how spatial organization within fracture networks influences seismic response. Clusters act as localized regions of high-fracture connectivity and high-fracture density, promoting fracture-to-fracture flow within the clusters, while the overall network still experiences significant fracture-to-background flow between clusters. This dual behaviour has important implications for interpreting seismic data and inferring hydraulic properties in heterogeneous fractured media. In other words, clustering significantly affects the attenuation and dispersion curves, promoting the non-invasive identification of fracture clusters. Such configurations require further investigation.
6.4 A minimum number of voxel elements per fracture aperture
In this study, most of the simulations (i.e. simulations 1–5 in Section 5) were performed at a resolution of |$N=511^3 \approx 133\,000\,000$| voxel elements. Given that each simulation included 54 fractures, it was chosen to have six voxel elements per fracture aperture to ensure proper modelling of fracture-to-fracture flow. To check the accuracy of the solver, we performed an additional simulation for the model of four clusters of fractures but with a lower resolution of |$N=255^3 \approx 17\,000\,000$| voxel elements. This resolution is eight times smaller than the standard resolution employed in this study. Fig. 15 shows the numerical results calculated using the present 3-D GPU solver FastBiot_QS for two different spatial resolutions: |$N=511^3$| (black curve) and |$N=255^3$| (green curve). It can be seen that slight discrepancies are present, but the overall behaviour is similar; that is two attenuation peaks are visible in both simulations. This comparison validates the accuracy of the present high-resolution simulations and also suggests that fewer voxel elements per fracture aperture can be chosen for future studies.
Numerical results calculated by using the present 3-D GPU solver FastBiot_QS for a fracture network distribution consisting of a four clusters of fractures for two different spatial resolutions |$N=511^3$| and |$N=255$|: (a) the real part of the |$H(f_i)$| component (the P-wave modulus) and (b) corresponding dimensionless attenuation |$1/Q(f_i)$|.
6.5 Performance
The present 3-D GPU solver, FastBiot_QS, demonstrates exceptional performance. Below, we explore the performance metrics of GPU acceleration on different hardware.
6.5.1 Nvidia A100 GPU
For 3-D simulations performed on an Nvidia A100 GPU with a resolution of |$N=511^3 \approx 133\,000\,000$| voxel elements, the quasi-static solution of Biot’s poroelastic equations for a single frequency |$f_i$| is achieved in 9 min for the stiffness modulus and 18 min for calculating attenuation |$1/Q$|. These simulations correspond to a travelling wave over a distance of |$20 \times n_x$|, where |$n_x$| is the number of grid cells in the x-dimension. The memory required for this setup is approximately 50.5 GB of DRAM (Dynamic Random-Access Memory).
6.5.2 Nvidia laptop GeForce RTX 4090 GPU
When using an Nvidia GeForce RTX 4090 laptop GPU with 16 GB of DRAM, the performance is slower by a factor of approximately 3. Due to the DRAM limitation, the maximum resolution that can be accommodated is |$N=319^3 \approx 32\,000\,000$| voxel elements, which requires 13 GB of DRAM. For this resolution, the computational time for a single frequency |$f_i$| is 6.7 min for the stiffness modulus and 13.4 min for calculating attenuation |$1/Q$|. Despite the reduced resolution, the laptop GPU is sufficient to fully reproduce the results of this study.
6.5.3 GPU acceleration compared to CPU implementation
The A100 GPU achieves an acceleration factor of approximately |$520\times$| compared to a CPU implementation (a fully vectorized MATLAB version comparable in computational time to a C implementation).
The Laptop GeForce RTX 4090 GPU achieves an acceleration factor of approximately |$170\times$| over the CPU.
Direct comparison with FEM solvers, such as COMSOL with a Pardiso direct solver, is challenging due to the fundamentally different approaches between FEM and the voxel-based GPU solver. However, for the same number of elements in COMSOL and voxels in the FastBiot_QS solver running on a single A100 GPU, the acceleration exceeds approximately |$500\times$|.
The key advantage of the FastBiot_QS solver lies in its scalability. It can handle resolutions exceeding 32 million voxel elements even on a laptop GPU, a capability that COMSOL with a direct solver cannot achieve.
6.6 Future research directions
The future research directions proposed in this study are centered on the application and further modification of the present 3-D solver. While some foundational concepts related to fracture networks and poroelastic processes have been explored in 2-D studies, they do not fully capture the complexities of 3-D fracture networks and the associated physical phenomena. Our future work aims to extend these concepts into 3-D, leveraging the unique capabilities of the proposed solver to address unresolved challenges and provide novel insights into fracture dynamics. This approach represents a significant step forward in understanding the behaviour of 3-D fracture networks, beyond what has been achieved in existing 2-D studies. Several avenues for future research can further enhance the capabilities and applications of the developed methodologies.
6.6.1 Incorporation of complex fracture geometries
Future studies should focus on extending the current model to incorporate more complex and realistic fracture geometries. This includes irregularly shaped fractures, varying aperture sizes, and multiscale fracture networks. Such advancements will improve the model’s applicability to a wider range of geological formations and provide more accurate predictions of their behaviour under various stress and fluid flow conditions.
6.6.2 Wave propagation
The current model addresses quasi-static equations. Incorporating dynamic conditions, such as seismic waves in poroelastic media, will be crucial for simulating real-world scenarios more accurately. Such algorithms already exist (Alkhimenkov et al. 2021) but have never been applied to study stochastic fracture networks in 3-D poroelastic formations.
6.6.3 Anisotropy
Investigating the impact of anisotropy in the mechanical and hydraulic properties of the porous media can provide deeper insights into the behaviour of fractured rocks. Future research should aim to integrate anisotropic material properties into the model, allowing for a more comprehensive understanding of directional dependencies in seismic responses and fluid flow patterns.
6.6.4 Real-field application
To enhance the practical applicability of the model, future research should focus on integrating real-field data into the simulations. This includes using data from field measurements and laboratory experiments to validate and calibrate the model.
6.6.5 Machine learning and data-driven approaches
Incorporating machine learning and data-driven approaches can significantly enhance the efficiency and predictive capabilities of the model. Future studies can explore the use of machine learning algorithms to identify patterns in simulation data, optimize computational resources, and predict the behaviour of fractured media under various conditions.
6.6.6 Multiphysics coupling
Extending the model to include multiphysics coupling, such as thermo-hydromechanical processes, will provide a more holistic understanding of the interactions between different physical phenomena in fractured rocks. This will be particularly useful for applications in geothermal energy extraction, where temperature changes play a critical role in fluid flow and rock mechanics.
7 CONCLUSIONS
This study addresses the complex challenge of characterizing 3-D fractured poroelastic media at the mesoscopic scale by leveraging advanced numerical methods and GPU-accelerated computing. A novel 3-D numerical approach has been developed to model the effective viscoelastic properties of hydro-mechanically coupled systems using quasi-static Biot’s poroelastic equations on GPUs. This method enables high-resolution simulations of up to 133 million voxel elements (|$512^3$|) on a single GPU, providing unprecedented detail in analysing fluid pressure diffusion and its effects on seismic attenuation and modulus dispersion within stochastic fracture networks.
The GPU-accelerated solver, FastBiot_QS, delivers exceptional performance, achieving a computational speed-up of approximately 520 times compared to CPU-based methods. Validation against 1-D and 3-D analytical solutions, as well as, against a 3-D finite element solver, confirms the solver’s accuracy in capturing complex poroelastic behaviours. Even at lower resolutions, the solver maintains high accuracy, suggesting a practical pathway for balancing computational cost and resolution in future studies without compromising fidelity.
Simulations reveal that fracture spatial arrangements significantly influence seismic attenuation and modulus dispersion. Densely clustered fractures enhance attenuation at higher frequencies due to localized fluid pressure equilibration, while sparsely distributed fractures contribute to lower frequency attenuation, reflecting increased fracture-to-background flow. Intermediate behaviours observed for log-normal and uniform fracture distributions highlight the impact of fracture density and connectivity on seismic responses.
The non-dimensional framework used in this study allows the results to be generalized across various rock types and fracture geometries. This adaptability makes the findings relevant for characterizing hydraulic properties in diverse geological settings. These findings provide valuable insights for applications such as geothermal energy, hydrocarbon exploration, nuclear waste disposal and |$\mathrm{CO}_2$| storage. The GPU-accelerated framework establishes a new standard for high-resolution modelling, enabling future research on larger and more complex fracture networks.
ACKNOWLEDGMENTS
Yury Alkhimenkov gratefully acknowledges support from the Swiss National Science Foundation, project number P500PN_206722. I thank Beatriz Quintal for the fruitful discussions that inspired this work and for providing the 1-D analytical solution used to validate the numerical solver.
DATA AVAILABILITY
Routines that can be used to reproduce the presented results are available from a permanent DOI repository (Zenodo) (Alkhimenkov 2024b) at 10.5281/zenodo.14357095.
