SUMMARY
Carbonates are highly heterogeneous and commonly develop large-scale fracture-cavity systems due to the strong diagenesis effect of dissolution, faulting and karstification. Understanding the seismic wave velocity and attenuation signatures is crucial for hydrocarbon exploration and recovery, carbon capture, geothermal exploitation and groundwater management using seismic methods. Nevertheless, due to the sample size limitation, traditional core-scale experiments and sonic logs fail to effectively unravel the essential factors controlling elastic and anelastic responses of those carbonates with large-scale heterogeneity. We developed a 2-D geology-constrained method for constructing digital rock models of carbonates with fracture-cavity systems, integrating outcrop data and statistical analyses. Using finite element simulations based on Biot's quasi-static equations of poroelastic consolidation, we examine the effects of cavity shape, size and number, as well as fracture number, length, width and angle, on wave dispersion and attenuation characteristics. Deep carbonate rocks with fracture-cavity systems exhibit wideband attenuation spanning over three orders of magnitude (from < 0.01 to >10 Hz), driven by coupled pore pressure relaxation mechanisms in the heterogeneous system. Horizontal attenuation was found to be roughly five times greater than vertical attenuation (as fractures are predominantly in a semi-vertical direction), highlighting significant attenuation anisotropy. While cavities have minimal effects on attenuation, fracture density and geometry—especially longer and narrower fractures—significantly influence wave dispersion and attenuation behaviour. These results underscore the importance of integrating both seismic velocity and attenuation data to enhance reservoir characterization reliability, and highlight the potential of using comprehensive broad-band seismic data and large-offset seismic data to improve geophysical interpretations of complex deep carbonate reservoirs.
1 INTRODUCTION
Carbonate rocks are widely distributed across the globe (White 1988; Palmer 1991; Ford & Williams 2007; Willimas 2008). Characterized by significant heterogeneity, carbonate reservoirs in the subsurface often form complex fracture-cavity systems as a result of diagenetic processes such as dissolution, faulting and karstification (Akbar et al. 2000; Lucia et al. 2003; Ahr 2011; Dreybrodt 2012). For instance, carbonate reservoirs with fracture-cavity systems serve as the primary hydrocarbon-bearing spaces at depths exceeding 6 km in the Tarim Basin, Northwestern China (Zhu et al. 2015; Li et al. 2016; Jun & Jibiao 2021; Wang et al. 2021). Beyond that, carbonate reservoirs with fracture-cavity systems also play a critical role in carbon capture and storage (Boot-Handford et al. 2014), contribute to the global carbon cycle (Gattuso et al. 1998; Ridgwell & Zeebe 2005), and are essential to efforts in karst rocky desertification control (Wang et al. 2019). Therefore, understanding the impact of the fracture-cavity systems on the geophysical response of carbonate rocks is crucial for assessing hydrocarbon reservoir potential (Ahr 2011), reducing production risks (Li et al. 2024), enhancing recovery (Zheng et al. 2019) and promoting sustainable development (Zhao et al. 2013).
A range of researches have demonstrated that carbonate reservoirs might exhibit significant and complex dispersion and attenuation behaviours, primarily driven by wave-induced fluid flow mechanisms across multiple scales, especially within heterogeneous rock formations and in the presence of fractures (Adam et al. 2009; Borgomano et al. 2019; Li et al. 2020). The dispersion and attenuation of carbonate rocks have been investigated through laboratory measurements (Batzle et al. 2006; Mikhaltsevitch et al. 2016; Borgomano et al. 2017; Li et al. 2020); while other studies have extracted attenuation from field seismic and well log data (Cheng et al. 1982; Bouchaala et al. 2016). The analysis of dispersion and attenuation characteristics offers valuable insights into various reservoir properties, including fluid distribution (Agersborg et al. 2008; Li et al. 2020), porosity (Pang et al. 2019), connectivity (Conti et al. 2019) and infills (Kong et al. 2013), rock composition (Verwer et al. 2008), presence of fractures (Quintal et al. 2014; Lissa et al. 2019).
The exploration of carbonate reservoirs with fracture-cavity systems can be conducted across a broad frequency range (Fig. 1), utilizing both active sources such as exploration seismology (Wang et al. 1998; Warrlich et al. 2019; Hendry et al. 2021), vertical seismic profiling (VSP) (Akbar et al. 2000; Khoshbakht et al. 2012) and distributed acoustic sensing (DAS), as well as passive sources like earthquake seismology and microseismology. Each of these frequency bands provides valuable insights and contributes to a comprehensive understanding of the carbonate reservoirs.
Schematic diagram of broad-band geophysical exploration methods for deep carbonate rocks. These methods include distributed acoustic sensing (DAS), which covers frequencies from below 10−4 Hz to above 103 Hz; earthquake seismology from around 10−4 to 100 Hz; microseismology from approximately 10−1 to 101 Hz; exploration seismology from below 101 Hz to over 102 Hz; and vertical seismic profiling (VSP), spanning frequencies from below 102 Hz to over 103 Hz. The REV (representative elementary volume) illustrates typical fractured-cavity systems within deep carbonate reservoirs. DAS uses fibre-optic cables to monitor subsurface acoustic changes by capturing backscattered laser signals, while both active and passive seismic sources generate waves that interact with subsurface geological structures.
To effectively characterize fracture-cavity systems using broad-band geophysical data, it is essential to apply rock physics concepts to bridge the gap between geological features and the corresponding elastic and attenuation properties. Nevertheless, the traditional rock physics approach based on logging and laboratory measurement is nearly infeasible, due to the size of the representative elementary volume (REV) of fracture-cavity carbonate (as shown in Fig. 1), which often extends to metres or even larger. Digital rock physics models offer a breakthrough solution by enabling us to model larger scales and systematically explore the elastic and attenuation properties associated with fracture-cavity interactions, allowing for a detailed decoupling of their effects and providing clearer insights into their complex behaviours (Popov et al. 2009; Dvorkin et al. 2011; Kalam 2012; Andrä et al. 2013; Sell et al. 2016; Saad et al. 2018; Reinhardt et al. 2022). Another key advantage of digital rock physics is its ability to combine actual outcrop images with mathematical methods to generate a large amount of reasonable digital models for statistical analysis. This capability complements traditional methods, enabling the analysis of elastic behaviour through broader and more generalized modes.
Several numerical simulation methods have been developed to model the elastic and attenuation properties of rocks. The finite element static simulation method excels in calculating static elastic moduli by minimizing elastic potential energy (Garboczi 1998; Arns et al. 2002; Grechka & Kachanov 2006; Zhang & Toksöz 2012b; Alkhimenkov 2024). In contrast, the transmitted wave simulation measures the time delays of transmitted waves to derive effective velocities (Saenger & Shapiro 2002; Madonna et al. 2012). The dynamic stress–strain simulation takes a more comprehensive approach by solving coupled wave equations with constitutive equations to simulate stress–strain evolution, allowing for the computation of velocity and attenuation from amplitude and phase differences (Das et al. 2019; Zhu et al. 2023; Yang et al. 2024). Other studies have enhanced the understanding of squirt flow effects in porous and fractured rocks through 3-D numerical modelling (Zhang & Toksöz 2012a; Alkhimenkov et al. 2020a, b) and GPU-based simulations (Alkhimenkov 2024). Additionally, the simulations based on quasi-static Biot's equations are used to analyse the seismic attenuation and dispersion in heterogeneous porous media (Masson & Pride 2007; Quintal et al. 2011; Lissa et al. 2019, 2021), offering detailed insights into complex geological systems. Each method addresses specific aspects of elastic and attenuation properties in rocks, and selecting an appropriate simulation method requires careful consideration of the specific geological context, as well as the properties and physical mechanisms of interest. We chose to use the method based on Biot's quasi-static equations of consolidation (Quintal et al. 2011), which considers only the fluid pressure relaxation effect and neglects the inertia term. By translating microscale features into equivalent rock physical parameters, Biot's equations are suitable for large-scale models. Using Biot's equations also allows us to account for pressure relaxation within both the complex fracture-cavity system and the background pores. Additionally, this approach enhances computational efficiency and stability, enabling the generation of sufficient data for statistical analysis.
In this study, we develop a geophysical characterization method for carbonate reservoirs with large-scale fracture-cavity systems using digital rock physics. Based on outcrop data combined with mathematical methods, we constructed numerous digital rock models and applied finite element simulations to calculate the equivalent elastic and attenuation properties. Through comprehensive statistical analysis of the numerical results, we evaluated the effects of cavity shape, size and number, as well as fracture density, length, width and angle, on the wave dispersion and attenuation characteristics.
2 METHODOLOGY
2.1 Geology-constrained digital rock model construction
As shown in the left corner of Fig. 2, the deep carbonate rocks with the fracture-cavity system are commonly delineated into two zones: the slip fault zone (intensely deformed damage zone, here referred as Zone 1), and the induced fracture zone (weakly deformed damage zone, here referred as Zone 2) (Micarelli et al. 2006). The presence of dense fractures in the slip fault zone (Zone 1) significantly increases the probability of dissolution, leading to the formation of a varying number of dissolution cavities.
Workflow for constructing the digital carbonate model with a fracture-cavity system. Cavities and fractures are digitized from outcrop data and geological information. Cavities are defined by their number, location and shape, while fractures in Zone 1 and Zone 2 are characterized by a fixed number, length and width, with small-range randomness applied to their location and angle to better represent natural variability. Outcrop images are after Wang 2015 and Wang et al. 2021.
In order to statistically map the relationship between the physical parameters of dissolution cavities and fractures to the velocity and attenuation characteristics of the rock, we developed a digital model construction method based on geological constraints. Based on the actual dissolution geological outcrops and additional geological information, we generated the digital rock model with the fracture-cavity system to portray the intricate characteristics of the deep carbonate reservoirs, with the workflow illustrated in Fig. 2. The morphology of the dissolution cavities is based on real geological outcrops, such as carbonate outcrops in the Tarim Basin, China (Wang et al. 2021). Parameters such as the number of dissolved cavities, as well as the number, length and width of fractures within the slip fault zone (Zone 1) and induced fracture zone (Zone 2) are determined. The digital rock model measures 12 × 10 m, and the basic parameters for constructing the models are shown in Table 1 without further explanation.
| Cave | Fractures | |||
|---|---|---|---|---|
| Slip fault zone (Zone 1) | Induced fracture zone (Zone 2) | |||
| Area | |$\approx 4.1620{\rm{\ }}{{{\rm{m}}}^2}$| | Number | 50 | 15 |
| Length | 4.5 m | Length | 0.6 (m) |$\ \times $| 2 | 1 (m) |$\ \times $| 2 |
| Width | 2 m | Width | 0.02 (m) |$\ \times $| 2 | 0.02 (m) |$\ \times $| 2 |
| Angle | |${{\phi }_1} \approx 70^\circ $| |$\Delta {{\phi }_1} \approx 5^\circ $| | |${{\phi }_2} \approx 90^\circ $| |$\Delta {{\phi }_2} \approx 2^\circ $| | ||
| Cave | Fractures | |||
|---|---|---|---|---|
| Slip fault zone (Zone 1) | Induced fracture zone (Zone 2) | |||
| Area | |$\approx 4.1620{\rm{\ }}{{{\rm{m}}}^2}$| | Number | 50 | 15 |
| Length | 4.5 m | Length | 0.6 (m) |$\ \times $| 2 | 1 (m) |$\ \times $| 2 |
| Width | 2 m | Width | 0.02 (m) |$\ \times $| 2 | 0.02 (m) |$\ \times $| 2 |
| Angle | |${{\phi }_1} \approx 70^\circ $| |$\Delta {{\phi }_1} \approx 5^\circ $| | |${{\phi }_2} \approx 90^\circ $| |$\Delta {{\phi }_2} \approx 2^\circ $| | ||
| Cave | Fractures | |||
|---|---|---|---|---|
| Slip fault zone (Zone 1) | Induced fracture zone (Zone 2) | |||
| Area | |$\approx 4.1620{\rm{\ }}{{{\rm{m}}}^2}$| | Number | 50 | 15 |
| Length | 4.5 m | Length | 0.6 (m) |$\ \times $| 2 | 1 (m) |$\ \times $| 2 |
| Width | 2 m | Width | 0.02 (m) |$\ \times $| 2 | 0.02 (m) |$\ \times $| 2 |
| Angle | |${{\phi }_1} \approx 70^\circ $| |$\Delta {{\phi }_1} \approx 5^\circ $| | |${{\phi }_2} \approx 90^\circ $| |$\Delta {{\phi }_2} \approx 2^\circ $| | ||
| Cave | Fractures | |||
|---|---|---|---|---|
| Slip fault zone (Zone 1) | Induced fracture zone (Zone 2) | |||
| Area | |$\approx 4.1620{\rm{\ }}{{{\rm{m}}}^2}$| | Number | 50 | 15 |
| Length | 4.5 m | Length | 0.6 (m) |$\ \times $| 2 | 1 (m) |$\ \times $| 2 |
| Width | 2 m | Width | 0.02 (m) |$\ \times $| 2 | 0.02 (m) |$\ \times $| 2 |
| Angle | |${{\phi }_1} \approx 70^\circ $| |$\Delta {{\phi }_1} \approx 5^\circ $| | |${{\phi }_2} \approx 90^\circ $| |$\Delta {{\phi }_2} \approx 2^\circ $| | ||
The fracture locations |$( {x,y} )$| are generated based on defined distribution rules. The normal distribution N with the mean of the model centre is used to obtain the x-coordinates of the fractures in the two zones, while a uniform distribution U is employed to determine their y-coordinates:
Here, |${{L}_x}$| and |${{L}_y}$| are the horizontal and vertical size of the model, respectively; |$\sigma $| denotes the horizontal variance and can be adjusted as needed.
The angle of the fractures can be introduced following a specific distribution with some degree of randomness. In the slip fault zone (Zone 1), there are usually complex conjugate fracture groups (Storti et al. 2003; Agosta & Aydin 2006; Agosta et al. 2012; Wu et al. 2019; Wu et al. 2020), and we assume that the mean angle of one set of conjugate fractures is |${{\phi }_1}$|, while the mean angle of another set of conjugate fractures is |$( {\pi - {{\phi }_1}} )$|. By constructing a mixed distribution of two uniformly distributions, we obtain the fracture angle |${{{\rm{\Phi }}}_{\textit{zone}1}}$| distribution centred on either of the two angles with a deviation of |$\Delta {{\phi }_1}$|:
For fractures within the induced fracture zone (Zone 2), the fracture angle |${{{\rm{\Phi }}}_{\textit{zone}2}}$| is uniformly distributed within a range of |$\Delta {{\phi }_2}$| around |${{\phi }_2}$|:
The digital rock construction approach not only maintains geological relevance but also expands the application potential and range through geostatistical considerations, contributing to the statistical analysis of elasticity and attenuation characteristics in deep carbonate reservoirs featuring fracture-cavity systems. Additionally, the metre-scale models developed using this method have the capability to characterize large-scale REVs, addressing the limitations of traditional geophysical methods.
2.2 Numerical simulation based on Biot's poroelastic equation
We employ the quasi-static Biot's equations of poroelastic consolidation in the displacement–pressure (u − p) formulation to compute the equivalent properties featuring fractures, caverns and background pores (Biot 1941, 1962; Quintal et al. 2011). Disregarding body forces, these equations in the space-frequency domain can be expressed as follows:
where |$\nabla $| represents the Nabla operator, |$\omega $| represents the angular frequency and i denotes the imaginary unit. The parameters p, |$\kappa $| and |$\eta $| denote fluid pressure, permeability and fluid viscosity, respectively. |${\boldsymbol{u}}$| represents the solid displacement vector, with |${{u}_i}$| denoting the displacement of the solid in the i-th direction. |${\boldsymbol{\sigma }}$| signifies the total stress tensor, with its components |${{\sigma }_{ij}}$| expressed as (Biot 1941, 1962):
where |${{\delta }_{ij}}$| is the Kronecker delta. |$\lambda $| and |$\mu $| are the Lamé parameters. The component of the strain tensor |${{\varepsilon }_{ij}}$| and the trace of the strain tensor e can be expressed as:
The Biot's poroelastic parameters |$\alpha $| and M are defined as:
Here, |${{K}_s}$|, |${{K}_m}$| and |${{K}_f}$| denote the bulk moduli of the solid grains, dry frame and fluids, respectively, with |$\varphi $| representing porosity.
We use a finite element approach based on the discrete weak form of eqs (4), accounting for natural undrained boundaries. The weak form, detailed in appendix B of Quintal et al. (2011), is directly implemented in the finite element software, COMSOL Multiphysics, where displacement boundary conditions for relaxation tests are applied to the model to solve the governing equations.
For a compressional relaxation test in the vertical direction, a harmonic solid displacement is imposed perpendicular to the upper boundary of the model, while the solid displacement perpendicular to the other three boundaries is held at zero (Fig. 3a). Horizontal tests are performed by rotating the boundary conditions by 90° (Fig. 3b). The prescribed boundary conditions embody approximately the periodic assumption which is essential for scaling up, implying an infinitely long repetition of the numerical model considered to be a representative elementary volume (REV).
Schematic diagram of boundary conditions for (a) vertical and (b) horizontal relaxation tests.
The numerical solution of the aforementioned equations, along with the described boundary conditions, yields the variables |${\boldsymbol{u}}$| and p as functions of frequency, and consequently the stress and strain fields [equations (5) and (6)]. Averaging the stress and strain fields over the spatial domain of the numerical model, allows for computing the P-wave modulus (H) in the vertical direction along with the corresponding P-wave attenuation (O'connell & Budiansky 1978; Quintal et al. 2011):
where |$\langle \cdot \rangle $| denotes the average of variables, while |${\rm{Im}}$| and |${\rm{Re}}$| signify the imaginary and real parts, respectively. Finally, the complex P-wave velocity can be formulated as follows:
where |$\rho $| represents the weighted density of the model.
Although typically applied to lower frequency problems, Biot's quasi-static equations are appropriate for the broad-band range in our simulations, as strong attenuation occurs in the 0.1–10 Hz range, satisfying its low-frequency requirements. The highest simulation frequency, 103 Hz, lies within a range where viscous and inertial effects are negligible. Moreover, prior studies confirm the accuracy and stability of quasi-static finite element simulations across broad frequency ranges (Quintal et al. 2011; 2019).
It is worth mentioning that the current study is limited to the 2-D numerical simulation, as the 3-D connectivity and morphology of fractures and cavities could influence dispersion and attenuation characteristics. However, constructing digital rock models of 3-D fracture-cavity systems is inherently more complex, requiring methodologies beyond the use of outcrop images employed in 2-D modelling. Another challenge of future research on the simulations of seismic wave dispersion and attenuation in 3-D digital rock models is the computational demands, which may need the development of more efficient algorithms or the deployment of more powerful hardware. Future research regarding numerical simulation on 3-D digital rock models is important to overcome these challenges and advance the understanding of wave dispersion and attenuation signatures of carbonate fracture-cavity systems.
2.3 Parameter setting for the numerical simulations
The solid and fluid parameters for the numerical simulations are presented in Table 2. A porosity of 0.99 is assigned to the near-void fractures and cavities (Quintal et al. 2019), while the bulk and shear moduli of the rock skeleton are set to 5 |$\ \times $| 106 Pa after verifying the numerical stability of the model.
| Background | Fractures and caves | |
|---|---|---|
| Mineral bulk modulus (GPa) | 76.8 | 76.8 |
| Mineral shear modulus (GPa) | 32 | 32 |
| Mineral density (kg m−3) | 2710 | 2710 |
| Permeability (m−2) | 5 |$\ \times $| 10–15 | 10–11 |
| Porosity (per cent) | 3 | 99 |
| Dry bulk modulus (GPa) | 70.6 | 0.005 |
| Dry shear modulus (GPa) | 30.2 | 0.005 |
| Properties of fluid | ||
| Bulk Modulus (GPa) | 2.2 | |
| Density (kg m−3) | 1010 | |
| Viscosity (cP) | 1 | |
| Background | Fractures and caves | |
|---|---|---|
| Mineral bulk modulus (GPa) | 76.8 | 76.8 |
| Mineral shear modulus (GPa) | 32 | 32 |
| Mineral density (kg m−3) | 2710 | 2710 |
| Permeability (m−2) | 5 |$\ \times $| 10–15 | 10–11 |
| Porosity (per cent) | 3 | 99 |
| Dry bulk modulus (GPa) | 70.6 | 0.005 |
| Dry shear modulus (GPa) | 30.2 | 0.005 |
| Properties of fluid | ||
| Bulk Modulus (GPa) | 2.2 | |
| Density (kg m−3) | 1010 | |
| Viscosity (cP) | 1 | |
| Background | Fractures and caves | |
|---|---|---|
| Mineral bulk modulus (GPa) | 76.8 | 76.8 |
| Mineral shear modulus (GPa) | 32 | 32 |
| Mineral density (kg m−3) | 2710 | 2710 |
| Permeability (m−2) | 5 |$\ \times $| 10–15 | 10–11 |
| Porosity (per cent) | 3 | 99 |
| Dry bulk modulus (GPa) | 70.6 | 0.005 |
| Dry shear modulus (GPa) | 30.2 | 0.005 |
| Properties of fluid | ||
| Bulk Modulus (GPa) | 2.2 | |
| Density (kg m−3) | 1010 | |
| Viscosity (cP) | 1 | |
| Background | Fractures and caves | |
|---|---|---|
| Mineral bulk modulus (GPa) | 76.8 | 76.8 |
| Mineral shear modulus (GPa) | 32 | 32 |
| Mineral density (kg m−3) | 2710 | 2710 |
| Permeability (m−2) | 5 |$\ \times $| 10–15 | 10–11 |
| Porosity (per cent) | 3 | 99 |
| Dry bulk modulus (GPa) | 70.6 | 0.005 |
| Dry shear modulus (GPa) | 30.2 | 0.005 |
| Properties of fluid | ||
| Bulk Modulus (GPa) | 2.2 | |
| Density (kg m−3) | 1010 | |
| Viscosity (cP) | 1 | |
It should be emphasized that when employing Biot's theory to do the simulation, it is crucial to be cautious about the properties of rock skeletons inside fractures and cavities. An excessively low skeleton's bulk modulus (dry bulk modulus) can induce numerical instability and yield anomalous results, whereas a high modulus may not reflect the near-void conditions of fractures and cavities. Hence, it is essential to conduct tests prior to simulation endeavours. Quintal et al. (2019) have corroborated that when setting the dry bulk modulus to 106 Pa, the simulation outcomes when using Biot's equations align with those of the coupled Lame–Navier and Navier–Stokes equations. Similarly, He et al. (2022) also validated that the Biot's equation is applicable for describing fractures within the dry bulk modulus range of approximately 3 |$\ \times $| 104 to 3 |$\ \times $| 106 Pa.
The numerical model is discretized using an unstructured triangular mesh, with the 0.8 m fracture divided into approximately 50 elements, and a finer mesh applied near the fracture intersections, where most of the fluid flow occurs. The meshes consist of approximately 210 000 domain elements and 13 000 boundary elements, ensuring a high resolution that accurately captures the wave-induced fluid flow generated by the complex fracture-cavity system. The problem is solved in the frequency domain using the direct solver MUMPS (multifrontal massively parallel sparse direct solver), with each frequency taking approximately 45 s to compute.
3 RESULTS
The objective of our work is to study the influence of karstified cavities and fractures on the attenuation properties of deep carbonate reservoirs. More precisely, we aim at understanding the effect of the shape, size and number of cavities, as well as the number, length, width and angle of fractures. In this section, we will analyse the effects of cavities and fractures on the anisotropic velocity and attenuation characteristics of carbonate reservoirs based on the results of numerical simulations.
3.1 Influences of cavity characteristics
In this section, we will comprehensively analyse how cavities characteristics affect wave dispersion and attenuation signatures, while keeping the number and locations of fractures in the two zones unchanged.
3.1.1 Cavity shape
Fig. 4 presents the simulation results for various cavity geometries. In addition to the irregular karstified cavity model, which was based on actual outcrops, we also constructed an elliptical cavity with the same aspect ratio and area, as well as a circle cavity with the same area.
Simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions for (c) cavity, (d) ellipse and (e) circle models with the same area.
The results indicate that the characteristic frequencies of metre-scale fracture-cavity systems typically occur at lower frequency bands than those of seismic exploration, highlighting the importance of employing broad-band geophysical detection methods. The influences of cavity shape on P-wave velocity and attenuation are minimal, particularly in the vertical direction, where the velocity and attenuation across different cavity shape models are nearly identical. In the horizontal direction, the circular cavity exhibits the highest velocity and strongest attenuation, followed by the model with the elliptical cavity. This is mainly because the circular cavity is much less sensitive to stress variations. The attenuation amplitude gradually increases as the cavity becomes more regular—from irregular karstified cavity to ellipse to circle—while the fracture network remains unchanged. One explanation is that the circular cavity exhibits the largest area-to-volume ratio (or contour-to-area ratio in two dimensions), followed by elliptical and irregular shapes. This increased ratio promotes more significant fluid pressure relaxation with the surrounding medium in the model, leading to stronger attenuation. Alternatively, the circular cavity may enhance the elastic modulus in certain regions, creating a greater contrast with the adjacent fractured regions, which increases the pore pressure relaxation between these areas.
Through the distributions from equations (1)–(3), we can systematically generate digital models with the same statistical parameters, which allows us to perform extensive simulations and obtain a comprehensive analysis of the impact of cavity shape. Fig. 5 illustrates the statistical results for 10 models, each with identical statistical parameters corresponding to the models depicted in Figs 4(c)–(e). The results encompass P-wave velocity and attenuation within the exploration seismic frequency band (1–100 Hz), along with the median cumulative attenuation values across the entire frequency band of the 10 models.
Statistical results for P-wave (a) velocity, (b) attenuation in the seismic frequency band (1–100 Hz) and (c) cumulative attenuation for cavity, ellipse and circle models with the same area.
In the vertical direction, the influence of cavity shape on velocity and attenuation is minimal even with the perturbations in the geometric properties of the models. This characteristic is evident in the analysis of the single model presented in Fig. 4 and is also reflected in the statistical results. The horizontal results indicate that the circle cavity model exhibits higher velocity and greater attenuation within the seismic frequency band, a trend that aligns with the analysis of the single model in Fig. 4. Additionally, the circular karst cavity model exhibits the highest cumulative attenuation across the full frequency band, suggesting the trends to dissipate more energy, which may hold significant implications for understanding the dissipation of seismic wave energy in deep carbonate reservoir exploration.
3.1.2 Cavity size
Fig. 6 illustrates the impact of cavity size on dispersion and attenuation. First, the size of the cavity exhibits a significant impact on velocity, with a larger cavity resulting in reduced velocities in both vertical and horizontal directions, as we would expect.
Simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions for the models with (c) small, (d) medium and (e) large cavity.
The attenuations associated with different cavity sizes are similar, because the attenuation is primarily governed by the complex fracture network within the fracture-cavity system, which partially obscures the influence of the cavities. However, the influence of cavity size on the characteristic frequency is noticeable, particularly in the horizontal direction, with a larger cavity causing the frequency to shift to a lower band. This is likely because the larger cavity corresponds to a longer equivalent characteristic length of the heterogeneities, which results in a longer relaxation time and a lower characteristic frequency.
Fig. 7 presents the statistical result for 10 models of different cavity sizes on dispersion and attenuation within the exploration seismic frequency band (1–100 Hz) and illustrates the changes in total cumulative attenuation. The impact of increasing size on velocity reduction is depicted in Fig. 6 and is further confirmed by the statistical results for the seismic frequency band shown in Fig. 7(a). The attenuation, primarily controlled by the dense fracture network, makes the increasing cavity size have a minimal impact on attenuation within the exploration seismic frequency band (Fig. 7b) and on cumulative attenuation (Fig. 7c), which illustrates the difficulty of determining the cavity size based solely on attenuation.
Statistical results for P-wave (a) velocity, (b) attenuation in the seismic frequency band (1–100 Hz) and (c) cumulative attenuation for different cavity sizes.
3.1.3 Cavity number
For a given area, one large cavity and several small cavities can lead to differences in rock properties, and Fig. 8 illustrates this effect. To simplify the analysis, the cavity shape is assumed to be circular. A row of multiple small cavities creates a weakened area, resulting in lower velocity compared to a single large cavity, especially in the horizontal direction (Fig. 8a). This effect is less pronounced in the vertical direction. An increase in the number of cavities will broaden the attenuation and cause the characteristic frequency to shift towards higher frequencies, probably because a greater number of smaller cavities have more contact with the fractures in the background compared to a single large cavity, leading to quicker pore pressure relaxation.
Simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal direction in the models with (c) one, (d) three and (e) eight cavities.
Fig. 9 shows the statistical results regarding the influence of the cavity number for 10 models with identical statistical parameters. The horizontal velocity in the seismic frequency band decreases with an increasing number of cavities because multiple small cavities arranged in a row reduce the equivalent horizontal strength of the rock. However, the effect on the vertical direction is minimal. The influence of cavity number on attenuation is complex. In the seismic frequency band (Fig. 9b), the models with three cavities exhibit the strongest attenuation, as their characteristic frequencies fall right within the seismic frequency range, which also results in the highest cumulative attenuation within the exploration seismic frequency band (Fig. 9c). An increase in cavity number reduces cumulative attenuation, particularly in the horizontal direction (Fig. 9c). This is due to the presence of multiple small cavities, which leads to more uniform properties between the slip fault zone (Zone 1) and the induced fracture zone (Zone 2) in the horizontal direction, thereby reducing attenuation between these regions.
Statistical results for P-wave (a) velocity, (b) attenuation in the seismic frequency band (1–100 Hz) and (c) cumulative attenuation for different cavity numbers.
3.2 Influence of fractures
In this section, we will primarily analyse the influence of fracture characteristics including fracture number, length, width and angle on the wave dispersion and attenuation signatures.
3.2.1 Fracture number
Figs 10 and 11 illustrate the impact of increasing fracture density on the equivalent properties of the carbonate reservoir containing fracture and cavities, where the colour blocks represent the range of velocities and attenuation of all 10 models generated based on the same statistical parameters.
Statistical simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions when the numbers of fractures in Zone 1 are 15, 25 and 50, while the fracture number in Zone 2 is 10. The three models displayed in Figs 12(d)–(f) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Statistical simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions when the numbers of fractures in Zone 2 are 10, 15 and 25, while the fracture number in Zone 1 is 50. The three models displayed in Fig. 12 (f)–(h) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
As shown in Fig. 10, when the number of fractures in Zone 2 is held constant at 10, increasing the number of fractures in Zone 1 leads to a significant reduction in velocity in both the horizontal and vertical directions. The characteristic frequency shifts towards higher frequencies as the fracture number increases. This is due to a higher number of interconnected fractures, which allow for greater fluid pressure diffusion within the fractures, leading to attenuation at higher frequencies compared to fluid pressure diffusion from the fractures into the surrounding background. Furthermore, an increased number of interconnections between fractures effectively reduces their characteristic lengths, causing the characteristic frequency to shift to a higher range, particularly when associated with fluid pressure diffusion within the fractures.
When the number of fractures in Zone 1 is fixed at 50 and the number of fractures in Zone 2 is increased, the velocity shows a slight decrease (Fig. 11a), while the change in attenuation remains minimal (Fig. 11b). This may be attributed to the geological constraints that naturally limit the number of fractures in Zone 2 compared to Zone 1, thereby reducing the influence of fractures in Zone 2 on the overall properties.
Fig. 12 illustrates the impact of the number of fractures on velocity and attenuation within the seismic frequency band, as well as cumulative attenuation across the entire frequency spectrum. As expected, an increase in the fracture number results in a decrease in velocity and an increase in attenuation (Fig. 12a). The effect of the fracture number in the slip fault zone (Zone 1) on both velocity and attenuation is more pronounced than that in the induced fracture zone (Zone 2). This is because the slip fault zone naturally contains denser fractures than the induced fracture zone, as it is closer to the centre of the fracture-cavity system and experiences more intense deformation (Micarelli et al. 2006). The statistical results indicate that the geophysical characterization of the slip fault zone can be effectively conducted using both velocity and attenuation, as both exhibit significant detectable changes. However, the characterization of the induced fracture zone should rely more on velocity changes, as attenuation is less influenced by these fractures.
Statistical analysis of the influence of fracture number on P-wave (a) velocity, and (b) attenuation at seismic frequency band (1–100 Hz) and (c) cumulative attenuation in the slip fault zone and induced fracture zone, respectively. The digital models from each of the five scenarios are shown in (d)–(h).
3.2.2 Fracture length
Fig. 13 illustrates the effect of fracture length on P-wave velocity dispersion and attenuation signatures. For digital rock models containing short, medium and long fractures, the half-lengths of fractures in slip fault zone and induced fracture zone are 0.4 and 0.8 m, 0.6 and 1.0 m, and 0.8 and 1.2 m, respectively (Li et al. 2016; Wu et al. 2019). The three models displayed in Figs 14(d), (g) and (j) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Statistical simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions in the models with short, medium and long fractures. The three models displayed in Figs 14(d), (g) and (j) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Statistical results for P-wave (a) velocity and (b) attenuation in the seismic frequency band (1–100 Hz) and (c) cumulative attenuation with different fracture lengths. The digital models from each of the seven scenarios are shown in (d)–(j).
An increase in fracture length leads to a decrease in velocity and an increase in attenuation because longer fractures effectively increase fracture density and the number of interconnections between fractures while the total number of fractures remains constant. The effect of increased fracture length on the characteristic frequency is mainly attributed to the fracture intersections which are more predominant for longer fractures (Hunziker et al. 2018). The intersections segment the fractures in the slip fault zone and the presence of shorter segments causes an increase in the characteristic frequency of attenuation caused by fluid pressure diffusion within the fractures.
Fig. 14 illustrates the effect of fracture length on velocity, attenuation within the seismic frequency band and cumulative attenuation. Generally, an increase in fracture length statistically reduces velocity and enhances attenuation within the seismic bands, both vertically and horizontally. And the fractures in the slip fault zone exhibit a more pronounced impact. The increase in fracture length also enhances cumulative attenuation (Fig. 14c), but this enhancement gradually levels off, indicating that the fractures are already sufficiently complex and intersected.
3.2.3 Fracture width
Fig. 15 illustrates the effect of fracture width on dispersion and attenuation, with the half-widths of fractures being 1, 2 and 4 cm for digital rock models containing narrow, medium and wide fractures, respectively. The three models displayed in Figs 16(d), (g) and (j) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Statistical simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions in the models with narrow, medium and wide fractures. The three models displayed in Figs 16(d), (g) and (j) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Statistical results for P-wave (a) velocity, and (b) attenuation at seismic frequency band (1–100 Hz) and (c) cumulative attenuation under different fracture widths. The digital models from each of the seven scenarios are shown in (d)–(j).
As fractures in the slip fault zone and induced fracture zone widen simultaneously, both velocity and attenuation decrease. An increase in fracture width reduces the equivalent bulk modulus, leading to a decrease in velocity (Fig. 15a). On the other hand, narrow fractures are more prone to stress accumulation, resulting in a larger pore pressure difference between the fractures and the embedding background, which leads to stronger attenuation. In contrast, the stress within wide fractures is distributed over a larger volume, reducing stress concentration and resulting in less pronounced attenuation. Meanwhile, changes in characteristic frequency are minimal. First, the fracture permeability is artificially set to be independent of fracture width (Quintal et al. 2019); while the background permeability is unaffected by fracture width, resulting in minimal impact of fracture width on the characteristic frequency. Secondly, variations in fracture width have minimal impact on the characteristic length of the fractures, and these changes do not affect the interconnection between fractures, which explains their limited effect on the attenuation pattern.
Statistical results reveal a similar trend, where velocity and attenuation decrease with increasing fracture width (Fig. 16). The fracture width in the slip fault zone (Zone 1) has a stronger influence, highlighting the significant role of the slip fault zone in carbonate reservoirs with fracture-cavity systems. Additionally, the statistical results suggest that in the geophysical characterization of fracture-cavern systems, especially in the context of oil and gas production, attention should not be limited to areas of high attenuation. Wide fracture areas with low attenuation may also be significant, as they could indicate potential oil and gas reservoirs.
3.2.4 Fracture angles
As illustrated in Fig. 17, the angles of conjugate fractures in the slip fault zone (Zone 1) have a significant impact on both P-wave dispersion and attenuation. In the vertical direction, as the fractures become more vertical, the velocity increases because the fractures align more closely with the direction of seismic wave propagation, resulting in a mechanically stiffer rock. In contrast, in the horizontal direction, a large number of fractures perpendicular to the seismic wave direction reduces the wave velocity.
Statistical simulation results for P-wave (a) velocity and (b) attenuation in vertical and horizontal directions when the fractures angles in Zone 1 are approximately 45|$^\circ $|, 70|$^\circ $| and 80|$^\circ $|. The three models displayed in Figs 18(d)–(f) represent one model from each of the three scenarios, with each scenario consisting of 10 models.
Regarding attenuation, more vertical fractures result in stronger attenuation in the horizontal direction, while the corresponding attenuation in the vertical direction is weaker. This occurs because a large number of fractures intersecting the seismic wave propagation at a high angle can generate stronger pore pressure relaxation with the surrounding background; however, the relaxation between intersecting fractures is highly complex and can behave oppositely. Similar behaviour was observed in the simulations by Quintal et al. (2014). Additionally, an anomalous attenuation peak appears in a very low-frequency range (around 10−3 Hz) for more vertical fractures (fractures angle at around 80|$^\circ $| in Fig. 17b), likely because a large number of near-vertical fractures lead to a coupled mechanical field similar to that of one long vertical fracture, equivalent to an increase in equivalent length, resulting in the lower frequency attenuation.
Statistical results indicate that fracture angles have a significant impact on the anisotropy of velocity and attenuation, as shown in Fig. 18. As the fractures become more vertically aligned, the vertical stiffness of the model increases significantly, which limits wave-induced fluid flow between the fractures and the surrounding matrix, thereby increasing vertical P-wave velocity and reducing vertical attenuation. Conversely, the changes are the opposite in the horizontal direction: velocity decreases and attenuation increases. Combining velocity anisotropy with attenuation anisotropy can reveal the dominant orientation of fractures, which is valuable for geophysical characterization.
Statistical results for P-wave (a) velocity, and (b) attenuation at seismic frequency band (1–100 Hz) and (c) cumulative attenuation under different fracture angles. The digital models from each of the three scenarios are shown in (d)–(f).
4 DISCUSSION
4.1 Mechanism of wave dispersion and attenuation in carbonates with complex fracture–cavity system
It is essential to understand how attenuation patterns might be affected by the cavities and fractures in the model. Fig. 19 illustrates a complex model with the fracture–cavity system (Fig. 19a), incorporating simplified models to illustrate the pore pressure diffusion mechanisms that exist within the system (Figs 19b–f). The mechanisms can be decomposed and analysed using simple models: a single cavity (Fig. 19b), a single fracture (Fig. 19c), intersecting fractures (Fig. 19d), a cavity with an unconnected fracture (Fig. 19e) and a cavity with connected intersecting fractures (Fig. 19f).
(a) Complex model of the fracture–cavity system. The simplified models include (b) one elliptical cavity, (c) one fracture, (d) intersecting fractures, (e) one elliptical cavity with a connected fracture and (f) one elliptical cavity with intersecting fractures, used to illustrate the pore pressure relaxation mechanisms within the system.
Fig. 20 illustrates the simulated attenuation results of the simple models shown in Figs 19(b)–(f) in responses to seismic waves at horizon directions. The models with only one type of heterogeneity (cavity or fracture) each exhibit a distinct characteristic frequency (models b and c), with pore pressure diffusion primarily occurring between them and the embedding background. The cavity model (model b) displays the smallest attenuation amplitude and the lowest characteristic frequency, whereas the fracture model (model c) exhibits both a higher characteristic frequency and greater attenuation amplitude. For models with multiple coexisting heterogeneities, the attenuation curves are more complex. Progressing from high to low frequency, pore pressure relaxation occurs between intersecting fractures, between the fracture and the cavity, between the fracture and the background and between the cavity and the background, as indicated in Fig. 20(b). These mechanisms may also interact with each other. For example, the presence of intersecting fractures contributes to the pressure diffusion within the fractures, resulting in differences in attenuation between models e and f. Model f exhibits two characteristic frequencies in the high-frequency range: one for the attenuation between the intersecting fractures and the other for the attenuation between the fractures and the cavity; while model e has only one characteristic frequency in the high-frequency range, corresponding to the attenuation between the fractures and the cavity. It is worth noting that models with complex fracture–cavity systems will exhibit more intricate and interconnected attenuation processes. Especially when the fracture number is large, the attenuations have different characteristic frequencies coupled with one another, forming a broad-spectrum continuous attenuation. Although the underlying attenuation mechanisms are the same, the attenuation in complex fracture-cavity systems can differ from the clear, isolated attenuation patterns observed in this section.
The simulated (a) P-wave velocity and (b) attenuation results of the models b–f shown in Fig. 19. And the pore pressure distribution at (c) f1 ≈ 0.0056 Hz, (d) f2 ≈ 1.7783 Hz and (e) f3 ≈ 17.7828 Hz of model f.
For the attenuation mechanism of model f incorporates all of the previously mentioned mechanisms simultaneously, we use model f as an example, and illustrate the pore pressure distributions at its three characteristic frequencies f1 (≈0.0056 Hz), f2 (≈1.7783 Hz) and f3 (≈17.7828 Hz), as shown in Figs 20(c)–(e). These results show that the relaxation process of pore pressure occurs initially between the intersecting fractures of high permeability, then diffuses to the connected caverns and ultimately spreads throughout the background region. At high frequency f3, pore pressure accumulates in vertical fracture, leading to a larger pore pressure differential compared to other areas (Fig. 20e). At frequency f2, the pore pressure in two intersecting fractures is relatively equilibrated, but it remains high compared to the cavity and the background area (Fig. 20d). At the lowest frequency f1, the pore pressure difference predominantly occurs between the heterogeneous bodies and the background pores, exhibiting lower amplitude and a broader range (Fig. 20c). Additionally, the attenuation between the intersecting fractures and the background does not exhibit a distinct characteristic frequency, because the fractures and cavity function as a cohesive unit, and perform pressure relaxation with the background (Fig. 20c).
4.2 Implications for geophysical characterization at broad frequency band
Based on the simulation results presented in Section 3, we observe that frequency dispersion and attenuation in deep carbonate reservoirs with complex fracture-cave systems are notably broad, with the primary attenuation typically occurring between approximately 10−2 to 101 Hz, which is slightly lower than the typical exploration seismic frequency band (e.g. 1–100 Hz). To enhance our understanding of the implications of broad-band attenuation for geophysical characterization, we calculate the mean attenuation across various frequency bands, as illustrated in Fig. 21.
Mean attenuation in different frequency bands.
The results indicate that attenuation can be detected across a range of five orders of magnitude, from approximately 10−3 to 102 Hz, which is attributed to the multiple mechanisms of pore pressure diffusion within the fracture-cave systems. Additionally, the highest attenuations are observed within the frequency band of approximately 10−1 to 100 Hz, which is lower than the traditional active-source exploration seismic frequency band. The broad- and low-frequency attenuation characteristics of deep, large-scale carbonate reservoirs with fracture-cavity systems necessitate the integration of diverse seismic detection methods, suggesting significant potential for geophysical characterization through a combination of active and passive sources. In fact, combining active seismic methods with passive sources such as ambient noise or natural earthquake recordings can enhance subsurface imaging, especially in regions where active source data are limited or difficult to acquire. Several studies have demonstrated the effectiveness of using natural earthquake attenuation to infer subsurface properties (Maxwell & Urbancic 2002; Saenger et al. 2009).
Within the exploration seismic frequency, small and numerous regular cavities are more easily detectable [sections 3.1(1) and 3.1(3)], while significant attenuation in this band may also indicate a complex fracture network, with fractures more likely to be long [section 3.2(1)) and narrow (section 3.2(3)]. The significant attenuation in earthquake seismology is more likely associated with large karst cavities [section 3.1(2)] and sparse, short fractures [sections 3.2(1) and 3.2(1)]. It should be noted that geometric parameters and rock physical properties, such as permeability, porosity and other characteristics of the fracture-cavity system, might influence the frequency distribution of attenuation, necessitating specific parameters for different reservoirs.
In addition, deep carbonate reservoirs with fracture-cavity systems exhibit pronounced attenuation anisotropy, with horizontal attenuations reaching approximately 5 to 6 times greater than vertical attenuations. This observation provides the insight that to fully leverage attenuation information for enhanced geophysical characterization, the acquisition of large-offset data is essential.
4.3 Rock physics template for seismic interpretation
Fig. 22 illustrates the relationship between attenuation and P-wave velocity at 10 Hz, with specific models featuring different fracture and cavity properties shown on the right side. The colour denotes the total fracture density, which is calculated using the fracture area proportions as well as the lengths and widths of fractures. The label size qualitatively represents the total porosity, encompassing fractures, cavities and background pores. The arrows indicate the increasing direction of the parameter.
Rock physics templates at (a) 10 Hz, with the specific models have (b) narrow fractures, (c) long fractures, (d) big cavity and (e) with fractures. The label colour indicates the fracture density, while the label size qualitatively represents the total porosity. The arrows indicate the increasing direction of the indicated parameter.
As illustrated in the template for the exploration seismic frequency band at 10 Hz (Fig. 22a), the increase in fracture density (resulting from the increase in the number and length of fractures) leads to a decrease in velocity and an increase in attenuation, indicating that velocity and attenuation can serve as important indicators of fracture density. When fracture density is similar, an increase in fracture width leads to higher porosity, thereby reducing velocity ((Figs 22b and e). However, attenuation increases as fracture width decreases because narrower fractures are more conducive to stress accumulation, resulting in a greater pore pressure difference. Please see the detailed analysis in Section 3.2(3). The impact of cavities is mainly reflected in velocity, where an increase in cavity size reduces velocity, with only a minimal effect on attenuation (Fig. 22d).
In summary, low P-wave velocity and high attenuation generally indicate high fracture density. Low-velocity reservoirs are more likely to contain larger cavities. Reservoirs with high porosity but low attenuation are likely associated with fractures having large widths, while those with low porosity but high attenuation may be due to fractures having small widths. Understanding the relationships depicted by the templates is crucial for interpreting seismic data and optimizing reservoir exploration strategies in heterogeneous carbonate reservoirs with fracture–cavity systems.
5 CONCLUSIONS
We perform numerical simulations on 2-D large-scale digital rock models to evaluate the velocity and attenuation characteristics of deep carbonate rocks with fracture–cavity systems. We develop a 2-D digital rock model construction method incorporating the geological outcrop constraints with statistical considerations. Through finite element simulations based on Biot's equations of poroelastic consolidation, our study reveals the intricate interplay between cavities and fractures in influencing seismic wave behaviour for deep carbonate reservoirs.
Deep carbonate rocks with heterogeneous fracture–cavity systems exhibit broad-band attenuation, potentially spanning over three orders of magnitude (from lower than 0.01 Hz to higher than 10 Hz), due to their coupled pore pressure relaxation mechanisms, highlighting the necessity of employing broad-band seismic exploration approaches to characterize those rocks in the subsurface. Horizontal attenuation is approximately five times greater than vertical attenuation due to the predominantly semivertical orientation of fractures, emphasizing the need for large-offset seismic data acquisition to accurately capture anisotropic effects in geophysical characterization.
Although small and numerous cavities can enhance attenuation in the seismic exploration frequency band—especially in the horizontal direction—the shape, size and number of cavities have limited influence within fracture–cavity systems, where attenuation is primarily governed by the surrounding fracture network. Lower velocities and higher attenuations are often associated with greater fracture density, longer fracture lengths or narrower fractures, with those in the slip fault zone having more pronounced impacts. We summarized our observations in a rock physics template which demonstrates the importance of the joint use of velocity and attenuation information to improve seismic interpretation and optimize reservoir characterization in complex deep carbonate systems.
ACKNOWLEDGMENTS
This work was supported by National Natural Science Foundation of China (grant no. 42330805 and U23B20157), Fundamental Research Funds for the Central Universities, and China Scholarship Council program.
DATA AVAILABILITY
Part of the data set and codes associated with this paper are available online at https://zenodo.org/records/14224435.
