Elastic anisotropy differentiation of thin shale beds and fractures using a novel hybrid rock physics model

Elastic anisotropy is frequently used to characterize fracture distribution. However, sets of parallel horizontal fractures and thin shale beds in tight sand can both cause vertical transverse isotropy. Here, we are not referring to shale layers on the logging scale but rather to very thin shale beds, a few centimeters thick, within tight sand. To accurately differentiate the anisotropy caused by horizontal fractures or thin shale beds, we propose a hybrid rock physics model. This new model combines the Hudson model and the shale-compacting orientation distribution function (ODF) model, based on the anisotropic self-consistent approximation (SCA) and differential effective medium (DEM) theory. The new model's reliability is demonstrated by comparison to the well logs. The proposed model can characterize the elastic properties of both thin shale beds and horizontal fractures. Based on this model, the rock physical analysis reveals that thin shale beds and horizontal fractures exhibit distinct elastic anisotropy characteristics. Furthermore, we analyze the seismic response differences between horizontal fractures and thin shale beds using the anisotropic Ruger's approximation formula. The analysis indicates that the seismic response of tight sand containing thin shale beds interferes with the fracture's identification. On the other hand, there are identifiable differences between the fractured tight sand and the tight sand containing thin shale beds. Based on this difference, we develop a new seismic attribute to characterize the fracture distribution. These difference-based attributes can effectively eliminate the interference from thin shale beds, making the distribution of horizontal fractures more apparent.

1 Introduction

With the continuous growth of global energy demand and the gradual depletion of conventional oil and gas resources, the development of unconventional oil and gas resources has become increasingly important (Gharavi et al., 2023). Among these unconventional resources, fracture reservoirs have become a research focus due to their excellent storage and flow capacity (Zhang et al., 2022). The presence and development of fractures affect the permeability and porosity of reservoirs and are directly related to the occurrence and flow of oil and gas (Liu et al., 2019; Jiang et al., 2023). Therefore, accurately predicting fracture distribution is crucial for oil and gas exploration and development. In recent years, fracture prediction techniques based on seismic data have made significant progress, becoming an important tool for researching fractures containing gas and oil (Liu et al., 2018; Li et al., 2020). For example, the application of seismic bright spot technology has proven effective in practice for the prediction of fractures containing gas (Fawad et al., 2020). The main idea of this technology is that fractures exhibit a significant low impedance contrast with the adjacent rock formations. However, thin shale beds present in tight sand may also exhibit the same low impedance properties. These two situations can cause similar seismic responses, which can mislead the characterization of fracture distribution. Therefore, it is crucial to differentiate the elastic properties caused by the fractures or thin shale beds (Lin et al., 2022). It is worth noting that the thin shale beds we studied are very thin, with thicknesses of a few centimeters, which are far below the resolution of well logging. Therefore, it is difficult to describe these thin shale beds within tight sand using seismic and logging techniques. Consequently, we need to rely on rock physics elastic models to equivalently represent their microstructure and convert it into macroscopic responses at seismic and logging scales.

The rock physics elastic modeling process is categorized into three primary components: matrix, skeleton, and fluid (Mavko et al., 2009). The matrix model represents the amalgamation of the diverse minerals found in a rock based on their composition. For a homogeneous mineral matrix, various averaging methods can be used to synthesize an isotropic rock matrix. Voigt (1890) proposed the equivalent strain-averaging model. Reuss (1929) proposed the equivalent stress-averaging model. These models give the theoretical elastic parameter range of the rock. Hill (1952) gave the elastic parameters by averaging the upper and lower bounds. Wyllie and Gregory (1953) proposed a linear formula so that, when the rock has uniformly distributed intergranular pores, there is a linear relationship between porosity and acoustic transit time. Hashin and Shtrikman (1963) gave the lower bound of elastic combination parameters for the softest rock and the upper bound for the hardest rock of the mineral composition. By using the above methods, the mineral component equivalent medium is treated as the rock matrix (Alabbad et al., 2023).

The skeleton models are rock structure models with inhomogeneous phases inserted into the matrix background phase (Ma et al., 2024). Kuster and Toksöz (1974) gave a skeleton-equivalent model for different pore shapes in carbonate. Hudson (1980, 1981) proposed a flat coin-shaped equivalent crack model. Schoenberg and Sayers (1995) gave a linear-slip model. Xu and White (1995, 1996) combined the self-consistent model to optimize the Kuster and Toksöz model. Xu and Payne (2009) further gave an equivalent pore model for carbonate. Chapman (2003) gave a multi-scale equivalent fracture model. Lian et al. (2024) combined experimental measurement results to obtain the compacted optimized skeleton model using compaction coefficients.

The fluid models describe the elastic and anisotropic characteristics of various fluids in relation to the skeleton. For simple isotropic models, Gassmann (1951) gave the fluid replacement formula at low frequencies. Brown and Korringa (1975) proposed an anisotropic Gassmann formula (BK model) for anisotropic skeletons (Thomsen, 2023). Using a fluid substitution model, Guo et al. (2023) investigated the relationship between acoustic velocity and fluid saturation. However, the microstructure of thin shale beds is complicated (Zhou et al., 2023). Hence, it is necessary to give greater attention to the properties of elasticity and anisotropy. Therefore, thin shale beds should not be regarded as a component of a matrix model or as a straightforward skeletal inclusion phase. This research aims to integrate the anisotropic model of thin shale beds with a specific component of the rock skeleton. For the orientation of the thin shale beds, Roe (1965) defined the direction functions in the three-dimensional space of rock. These direction functions can be represented by the Legendre coefficients corresponding to a sequence of functions. For the thin shale beds by compaction, their normal is parallel to the third axis, so it can be equivalent to the vertical transverse isotropic (VTI) medium. For VTI media, Sayers (1995) simplified the Legendre coefficients into two, which can express the elastic stiffness matrix of the thin shale beds. On the other hand, Johansen et al. (2004) gave the self-consistent approximation (SCA) and differential effective medium (DEM) models to calculate the stiffness matrix of the single thin shale bed.

In this paper, we first propose a model to combine the thin shale beds and horizontal fracture skeleton. This model is verified by field data in a Sichuan Basin gas reservoir. The tight sand containing fractures and the thin shale beds are funded in this area (Ding et al., 2021; Yurikov et al., 2021). Based on the fracture orientations and dip angles within the area, we assume the fracture anisotropy to be VTI anisotropy, then we analyze the seismic response differences between fracture and thin shale beds using the anisotropic Ruger's approximation formula. Finally, based on the seismic response, we develop a new seismic attribute to explore the potential position of the horizontal fractures.

2 Methods

2.1 Rock physics modeling process

To study the VTI anisotropy of shale and horizontal fractures, we proposed the modeling workflow shown in Fig. 1. Vertical fractures, due to their orientation, exhibit horizontal transverse isotropic (HTI) anisotropy, which can be directly distinguished from shale without modeling. The modeling process is mainly divided into two parts: the fracture-skeleton and the thin-shale-bed rock physics modeling. Specifically, we calculated the stiffness matrix of the sand skeleton containing fractures and fluids in the sand model. On the other hand, we computed the stiffness matrix of the bed structure in the thin-shale-bed model. Finally, we combined the two stiffness matrices using the SCA and DEM models to obtain the rock stiffness matrix that includes horizontal fractures and thin shale beds. The methods of obtaining the parameters required for the whole modeling process and their meanings are shown in Table 1.

Table 1 The meanings and methods of the rock model parameters. The table outlines the methods for obtaining the parameters required in the technical process depicted in Fig. 1, along with their corresponding definitions. These fundamental parameters facilitate the utilization of rock physics elastic models, enabling the conversion between elastic properties and the other physical parameters.

Figure 1 The rock modeling flowchart. The modeling consists of three branches. The left branch simulates the elastic properties of shale using the orientation distribution function (ODF) and Hornby models. The middle branch focuses on the elastic properties of the fracture with the Hudson approach. The right branch simulates the elastic properties of fluids based on the Wood model.

2.1.1 The fracture-skeleton rock physics modeling process

For the fracture skeleton, we used the Voigt–Reuss–Hill (VRH) model (Hill, 1952) to build the sand matrix (Msand):

$$\begin{array}{}\text{(1)}& \mathrm{Msand}=\left({\sum }_{i=\mathrm{1}}^N{f}_i{M}_i+{\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^N\frac{f_i}{M_i}}}\right)/\mathrm{2},\end{array}$$

where f_i is the content of the ith mineral in the matrix and M_i is the isotropic modulus of the ith mineral in the matrix (bulk modulus K_i or shear modulus μ_i). Msand is the background isotropic modulus (bulk modulus K_m or shear modulus μ_m). The isotropic modulus can be converted into the sand isotropic stiffness of the matrix using the following equation:

$$\begin{array}{}\text{(2)}& C_{\mathrm{sandm}}=\left[\begin{array}{llllll}{K}_m+\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\mu }}_m& K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& K_m+\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\mu }}_m& K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& K_m-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_m& K_m+\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\mu }}_m& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{\mu }}_m& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{\mu }}_m& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{\mu }}_m\end{array}\right].\end{array}$$

The Hudson model (Hudson, 1980, 1981) is based on a scattering theory analysis of the mean wave field in an elastic solid with thin, penny-shaped ellipsoidal cracks or inclusions. Hudson proposed a Taylor expansion approximation to calculate the stiffness matrix for a fracture porosity composite system (Hudson et al., 1996):

$$\begin{array}{}\text{(3)}& C_{\mathrm{dsand}}=C_{\mathrm{sandm}}+C^{\mathrm{1}}+C^{\mathrm{2}},\end{array}$$

where C¹ is the first-order correction term for the anisotropy caused by the fracture and C² is the second-order correction term of the anisotropy caused by the mutual coupling between the directional fractures. Both C¹ and C² are calculated from the pore aspect ratio a, the rock matrix porosity ϕ_m, and the fracture porosity ϕ_f ($\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{m}}+{\mathit{\varphi }}_{\mathrm{f}}$). We assume that the percentage of fracture porosity to total porosity ϕ in the reservoir is a constant (3.2 %). These skeleton parameters were obtained from log curves and thin sections.

We introduced fluids into the sand skeleton using the BK model (Brown and Korringa, 1975) to obtain the saturated sand rock stiffness matrix:

$$\begin{array}{}\text{(4)}& s_{\mathrm{ijkl}}^{\mathrm{sand}}=s_{\mathrm{ijkl}}^{\mathrm{dsand}}-{\displaystyle \frac{\left(s_{\mathrm{ijaa}}^{\mathrm{sand}}-s_{\mathrm{ijaa}}^{\mathrm{sandm}}\right)\left(s_{\mathrm{bbkl}}^{\mathrm{sand}}-s_{\mathrm{bbkl}}^{\mathrm{sandm}}\right)}{\left(s_{\mathrm{ccdd}}^{\mathrm{sand}}-s_{\mathrm{ccdd}}^{\mathrm{sandm}}\right)-\mathit{\varphi }\left(\frac{\mathrm{1}}{K_{\mathrm{fl}}}-\frac{\mathrm{1}}{K_{\mathrm{sandm}}}\right)}},\end{array}$$

where the parameters $s_{\mathrm{ijkl}}^{\mathrm{dsand}}$ and $s_{\mathrm{ijkl}}^{\mathrm{sandm}}$ represent the flexibility of dry rock skeleton and rock matrix minerals, respectively. The stiffness matrix can be inverted from the flexibility matrix following C_ijklS_ijkl=I and vice versa. K_fl can be obtained by the Wood (1957) formula.

2.1.2 The thin-shale-bed rock physics modeling process

The thin shale beds in the tight sand can be regarded as being composed of clay domains (Bandyopadhyay, 2009). These clay domains exhibit laminar structures and strong anisotropy. To construct a single clay domain, we referenced Hornby's procedure (Hornby et al., 1995). Hornby's method effectively describes the complex structure of a single clay domain. After we obtain the anisotropy of a single clay domain, the equivalent elastic stiffness of shale beds' orientation is obtained by taking the Voigt average (Sayers, 1995):

$$\begin{array}{}\text{(5)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{11}}}=L+\mathrm{2}M+{\displaystyle \frac{\mathrm{4}\sqrt{\mathrm{2}}}{\mathrm{105}}}{\mathit{\pi }}^{\mathrm{2}}\left[\mathrm{2}\sqrt{\mathrm{5}}{a}_{\mathrm{3}}{W}_{\mathrm{200}}+\mathrm{3}{a}_{\mathrm{1}}{W}_{\mathrm{400}}\right],\text{(6)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{33}}}=L+\mathrm{2}M{\displaystyle \frac{\mathrm{16}\sqrt{\mathrm{2}}}{\mathrm{105}}}{\mathit{\pi }}^{\mathrm{2}}\left[\sqrt{\mathrm{5}}{a}_{\mathrm{3}}{W}_{\mathrm{200}}-\mathrm{2}{a}_{\mathrm{1}}{W}_{\mathrm{400}}\right],\text{(7)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{12}}}=L-{\displaystyle \frac{\mathrm{4}\sqrt{\mathrm{2}}}{\mathrm{315}}}{\mathit{\pi }}^{\mathrm{2}}\left[\mathrm{2}\sqrt{\mathrm{5}}\left(\mathrm{7}{a}_{\mathrm{2}}-a_{\mathrm{3}}\right)W_{\mathrm{200}}-\mathrm{3}{a}_{\mathrm{1}}{W}_{\mathrm{400}}\right],\text{(8)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{13}}}=L+{\displaystyle \frac{\mathrm{4}\sqrt{\mathrm{2}}}{\mathrm{315}}}{\mathit{\pi }}^{\mathrm{2}}\left[\sqrt{\mathrm{5}}\left(\mathrm{7}{a}_{\mathrm{2}}-a_{\mathrm{3}}\right)W_{\mathrm{200}}-\mathrm{12}{a}_{\mathrm{1}}{W}_{\mathrm{400}}\right],\text{(9)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{44}}}=M-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{2}}}{\mathrm{315}}}{\mathit{\pi }}^{\mathrm{2}}\left[\sqrt{\mathrm{5}}\left(\mathrm{7}{a}_{\mathrm{2}}+a_{\mathrm{3}}\right)W_{\mathrm{200}}+\mathrm{24}{a}_{\mathrm{1}}{W}_{\mathrm{400}}\right],\text{(10)}& {\displaystyle }{C}_{{\mathrm{shale}}_{\mathrm{66}}}={\displaystyle \frac{C_{{\mathrm{shale}}_{\mathrm{11}}}-C_{{\mathrm{shale}}_{\mathrm{12}}}}{\mathrm{2}}},\end{array}$$

where a₁, a₂, a₃, and L are the elastic parameters of a single clay domain ($a_{\mathrm{1}}=C_{{\mathrm{claym}}_{\mathrm{11}}}+C_{{\mathrm{claym}}_{\mathrm{33}}}-\mathrm{2}{C}_{{\mathrm{claym}}_{\mathrm{13}}}-\mathrm{4}{C}_{{\mathrm{claym}}_{\mathrm{44}}}$, $a_{\mathrm{2}}=C_{{\mathrm{claym}}_{\mathrm{11}}}-\mathrm{3}{C}_{{\mathrm{claym}}_{\mathrm{12}}}+\mathrm{2}{C}_{{\mathrm{claym}}_{\mathrm{13}}}-\mathrm{2}{C}_{{\mathrm{claym}}_{\mathrm{44}}}$, $a_{\mathrm{3}}=\mathrm{4}{C}_{{\mathrm{claym}}_{\mathrm{11}}}-\mathrm{3}{C}_{{\mathrm{claym}}_{\mathrm{33}}}-C_{{\mathrm{claym}}_{\mathrm{13}}}-\mathrm{2}{C}_{{\mathrm{claym}}_{\mathrm{44}}}$, $L=\frac{\mathrm{1}}{\mathrm{15}}\left(C_{{\mathrm{claym}}_{\mathrm{11}}}+C_{{\mathrm{claym}}_{\mathrm{33}}}+\mathrm{5}{C}_{{\mathrm{claym}}_{\mathrm{12}}}+\mathrm{8}{C}_{{\mathrm{claym}}_{\mathrm{13}}}-\mathrm{4}{C}_{{\mathrm{claym}}_{\mathrm{44}}}\right)$) and c_claym is the elastic matrix of a single clay domain. The coefficients W₂₀₀ and W₄₀₀ are the Legendre coefficients of the orientation distribution function (ODF) and can be obtained through X-ray diffraction (XRD) experiments on core samples. However, for shale, core measurement is challenging and expensive. Thus, we used the compaction distribution function W(ξ) derived by Johansen et al. (2004) to calculate these two parameters as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }{W}_{\mathrm{200}}=\sqrt{{\displaystyle \frac{\mathrm{5}}{\mathrm{2}}}}{\int }_{-\mathrm{1}}^{\mathrm{1}}W\left(\mathit{\xi }\right)P_{\mathrm{2}}\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi },\text{(12)}& {\displaystyle }{W}_{\mathrm{400}}=\sqrt{{\displaystyle \frac{\mathrm{9}}{\mathrm{2}}}}{\int }_{-\mathrm{1}}^{\mathrm{1}}W\left(\mathit{\xi }\right)P_{\mathrm{4}}\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi },\end{array}$$

where $P_{\mathrm{2}}\left(\mathit{\xi }\right)=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{3}{\mathit{\xi }}^{\mathrm{2}}-\mathrm{1}\right)$, $P_{\mathrm{4}}\left(\mathit{\xi }\right)=\frac{\mathrm{1}}{\mathrm{8}}\left(\mathrm{35}{\mathit{\xi }}^{\mathrm{4}}-\mathrm{30}{\mathit{\xi }}^{\mathrm{2}}+\mathrm{3}\right)$, and ξ=cos (θ), $W\left(\mathit{\xi }\right)=\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\frac{A^{\mathrm{2}}}{{\left({\mathit{\xi }}^{\mathrm{2}}+A^{\mathrm{2}}\left(\mathrm{1}-{\mathit{\xi }}^{\mathrm{2}}\right)\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}.$

The only parameter that needs to be inputted is A (compaction factor), which represents the ratio of shale thicknesses before and after compaction. It is worth noting that obtaining A is difficult. For more convenient applications, we will use the critical-porosity-based compacted ODF to characterize the bed orientation. Specifically, the relationship between A and the critical porosity is given by Bachrach (2011):

$$\begin{array}{}\text{(13)}& A={\left({\displaystyle \frac{\mathit{\varphi }}{{\mathit{\varphi }}_{\mathrm{0}}}}\right)}^k.\end{array}$$

The model parameter k controls the speed at which the pore space deforms (Bachrach, 2011). Thus, the value of k is calibrated based on the shale compaction curve, and, for normal compaction, k equals 1. In the entire process of shale modeling, the sole free variable, A, can be converted into the relationship between porosity and critical porosity. Critical porosity is defined as the point at which clay particles in suspension make contact, leading to a phase transition that results in a finite shear modulus.

2.1.3 Thin-shale-bed and fracture-skeleton hybrid rock physics model

The stiffness matrices of the thin shale beds and the fracture skeleton were calculated using the previous two sections, and the thin shale beds were inserted into the fracture skeleton model using the SCA and DEM models. To avoid having sand and shale as two isolated entities, the model constructs the structure in two steps: establishing a mixed medium with half sand and half shale, then randomly inserting parts with actual sand or shale content exceeding half. The specific steps are as follows.

Firstly, the anisotropic SCA model (Hornby et al., 1995) is used to construct a shale and sand bi-connected equivalent structure. The bi-connected equivalent structure is

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& C_{\text{shale-sand}}^{\mathrm{bio}}={\sum }_{n=\mathrm{1}}^{\mathrm{2}}\mathrm{0.5}{C}_n{\left(I+\widehat{G}\left(C_n-C^{\mathrm{bio}}\right)\right)}^{-\mathrm{1}}\\ & {\left\{{\sum }_{p=\mathrm{1}}^{\mathrm{2}}\mathrm{0.5}{\left(I+{\widehat{G}}_{ijkl}\left(C_n-C^{\mathrm{SCA}}\right)\right)}^{-\mathrm{1}}\right\}}^{-\mathrm{1}}.\end{array}\end{array}$$

Here, when n and p are equal to 1, it refers to the shale stiffness matrix, and, when they are equal to 2, it refers to the sand stiffness matrix; $C_{\text{shale-sand}}^{\mathrm{bio}}$ is the stiffness matrix of the shale and sand bi-connected equivalent structure. The anisotropic DEM model is then used to complete the remaining structure. ${\widehat{G}}_{\mathrm{ijkl}}$ represents the geometric parameters of the inclusions, and its calculation process is provided by Mura (2013). The insufficient components are added gradually in equal amounts to the shale and sand bi-connected equivalent structure until the actual shale content of the rock is reached. Thus, the calculation of the DEM model is an iterative process. Specifically, it involves subdividing the inclusions into n parts. With each addition of a part, the stiffness matrix of the background phase is updated. The result of the ith iteration is as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left(C^{\mathrm{DEM}}\left(v_i\right)\right)}{d{v}_i}}& ={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-v_i\right)}}\left(C^{\mathrm{inclusion}}\left(i\right)-C^{\mathrm{DEM}}\left(v_i\right)\right)\\ & {\left[\underset{\mathrm{¯}}{I}+{\widehat{G}}_{\mathrm{ijkl}}\left(C^{\mathrm{inclusion}}-C^{\mathrm{DEM}}\left(v_i\right)\right)\right]}^{-\mathrm{1}},\end{array}\text{(16)}& {\displaystyle }{C}^{\mathrm{DEM}}\left(v_{\mathrm{1}}\right)=C_{\text{shale-sand}}^{\mathrm{bio}},\text{(17)}& {\displaystyle }{C}^{\mathrm{DEM}}\left(v_{i+\mathrm{1}}\right)=C^{\mathrm{DEM}}\left(v_i\right)+\mathrm{d}\left(C^{\mathrm{DEM}}\left(v_i\right)\right),\text{(18)}& {\displaystyle }{C}^{\mathrm{DEM}}\left(v_n\right)=C^{\mathrm{rock}},\end{array}$$

where v_i and C^DEM(v_i) are, respectively, the corresponding inclusion volume content and the elastic matrix of the inclusion component before the ith insertion. C^inclusion is the stiffness matrix of the inclusion component, and C^rock is the stiffness matrix of the result.

The final rock stiffness matrix can be used to obtain the rock acoustic and anisotropic parameters through the velocity–elasticity relationship and the anisotropy model (Thomsen, 1986):

$$\begin{array}{}\text{(19)}& {\displaystyle }{V}_p=\sqrt{{\displaystyle \frac{c_{\mathrm{33}}^{\mathrm{rock}}}{\mathit{\rho }}}},\text{(20)}& {\displaystyle }{V}_{\mathrm{s}}=\sqrt{{\displaystyle \frac{c_{\mathrm{44}}^{\mathrm{rock}}}{\mathit{\rho }}}},\text{(21)}& {\displaystyle }\mathit{\epsilon }={\displaystyle \frac{c_{\mathrm{11}}^{\mathrm{rock}}-c_{\mathrm{33}}^{\mathrm{rock}}}{\mathrm{2}{c}_{\mathrm{33}}^{\mathrm{sat}}}},\text{(22)}& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{c_{\mathrm{66}}^{\mathrm{rcok}}-c_{\mathrm{44}}^{\mathrm{rock}}}{\mathrm{2}{c}_{\mathrm{44}}^{\mathrm{rock}}}},\text{(23)}& {\displaystyle }\mathit{\delta }={\displaystyle \frac{{\left(c_{\mathrm{13}}^{\mathrm{rock}}+c_{\mathrm{44}}^{\mathrm{rock}}\right)}^{\mathrm{2}}-{\left(c_{\mathrm{13}}^{\mathrm{rock}}-c_{\mathrm{44}}^{\mathrm{rock}}\right)}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{33}}^{\mathrm{rock}}\left(c_{\mathrm{33}}^{\mathrm{rock}}-c_{\mathrm{44}}^{\mathrm{rock}}\right)}},\end{array}$$

where ρ, V_p, and V_s are the equivalent density, the velocities of the longitudinal wave, and the shear wave for the rocks, respectively. ε, γ, and δ are the Thomsen anisotropy parameters, and the equivalent density can easily be calculated by Voigt's averaging.

3 Background and model calibration

3.1 Geological background

The research focuses on the tight sand gas in the Xujiahe Formation in the Sichuan Basin. The region has dense lithology, with fractures serving as the primary migration pathways (Huang et al., 2022). Horizontal fractures can help identify tight gas reservoirs in this area (Yue et al., 2018; Zhang, 2021). The developed horizontal fractures can also effectively assist in the water injection development of tight gas (Zhao et al., 2021). Previous Chinese scholars have analyzed horizontal fractures, which exhibit preferential distribution and structural features of VTI anisotropy (Su, 2011). The core samples and imaging logging showed horizontal fractures in this study region (Fig. 2). We also statistically analyzed the fracture distribution primarily based on the imaging logging and core samples.

Figure 2 The evidence for the horizontal fractures in the area. Panels (a) and (b) show cores with developed horizontal fractures; the red arrows point to the horizontal fractures. Panel (c) shows the imaging logging of the horizontal fracture formation.

Through the core samples, we have analyzed the width and length of these horizontal fractures (Fig. 3a, b). Most fractures have a width of less than 1 mm, with lengths primarily ranging from 100 to 300 mm. Therefore, we set the fracture aspect ratio a in the Hudson model to 0.01. We also calculated the orientation of these horizontal fractures though imaging logging (Fig. 3c). Most fractures were concentrated around the 60° direction. This revealed that these horizontal fractures exhibited VTI anisotropy. Therefore, we used the VTI equivalent model to study the anisotropy of fractures in the area. Han et al. (2022) analyzed the fracture VTI anisotropy in the area but did not consider the VTI anisotropy caused by thin shale beds.

Figure 3 Fracture parament statistics. Panels (a) and (b) show the statistics of fracture length and width, respectively. Most fractures have a width of less than 1 mm, with lengths primarily ranging from 100 to 300 mm. Panel (c) is a rose diagram of fracture orientations, showing that most fractures are concentrated around the 60° direction.

We conducted statistical analysis on the average fracture porosity and total porosity interpretation results from 10 wells in the study region. The results are shown in Fig. 4. The fracture porosity in the region ranges from 0.04 % to 0.39 %, accounting for 3.2 % of the total porosity.

Figure 4 Fracture porosity and matrix porosity analysis. The bar chart represents the average values of porosity and fracture porosity within the target interval from logs in the area. The line chart indicates the percentage of fracture porosity relative to total porosity.

The compliance effectively characterizes a material's ability to deform under stress, making it particularly suitable for representing fractures. The linear-slip deformation (LSD) model proposed by Schoenberg and Sayers (1995) assumes a linear relationship between the displacement discontinuity across a fracture surface and the applied stress, enabling the compliance contributions from multiple fractures to be directly summed. In contrast, Hudson's model treats fractures as small perturbations to the stiffness tensor, based on the assumption that fractures introduce only minor modifications to the medium. Specifically, it uses a first-order Taylor expansion to approximate the effects of fractures on the stiffness. This approach makes stiffness perturbation a more suitable and computationally efficient framework for modeling low fracture densities. Furthermore, stiffness-based models offer a more direct approach for analyzing seismic wave propagation, particularly when fracture porosity is sufficiently low. When fracture porosity exceeds 0.45 % (fracture density = 0.1), the Hudson model becomes less effective in describing the elastic properties of fractures (Fig. 5). In this study, the fractures exhibit a maximum porosity of 0.39 %, which is well within the descriptive capabilities of the Hudson model. Therefore, we adopted the Hudson model for our analysis.

Figure 5 The solid line represents the first-order Hudson formula, the dashed line represents the second-order Hudson formula, and the points indicate well log samples. The first-order results of the Hudson model consider the effects of fractures, while the second-order results account for both the effects of fractures and their interactions. In this study, we used the second-order results of the Hudson model.

The area is characterized by the widespread development of tight sand containing thin shale beds (Fig. 6), which interferes with our prediction of fracture zones (Wu et al., 2022). These compacted shales exhibit a microstructure with a preferential orientation in the plane (Bandyopadhyay, 2009), which can be confused with fracture-induced anisotropy.

Figure 6 The tight sand containing thin shale beds. The light-colored sections of the core are tight sand, while the dark-striped sections are thin shale beds.

3.2 Model calibration

To verify and calibrate the background parameters of the model, we analyze the sample points within the well. We extracted P-wave and S-wave velocity under three different control variable environments, as shown in Fig. 7. Figure 7a and b show the effect of shale content on the velocity of rock for saturation between 0.6 and 0.8 and porosity between 0 and 0.02. Figure 7c and d show the effect of saturation on the velocity of the rock when the shale content is between 0.09 and 0.11 and the porosity is between 0 and 0.02. Figure 7e and f show the influence of porosity on the velocity of rock when the saturation is between 0.6 and 0.8 and the shale content is between 0.09 and 0.11.

Figure 7 The model results compared with actual well logging measurements. Panels (a) and (b) show the effect of shale content on velocity. Panels (c) and (d) show the effect of saturation on velocity. Panels (e) and (f) show the influence of porosity on velocity.

As previously mentioned, the physical properties in this area are complex, so sample points exhibit some divergence with the hybrid model. However, the variation trends in physical properties are consistently with the model under different conditions. We can utilize these trends to calibrate the model's background parameters (Table 2).

Table 2 Table of background parameters for rock physics modeling. These data primarily originate from the Appendix of The Rock Physics Handbook (Mavko et al., 2009) and represent the most commonly used fundamental rock parameters in the field.

4 Analysis and applications

In this chapter, we focused on the impact of the thin shale beds and fractures on the elastic anisotropy. Based on the petrophysical characteristics of tight sand in the area, three theoretical models were established (Table 3). The fracture parameters were set using the parameters defined in Sect. 3.1, while the background parameters were derived from Table 2. These models were used to verify the elastic anisotropy of different types of tight sand. Model 1 represents the tight sand containing fractures. Model 2 represents tight sand containing fractures and thin shale beds. Model 3 refers to the tight sand containing thin shale beds.

Table 3 The parameters of the three models. The model parameters listed in the table are established based on statistical analysis of well log and geological data. Model 1 corresponds to tight sand containing fracture, Model 2 represents hybrid sand, and Model 3 pertains to sand containing thin shale beds.

The main physical property of a rock is porosity, which indicates whether the rock has enough space to collect and migrate fluids. Therefore, we analyzed the effect of porosity on the acoustic velocity and Thomsen parameters of the three models, as shown in Fig. 8. Figure 8a, b, and c depict the effect of porosity on rock acoustic velocity. The elastic parameters of fractures and thin shale beds decrease with increasing porosity, with Model 1 (fractures) being more sensitive to velocity changes than the other models. On the other hand, Fig. 8d, e, and f depict the anisotropic characteristics of the rock. The anisotropy parameter increases with increasing porosity, especially in Model 1. From the analysis of porosity, we can see that thin shale beds and fractures have a similar trend in their influence on acoustic velocity and anisotropy, but the sensitivity of these elastic characteristics is greater in fractures than in thin shale layers.

Figure 8 The impact of porosity on the three models. Panels (a) to (c) show the impact of porosity on rock acoustic velocity. Panels (d) to (f) show the anisotropic characteristics of the rock influenced by porosity. In the plots, the solid black line represents Model 1 (fractures model), the dashed line represents Model 2 (hybrid model), and the dotted line represents Model 3 (shale model).

Figure 9 The anisotropic results of sample points in the working area. All logging samples were collected from the target and adjacent sections, with classifications based on logging data. In the figure, the black triangles represent tight sand containing fractures, the dark-gray circles represent tight sand without anything, and the light-gray circles represent tight sand containing thin shale beds. The dashed box highlights samples where the two types of tight sand have similar acoustic impedance.

Therefore, we further analyze the Thomsen anisotropy parameter variation with acoustic impedance of the three corresponding theoretical models existing in the well logs. Tight sand containing fractures have different Thomsen anisotropy from tight sand without fractures. On the other hand, tight sand containing thin shale beds and fractures has the same acoustic impedance within the dashed box but varies in anisotropy parameters (Fig. 9). This means that the best way to distinguish them is pre-stack inversion. Applying most methods based on post-stack seismic data remains challenging. To facilitate the subsequent description, we define the portion of tight sand containing fractures that have the same response as tight sand containing thin shale beds as “complex sand”.

To further investigate the pre-stack seismic response characteristics of the sand containing them, we utilize the anisotropic Ruger approximation formula proposed by Wang (2024) to analyze the amplitude variation with the incident angle for both. We apply a three-layer model with different elastic and anisotropy parameters listed in Table 4.

Table 4 The parameters in the table are derived from calibrated rock physics models. The first two rows represent the parameters for a background layer, while the remaining rows represent the calculated elastic parameters under various degrees of fracture development in the target layer.

Figure 10 The pre-stack seismic angle response characteristics in the dashed box. Panel (a) shows the relationship between reflection coefficient and incidence angle. The solid line represents the results for tight sand containing fractures (from the dashed box in Fig. 9), and the dashed line represents the results for tight sand containing thin shale layers. At smaller incidence angles, the reflection coefficients of both tight sands are similar. Panel (b) shows the synthetic seismic records obtained from the convolution of the reflection coefficients.

The pre-stack seismic angle gather and corresponding reflectivity coefficient are shown in Fig. 10. The forward analysis shows that, when the incident angles are small, the reflection coefficients are close and the waveforms are similar. However, as the incident angle increases, the complex sand reflectivity decays faster and the waveform is weakened more significantly compared with tight sand containing thin shale beds.

Figure 11 The response corresponding to different porosities in tight sand (right) compared with the response of tight sand with shale beds (left).

In this paper, considering that post-seismic data are easily obtained and with small data size, we propose the fractures and thin shale beds distinguish method based on the post-seismic data and previous hybrid rock physics model, and we analyze the post-stack seismic response of tight sand based on their physical properties. Excitingly, complex sand corresponds to the low-porosity portion of tight sand containing fractures. This means that tight sand with high porosity is not considered complex sand and significantly differs from both complex sand and tight sand containing thin shale beds (Fig. 11).

Figure 12 The fracture and shale ratio of waveforms in tight sand varying with porosities. The attribute value of 1 corresponds to the complex sand within the dashed box in Fig. 8 with lower attribute values indicating higher porosity (more developed fractures).

To further quantify the waveform differences, we extracted the ratio r₁ of the maximum peak amplitude and peak travel time between 60 and 90 ms as waveform shape attributes of the tight sand containing fractures. Similarly, we extracted the ratio r₂ of the maximum peak amplitude and peak travel time between 35 and 50 ms as waveform shape attributes of the tight sand containing thin shale beds. We define the ratio parameter $\frac{r_{\mathrm{1}}}{r_{\mathrm{2}}}$ (namely fracture and shale ratio) to describe the similarity between the tight sand containing thin shale beds and containing fractures. As shown in Fig. 12, when the porosity is between 0.025 and 0.03, the ratio value approaches 1, which means that the response of the tight sand containing thin shale beds and containing fractures is similar, indicating complex sand. As mentioned earlier, this difference-based new attribute (ratio parameter) can effectively identify the tight sand containing fractures. The quantified results show that ratio values are less than 0.8. Thus, the differences in waveforms between the two tight sands can be used to describe the potential range of the tight sand containing fractures. We calculate the attribute of anisotropic aspect ratio of the target layer. We first select a layer of tight sand containing shale beds as the reference layer to facilitate the evaluation of differences between the target layer and the layer of tight sand containing shale beds. The dashed line in Fig. 13 represents the layer of tight sand containing shale beds (reference layer), and the solid black line represents the target layer. For analysis, we extracted anisotropic aspect ratio attributes from the two layers.

Figure 13 Profile extracted from two tight sand layers. The solid line represents the target layer where fracture evaluation is required, and the dashed line represents the reference layer containing thin shale beds, identified through logging and geological analysis.

The results along the target layer are shown in Fig. 13. As previously discussed, the larger the difference, the smaller the ratio and the more developed the fractures are. Conversely, ratios close to 1 or greater than 1 indicate dense rocks or shale-containing dense sandstones. Therefore, in Fig. 14, black areas represent well-developed fractured dense sandstones, while white areas represent tight sand containing low-porosity or thin shale beds. The new seismic attribute shows a clearer correlation with fault distribution, which is depicted with red lines in Fig. 13.

Figure 14 Fracture and shale ratio map. In the figure, dark areas represent regions with high fracture content, while light and gray areas indicate tight sand or tight sand containing thin shale beds. Red lines mark fault lines.

5 Discussion

Simultaneously studying fractures and thin shale beds in tight sand presents significant challenges. Previous research primarily focused on either fractures or shale but not both together. As discussed in Sect. 3, current techniques for identifying fractures in tight sand using seismic data are limited, particularly in formations with thin shale beds. In Sect. 2, we emphasized methods that couple the elastic characteristics of thin shale beds with tight sand skeletons. Complex models, which include various experimental and derived approaches, often require more computational resources and input parameters than simpler methods like the Gassmann model, especially in processes such as forward and inverse modeling based on rock physics. Despite these complexities, complex models offer advantages in capturing more nuanced aspects of rock behavior and properties. To address these computational demands, we propose a novel approach. During our theoretical model analysis, we summarized a seismic attribute based on the differences between fractures and shale beds. In our subsequent discussions, we further evaluate both the reliability and the limitations of our proposed methodology compared to existing technologies.

5.1 Comparison with other shale models

Current shale models utilized in logging and seismic analysis predominantly rely on SCA and DEM models, yet these fail to consider the preferred orientation of shale plates. This limitation arises due to the inherent complexity in mathematically expressing shale plate orientation and the formidable challenge of measuring the corresponding model parameters. In tight sand containing thin shale beds, the orientation of shale plates significantly influences the elastic properties of the formation.

Figure 15 Stiffness matrix error analysis of the critical porosity. Panels (a) to (d) show the errors between the stiffness coefficients of the core and the stiffness coefficients of the different compaction states (critical porosity). Panel (e) shows the error between the stiffness matrix of the core sample and the stiffness matrices of the different compaction states (critical porosity).

To address this limitation, our approach employs a critical-porosity-based method inspired by Bachrach (2011), which effectively approximates shale plate orientation through shale compaction states. This method circumvents the complexities associated with preferred orientation models while calculating the parameters of shale plate orientation. We further validate our models by applying this approach to Bazhenov shale, with detailed shale core data presented in Table 5.

Table 5 The mineral elastic parameters of Bazhenov shale (Vernik and Liu, 1997). The core matrix primarily consists of quartz, feldspar, and clay minerals, while the pores are predominantly filled with organic matter such as kerogen. Note that n/a represents not applicable.

As shown in Fig. 15, the critical-porosity-based model has the smallest total error at the critical porosity of shale (0.4). This demonstrates that describing the optimal orientation of shale plates based on the shale compaction state is reliable. We further compare the model with other models. Our results and those of Qian et al. (2014) are shown in Table 6. The optimized shale model has a lower root-mean-square error than the other widely used shale physical models at the seismic and well logging scales.

Table 6 The errors of the models. The first row presents the stiffness coefficients calculated from the core measurement data. The last row contains results derived from our rock physics model. The intermediate rows reflect calculations from other scholars' models. The final column indicates the root-mean-square error between the calculated stiffness matrix and the measured values.

5.2 Comparison with Hudson models

The widely used model for tight sand is the Hudson model, which effectively describes the elastic characteristics of thin, coin-shaped fractures. Our method further couples the elastic characteristics of thin shale beds with the Hudson model. We used wells 5 and 202 from the work area to compare the effects before and after coupling.

Figure 16 The comparison of results from two models and logs. The left side of the figure shows the logging curve (black curves) and the Hudson model and new model results (red curves). The blue box highlights formations where thin shale beds are developed. The subplot on the right displays the fitting details of both models within the corresponding formations.

The results are shown on the left of Fig. 16a and b. The three blue boxes refer to the layers that contain both thin shale beds and fractures in the tight sand. Detailed comparison results for these three layers are shown on the right of Fig. 16a and b. The Hudson model struggles to accurately capture the velocity characteristics of both shale and fractures simultaneously. In contrast, our optimized model fully expresses both the low velocity of the thin shale beds and the fractures in the tight sand. This enables a more accurate representation of layers that simultaneously develop thin shale beds and fractures compared to the Hudson model.

5.3 Comparison with other seismic attributes

Before discussing our final seismic attribute results, we first extract some typical seismic attributes that describe tight sand fractures for analysis (Fig. 17). The main idea behind these seismic attributes is that fractures enhance seismic wave reflections, resulting in higher amplitude values. Therefore, attributes such as peak amplitude, dominant frequency, reflection energy, composite absolute amplitude, reflection strength energy in decibels (dB), and reflection strength slope can indicate the presence of fractures when high. Due to stress concentration near faults, there is considerable consistency between the distribution of faults and fracture development zones. As faults extend from north to south and from west to east within the work area, the fault density and consequently the fracture density increase. It is evident that the seismic response of thin shale beds, which have similar elastic characteristics to fractures, interferes with the identification of fractures using conventional seismic attributes, and conventional seismic attributes cannot effectively describe this difference.

Figure 17 Commonly used highlighted seismic attributes. In the figure, dark areas represent regions with high fracture content, while light and gray areas indicate tight sand or tight sand containing thin shale beds. Red lines mark fault lines.

When these seismic attributes for the target layer fail to identify fractures, we need to use other methods, such as rock physics model forward and inversion technology. For forward technology, we need to correct our model using a forward model by comparing it with the actual seismic data to obtain reliable fracture parameters, according to the general forward process based on rock physics models. Alternatively, the elastic parameters obtained through seismic inversion techniques may serve as fracture sensitivity parameters in rock physics. Regardless of the method used, the entire technical system requires microstructural parameters such as pore aspect ratio and mudstone plate aspect ratio, necessitating substantial data support. Additionally, unless there is a simple calculation formula like the Gassmann model, the computation is extensive, with most anisotropic models falling into the latter category. This is a key difficulty in current fracture prediction work.

The new technical process based on the rock physics model proposed in Sect. 4, which utilizes the response characteristics of rock physics analysis to study the target layer through a reference layer, effectively avoids this problem. By analyzing the differences between the two layers, we can effectively remove the interference caused by tight sand containing thin shale layers within the target layer.

Moreover, due to stress concentration near faults, there is significant consistency between fault distribution and fracture development zones. The results shown in Fig. 14 align more closely with geological laws compared to other conventional seismic attributes.

5.4 Limitations and future work

We discuss the method limitations and the potential optimizations from two perspectives. On one hand, there are the optimizations and limitations of the rock physics model for tight sand. On the other hand, there are optimizations and limitations of the seismic attribute extraction techniques based on the rock physics model.

The fracture characterization in this paper is simplified. To fully capture the anisotropy of fractures, a more detailed statistical analysis of fracture dip angles and orientations is required. Fractures with significant differences in dip and orientation should be treated as multiple fracture sets, which can then be superimposed. In these studies, fracture compliance (rather than stiffness) should be the primary focus. We recommend using the LSD model (Michael et al., 1995). Additionally, if fluid flow within fractures is to be studied, fractures and pores across multiple scales must be considered. We suggest adopting the Chapman model (Chapman, 2003). Our regional data show that the area mainly contains horizontal fractures with a common preferred orientation, a sufficiently small aspect ratio, and low fracture porosity and density (Han et al., 2022; Zhang et al., 2022). For our study, which aims to eliminate the VTI anisotropy interference from shale in horizontal fractures, simplifying fractures in the reservoir as a single set of horizontal fractures using the Hudson model is a suitable approach. And the Hudson model, compared to other models, is computationally simpler and requires fewer fracture parameters.

The Brown and Korringa (BK) formula is derived from the Gassmann model. However, this classical fluid substitution model still has limitations. Thomsen (2023) pointed out that the Gassmann model incorrectly applied an open-system theorem to a closed-system environment, violating the assumptions of undrained compressibility. This results in inconsistencies between Gassmann's results and Biot's original theory. The BK derivation partially repeats Gassmann's logical error, conflating open and closed systems, particularly when applying unjacketed compression test results to undrained conditions. To optimize fluid substitution models, we suggest two approaches: firstly, conducting more low-frequency measurement experiments; secondly, integrating digital rock techniques to further derive and summarize the elastic changes caused by fluid flow in complex seepage channels.

On the other hand, the technical process of seismic attribute analysis based on the rock physics model is an attempt to follow the theoretical response analysis of the rock physics model. The results show that this approach is better than conventional seismic attributes. However, this is not applicable to all work areas. This requires that the selected reference layer has stable physical properties. If the work area is large, this may not be achievable. Nevertheless, the corresponding research ideas can be further expanded. We can use the seismic response of tight sand forward modeling, incorporating specific factors, such as thin shale layers and cracks, as a waveforms dictionary. Then, we can use increasingly mature artificial intelligence technology to match the actual seismic response to more accurately and quantitatively explain the details behind the waveform.

6 Conclusions

This study focuses on the previously neglected thin shale beds within tight sand below the log observation scale. Our rock physics model analysis demonstrates that both fractures and thin shale beds within tight sand influence dynamic elastic parameters similarly. We found that porosity significantly affects the elastic and anisotropy parameters of fractures more than it does for thin shale beds. However, subsequent forward modeling further demonstrated that certain tight sands containing fractures exhibit similar responses to those containing thin shale beds. Therefore, relying solely on fracture models can result in errors for tight sands containing thin shale beds. To address this issue, our new model integrates two types of anisotropy: fracture anisotropy and shale anisotropy models. The new model achieves better results in tight sands with thin shale beds compared to the Hudson model.

In our rock physical analysis, we found that pre-stack seismic data could distinguish between fractures and thin shale beds. However, in practical workflows, post-stack seismic data are more widely used due to their convenience. To identify fractures using post-stack data, we further analyzed the physical properties of this complex sand, significantly simplifying the workflow. Fortunately, complex sand has insufficient porosity. We also found that, the greater the porosity, the more distinct the difference between fractures and thin shale beds. As a result, we regarded the aspect ratio of the waveform as a new seismic attribute and applied it to the region. The results show that the new attribute aligns closely with fault distribution and can effectively characterize fracture distribution in tight sand.

Our application and analysis demonstrate the effectiveness of the developed hybrid rock physics model in identifying thin shale beds and fractures, with advantages over conventional methods. However, our model has limitations, including the dip of fractures, the presence of organic matter within shale, and other microscopic factors that may affect its accuracy in describing the elastic parameters of specific fractures and shales. Future research should conduct microscopic experiments on this type of tight sand, optimize the model, enhance its adaptability to tight sands with different microstructures, and explore more practical application scenarios to enhance the potential and practical significance of the research results.