Modeling Method to Characterize the Pore Structure of Fractured Tight Reservoirs

Abstract

The study of unconventional reservoirs has gained increasing attention with the deepening of exploration and development especially in deep-buried tight sandstone reservoirs. We could not obtain the accurate elastic parameters of reservoirs using the conventional rock physics model, since tight sandstone reservoirs have the characteristics of strong heterogeneity, complex lithology and storage space. In order to better describe tight sandstone reservoirs, we improved the traditional tight sandstone rock physics model by combining the dual-connected pore model and the linear slip model. Since the combined modeling process subtly considers four characteristics including the diversity of tight sandstone matrix minerals, the irregularities of pores structure, the connectivity between pores, and the anisotropy caused by fractures, multiple reservoir characteristic parameters can be derived from the limited logging information by the improved model. These reservoir characteristic parameters could account for the difference of diagenesis, which are useful to distinguish different pore types and eliminate ineffective reservoirs. The practical application shows that the improved model can extract microscopic reservoir information hidden in the logging data more effectively than the traditional model. It provides a reliable tool for identifying effective reservoirs in tight sandstone.

1. Introduction

While conventional reservoirs are still a very important part of the global energy supply, the progress in horizontal drilling and hydraulic fracturing technologies makes it possible to exploit unconventional reservoirs [1]. As a typical unconventional oil and gas resource, tight sandstone reservoirs have become an important target for petroleum exploration in China [2]. Effective reservoir characterization is the key to their successful exploration. The rock physics model is the bridge between elastic parameters and subsurface rock physical properties, which have been proven to be a powerful tool for studying reservoir characteristics [3,4]. Deep-buried tight sandstone has often been affected by a series of diagenesis events, such as dissolution, mechanical compaction and cementation, forming a complex pore structure with low permeability, low porosity and abundant micro-fractures [5]. Moreover, the pore structure of rocks is usually characterized in terms of the characteristics of their geometry, connectivity and size of stiff pores and micro-fractures. However, due to the micro-scale and macro-scale heterogeneity of tight sandstone, there is no clear relationship between pore structure and physical properties of the rocks, which brings great difficulty to accurately estimating reservoir elastic parameters via rock physics models. Thus, in order to correctly explore the interrelations between different physical properties of rocks, it is essential to construct a practical rock physics models for tight sandstone.

Many classical rock physics models are empirical and typically assume a linear relationship between velocity and porosity [6,7]. Some parameters used in empirical models can be adjusted according to field data, so these models are suitable and can be applied to reservoirs with similar lithology [8]. However, due to the differences in the mineral composition, pore structure, and fluid properties of underground reservoirs in various regions, these empirical models cannot properly describe the complicated elastic behaviors of tight sandstone with heterogeneous pore structures.

Some scholars consider heterogeneous porous rocks as composite materials composed of mineral particles, pores and micro-fractures [9,10,11,12]. Pores and micro-fractures can be filled with different fluids: oil, gas, water, etc. In this way, it is possible to use the effective medium theory to estimate elastic properties according to the composition and microstructure of rocks at different scales [13,14]. Under long wavelength conditions, heterogeneous mineral particles, pores and micro-fractures can be represented by ellipsoidal-shaped inclusions [15]. Some scholars use self-consistent approximation (SCA) and differential effective medium (DEM) models to determine the elastic properties of isotropic saturated porous media iteratively [16,17]. However, the expressions of the SCA and DEM models are implicit and we need to calculate them iteratively, which reduces the computational efficiency of the model. The Kuster-Toksöz (KT) and Mori-Tanaka (MT) models are widely used in the rock physics modeling of conventional sandstone reservoirs due to their explicit expressions and high computational efficiency [18]. Inclusion-based models assume that each pore in a medium is isolated and distributed randomly, meaning that the KT and MT models may be not suitable for the tight sandstone reservoirs with complicated pore structures. Several scholars have concluded that the MT model can give more consistent results with measured samples than other equivalent medium models, through experiments. With the help of the MT model, they added dry connected pores and saturated isolated pores into a rock matrix to form a dual-connected pore model describing tight sandstone reservoirs, but this model ignores the anisotropic effects of fractures in rocks [19].

Underground fractures connecting isolated pores increase effective porosity, which can be treated as the storage space and migration pathway of tight sandstone reservoirs [20]. The presence of fractures in the subsurface affects the fluid flow and make the propagation of seismic waves azimuthally anisotropic in individual layers. Thus, the estimation of fracture anisotropy parameters is very important for tight reservoir description. There are two commonly used fracture rock physics models. One is the penny-shaped model for cracked media, which assumes that there are thin, penny-shaped ellipsoidal cracks or inclusions in elastic solids [21]. The other is the linear slip model of a fracture medium, which regards fractures as infinitely thin and highly compliant planes regardless of their shape, microstructure, and density, and assumes that the displacement across the fracture surface is linearly correlated with the stress acting on the fracture surface under low-frequency conditions [22]. These fracture rock physics models connect fracture micro-parameters (fracture density, weaknesses, aspect ratio) with rock macro-properties (elastic parameters) and are widely used in characterizing fractures by logging and seismic data [23]. However, it is difficult to characterize the influence of the rock mineral matrix and fluid on rock physical properties only by fracture rock physics models. Hence, some scholars established the anisotropic rock physical model of tight reservoirs by combining the penny-shaped model with the Xu–White model. This combined modeling method is more suitable for underground rocks [24,25].

Therefore, according to the characteristics of storage space in tight sandstone reservoirs, considering the effects of mineral matrix, pore structures, fractures and fluid substitution in anisotropic rocks, we referred to the modeling workflow of the Xu–White model, combined the dual-connected pore model with the linear slip model, and proposed a practical rock physics model scheme for fractured tight sandstone reservoirs. The model can accurately predict the elastic modulus and anisotropic parameters of a formation from limited logging information. The application results in the actual study area show that the variation characteristics of the fracture parameters predicted by the model can directly identify the location of fracture development, which can provide the basis for fine reservoir characterization.

2. Materials and Methods

2.1. Experimental Samples

As shown in Figure 1a, the study area was the Xinchang gas field, which is located in the western depression of Sichuan Basin, China. Large-scale tight sandstone reservoirs are widely developed in the Upper Triassic Xujiahe Formation, and play an important role in natural gas exploration in the study area. The tight sandstone of the Xu-2 Member reservoir is formed in the delta facies of fluvial and marine interaction, which belong to a braided river delta sedimentary system. The reservoir sandstone is mainly clastic sandstone, followed by quartz sandstone and a small amount of feldspar sandstone. Due to long-term geological processes, the Xu-2 Member reservoir is characterized by deep burial, low porosity–permeability and complex reservoir space [26].

Due to the complex tectonic-sedimentary evolution process, diagenesis plays a strong role in controlling the quality and heterogeneity of reservoirs, making the Xu-2 Member reservoir belong to the fracture-pore pattern reservoir. It can be seen from Figure 1b that the higher structural curvature corresponds to the intersection of multiple faults and the bending of strata, which are usually the positions of reservoirs with good pore structure. According to the FMI imaging logging identification results in Figure 1c, the faults in the study area are mainly E-W direction, followed by the N-S direction. Meanwhile, the horizontal maximum principal stress direction of Xu-2 Member is in the E-W direction, and the fractures along this direction have a good opening degree. Therefore, the evaluation of pore structure and natural fractures is the key issue to predict the effective reservoir of Xu-2 Member.

Sixteen rock samples were collected from the tight sandstone strata of Xu-2 Member. They were mainly composed of quartz, feldspar and debris. As shown in

Table 1. The physical properties and mineral composition information of rock samples.

Sample	Total Porosity (%)	P-Wave Velocities (km/s)	Quartz (%)	Feldspar (%)	Debris (%)
1	1.79	5.49	57	11	27
2	2.46	5.37	64	15	17
3	3.26	5.47	61	23	12
4	3.63	5.13	62	16	14
5	4.84	5.33	55	27	11
6	4.92	4.82	66	14	13
7	5.34	5.28	55	31	10
8	5.41	5.05	63	16	17
9	5.37	4.69	68	12	13
10	3.87	4.87	65	13	16
11	6.15	4.61	71	11	13
12	6.03	5.07	62	18	14
13	6.53	4.95	61	17	17
14	7.16	5.16	59	27	10
15	9.05	5.12	60	28	8
16	7.81	4.88	55	12	31
Mean values	5.23	5.08	61.52	18.19	15.19
Standard deviation	1.88	0.26	4.54	6.54	5.87

, quartz was in the range from 55% to 71% with an average of 61%, feldspar content varied from 11% to 31% with an average of 18%, and the rock debris was in the range from 8% to 31% with an average of 15%. High compositional maturity and feldspar content are favorable for the development of dissolution pores and micro-fractures in reservoirs.

At the same time, it can be found from Table 1 that the porosity distribution range of the rock samples was 1.79–9.05%, with an average of 5.23%. However, the P-wave velocity (V_P) range of the rock samples was 4.5–5.4 km/s, and the relationship between velocity and porosity was relatively dispersed. In the case of small changes in porosity of rock sample, the velocity difference of each sample was obvious, which was very unfavorable for reservoir evaluation. Therefore, it was necessary to analyze the microscopic pore structure of the rock samples by micro-CT scanning.

Detailed analysis of rocks at the core scale is essential because it helps to establish reliable models at the initial scale [27]. The micro-CT scanning images of rock samples were reconstructed using Avizo software to obtain a three-dimensional reconstruction map, and the three-dimensional reconstruction map of some samples was selected as shown in Figure 2. It can be seen from Figure 2 that the tight sandstone samples have the characteristics of different pore throat sizes and irregular pore shapes at the micron scale. The pores of rock samples mainly included contiguous pores and isolated pores, and the connectivity of the contiguous pores was better than that of the isolated pores. Due to the disordered distribution of the pores, the rock samples had obvious microscopic heterogeneity.

In order to observe the distribution of microscopic pore system, we sliced the three-dimensional digital core of the rock sample. The result of rock sample slicing is shown in Figure 3a. We found from the rock section that there were so many kinds of pores with different size, shape and connectivity, the tight sandstone can be regarded as a kind of medium with multiple pore structures under ellipsoid inclusion theory [28]. Therefore, the pores in the rock section were fitted to ellipses with different pore aspect ratios using an image processing algorithm [29]. It was calculated that the pore aspect ratio distribution range of rock samples was 0.01–0.65, with an average of 0.21 and a standard deviation of 0.16. Figure 3b shows the equivalent slicing of the rock pore structure after digital image processing. According to the equivalent slicing of the rock pore structure, different types of pores in the rock were represented by some ellipsoids with different aspect ratios. Thus, we set the initial pore aspect ratio of the primary intergranular pores, secondary dissolution pores and micro-fractures at 0.25, 0.45, and 0.05, respectively. This provided a reference for the parameter settings of the subsequent rock physical modeling.

Tight sandstone reservoirs tend to exhibit strong elastic anisotropy resulting from fluid-filled aligned fractures. Since the wavelength of a seismic wave is much larger than the fracture aperture, we can ignore the opening degree of the fracture aperture and the details of fracture spatial distribution, and treat the fracture equivalently as an anisotropic solid in the rock physics modeling [30]. The anisotropic characteristics of the rock will depend on the development direction and intensity of the fractures, the elastic properties of the material filling the fracture, and the elastic modulus of the surrounding rock.

The statistical results of 318 structural fractures in the study area are shown in Figure 4. It can be seen from the figure that vertical fractures and sub-vertical fractures were mainly developed, and the number of low angle fractures and horizontal fractures was less in the study area. Moreover, vertical or sub-vertical fractures were mostly effective fractures with good opening degrees and partial fracture width up to 0.1–0.5 cm. In contrast, the low-angle fractures were mostly ineffective and closed due to the gravity of the overlying strata.

Fractured reservoir rocks filled with vertical or sub-vertical arranged fractures may exhibit horizontal transverse isotropic (HTI) characteristics based on the fracture equivalent medium theory [31,32]. According to the statistical results, the Xu-2 Member reservoir can be assumed to be an HTI medium, which is formed by a group of vertically arranged non-interactive fractures embedded in a pure isotropic background medium. This assumption helps us to describe the propagation characteristics of seismic waves in tight sandstone of study area [33].

2.2. Modeling Method

In this section, we refer to the working process of the Xu–White model to establish an improved anisotropic fractured rock physics model [34,35].

2.2.1. Calculating the Matrix Mineral Elastic Modulus

The equivalent elastic modulus of a rock matrix is determined by the elastic modulus and the volume content of rock-forming minerals, of which the volume content of minerals can be obtained from core and log data [36]. Assuming that the rock matrix is a mixture of quartz, feldspar and debris, in which the debris is mainly composed of carbonate minerals, once the percentage of each mineral is determined, the Voigt–Reuss–Hill (VRH) mixing model is further used to calculate the elastic modulus of rock matrix, as follows [37]:

$$M_{VRH}=\frac{M_V+M_R}{2}$$

(1)

where, $M_v={\displaystyle \sum _{i=1}^n{f}_i}{M}_i$ and $\frac{1}{M_R}={\displaystyle \sum _{i=1}^n\frac{f_i}{M_i}}$ are the upper and lower limits of the matrix mineral elastic modulus calculated by the Voigt and Reuss models, respectively. $f_i$ and $M_i$ are the volume fraction and modulus of mineral component, and $M_{VRH}$ is the equivalent elastic modulus of the matrix minerals.

2.2.2. Calculating the Rock Matrix Elastic Modulus

According to the observation of the thin sections, it was found that the pore types of tight sandstone reservoirs mainly include intergranular pores, dissolution pores and micro-fractures. Thus, the total porosity ${\varphi }_t$ of rock can be expressed by the following equation:

$${\varphi }_p+{\varphi }_s+{\varphi }_f={\varphi }_t$$

(2)

where ${\varphi }_p,{\varphi }_s,{\varphi }_f$ denote intergranular pores, dissolution pores and micro-fractures, respectively. Each type of porosity can be calculated using logging curves [38], and the detailed calculation formula is shown in Appendix A.

Assuming that dissolution pores and micro-fractures form connected pores of rocks, and intergranular pores form isolated pores of rocks, the connectivity between pores can be calculated by the following equation:

$$v_i=\frac{{\varphi }_i}{{\varphi }_t}$$

(3)

$$\xi =\frac{{\varphi }_{con}}{{\varphi }_t}$$

(4)

The scale factor $v$ is the ratio defined by the Equation (3), where ${\varphi }_i$ is the porosity of the ith type of pore. According to the connectivity between pores, they can be divided into connected pores ${\varphi }_{con}$ and isolated pores ${\varphi }_{iso}$, and the connectivity of pores can be described by the connectivity coefficient $\xi$ based on the Equation (4).

Although dissolution and tectonic processes increase the connectivity between pores, there were still a certain proportion of disconnected pores in the rock. Since fluid in disconnected pores has difficulty flowing and substituting, when elastic wave propagates in this kind of porous medium, fluids in pores with different shapes and connectivity will produce different responses. Therefore, in order to satisfy the assumption that all pores are connected in the classical fluid substitution theory, the disconnected pores are regarded as a part of the rock–solid matrix, and the fluid substitution only occurs in the connected pores.

Due to the inclusion-based effective medium theory and fluid substitution equation being effectively able to characterize the influence of different pore shapes on the frame modulus, we assumed that the rock matrix is composed of matrix minerals and isolated pores, and adopted the MT model by considering the pore shape and the interaction of adjacent inclusions to couple the matrix minerals with the saturated isolated pores. Thus, the rock matrix modulus can be calculated by the following equation:

$${\varphi }_{mat}=\frac{{\varphi }_t-{\varphi }_{con}}{1-{\varphi }_{con}}$$

(5)

$$\{\begin{array}{l}{K}_o=\frac{(1-{\varphi }_{mat})K_{ma}+{\varphi }_{mat}{\displaystyle \sum _{i=1}^n{\upsilon }_i{K}_w{P}_i}}{(1-{\varphi }_{mat})+{\varphi }_{mat}{\displaystyle \sum _{i=1}^n{\upsilon }_i{P}_i}}\\ {\mu }_o=\frac{(1-{\varphi }_{mat}){\mu }_{ma}}{(1-{\varphi }_{mat})+{\varphi }_{mat}{\displaystyle \sum _{i=1}^n{\upsilon }_i{Q}_i}}\end{array}$$

(6)

where the matrix porosity ${\varphi }_{mat}$ can be calculated by Equation (5). $P_i$ and $Q_i$ are the geometric factors of the ith type of saturated isolated pores added to the rock matrix, which is the function of the fluid modulus and the pore aspect ratio [39]. $K_w$ is the bulk modulus of the bound water in the matrix pores, $K_o$ and ${\mu }_o$ are the rock matrix elastic modulus, and $K_{ma}$ and ${\mu }_{ma}$ are the mixed minerals’ elastic modulus.

2.2.3. Calculating the Elastic Modulus of Porous Dry Rock

The elastic modulus of dry rock skeleton is calculated using the MT model by adding dry connected pores into the rock matrix; the calculation equation is as follows:

$$\{\begin{array}{l}{K}_{dry}=\frac{\left(1-{\varphi }_t\right)K_{ma}+\left({\varphi }_t-{\varphi }_{con}\right){\displaystyle \sum _{i=1}^n{\upsilon }_i{K}_w{P}_i}}{\left(1-{\varphi }_t\right)+\left({\varphi }_t-{\varphi }_{con}\right){\displaystyle \sum _{i=1}^n{\upsilon }_i{P}_i+{\varphi }_{con}{\displaystyle \sum _{i=1}^n{\upsilon }_i{\widehat{P}}_i}}}\\ {\mu }_{dry}=\frac{\left(1-{\varphi }_{mat}\right){\mu }_{ma}}{\left(1-{\varphi }_t\right)+\left({\varphi }_t-{\varphi }_{con}\right){\displaystyle \sum _{i=1}^n{\upsilon }_i{Q}_i+{\varphi }_{con}{\displaystyle \sum _{i=1}^n{\upsilon }_i{\widehat{Q}}_i}}}\end{array}$$

(7)

where ${\widehat{P}}_i$ and ${\widehat{Q}}_i$ are geometric factors of the ith type of empty connected pores added to the rock skeleton, $K_{dry}$ and ${\mu }_{dry}$ are the elastic modulus of dry rock skeleton.

2.2.4. Calculating the Stiffness Matrix of Fractured Dry Rock

Assuming that the rock matrix is an isotropic medium, the stiffness matrix of dry rock skeleton $C_0$ is as follows [40]:

$$C_0=\left[\begin{array}{cccccc}{\lambda }_{dry}+2{\mu }_{dry}& {\lambda }_{dry}& {\lambda }_{dry}& 0& 0& 0\\ {\lambda }_{dry}& {\lambda }_{dry}+2{\mu }_{dry}& {\lambda }_{dry}& 0& 0& 0\\ {\lambda }_{dry}& {\lambda }_{dry}& {\lambda }_{dry}+2{\mu }_{dry}& 0& 0& 0\\ 0& 0& 0& {\mu }_{dry}& 0& 0\\ 0& 0& 0& 0& {\mu }_{dry}& 0\\ 0& 0& 0& 0& 0& {\mu }_{dry}\end{array}\right]$$

(8)

where $C_{ij}$ is the stiffness element of dry rock skeleton, and the Lamé parameter ${\lambda }_{dry}=K_{dry}-\frac{2}{3}{\mu }_{dry}$, where $K_{dry}$ and ${\mu }_{dry}$ can be calculated by Equation (7).

As discussed in the previous section, the Xu-2 Member reservoir can be assumed to be an HTI medium. An HTI medium can be considered as a simple combination of fractures and isotropic background rocks. The isotropic background rocks are characterized by the Lamé parameter $\lambda$ and the shear modulus $\mu$, while the imbedded fractures are characterized by the dimensionless normal and tangential fracture weakness parameters ${\delta }_N$ and ${\delta }_T$, respectively [41]. By using these four parameters, expressions for the stiffness matrix of fractured dry rock $C^{dry}$ can be written as

$$C^{dry}=\left[\begin{array}{cccccc}{M}_{dry}\left(1-{\delta }_N\right)& {\lambda }_{dry}\left(1-{\delta }_N\right)& {\lambda }_{dry}\left(1-{\delta }_N\right)& 0& 0& 0\\ {\lambda }_{dry}\left(1-{\delta }_N\right)& M_{dry}\left(1-{\left(\frac{{\lambda }_{dry}}{M_{dry}}\right)}^2{\delta }_N\right)& {\lambda }_{dry}\left(1-\frac{{\lambda }_{dry}}{M_{dry}}{\delta }_N\right)& 0& 0& 0\\ {\lambda }_{dry}\left(1-{\delta }_N\right)& {\lambda }_{dry}\left(1-\frac{{\lambda }_{dry}}{M_{dry}}{\delta }_N\right)& M_{dry}\left(1-{\left(\frac{{\lambda }_{dry}}{M_{dry}}\right)}^2{\delta }_N\right)& 0& 0& 0\\ 0& 0& 0& {\mu }_{dry}& 0& 0\\ 0& 0& 0& 0& {\mu }_{dry}\left(1-{\delta }_T\right)& 0\\ 0& 0& 0& 0& 0& {\mu }_{dry}\left(1-{\delta }_T\right)\end{array}\right]$$

(9)

where the normal and tangential weaknesses parameters ${\delta }_N$ and ${\delta }_T$ can be calculated from the normal and tangential fracture compliances of the fractured rocks [42]. The P-wave modulus $M_{dry}=K_{dry}+\frac{4}{3}{\mu }_{dry}$, can also be calculated by Equation (7).

2.2.5. Calculating the Stiffness Matrix of Fractured Saturated Rock

Assuming that the dry isotropic porous rocks with interconnected vertical fractures contribute little to the porosity of rocks [43], according to the stiffness matrix of fractured dry-rock, we can use the anisotropic fluid substitution equation to substitute a mixture of hydrocarbon and water into the pore system of rock matrix. The anisotropic fluid substitution equation is as follows [44]:

$$C_{ij}^{sat}=C_{ij}^{dry}+a_i{a}_j\frac{K_o}{\left(1-\frac{K^{*}}{K_o}\right)-{\varphi }_t\left(1-\frac{K_o}{K_f}\right)},\hspace{1em}i,j=1,2,\cdots ,6$$

(10)

$$K^{*}=\frac{1}{9}{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^3{C}_{ij}^{dry}}}$$

(11)

$$\frac{1}{K_f}={\displaystyle \sum _{i=1}^N\frac{f_i}{K_i}}$$

(12)

$$a_m=b_m-\frac{{\displaystyle \sum _{n=1}^3{C}_{mn}^{dry}}}{3K_o},m=1,2,\cdots ,6$$

(13)

$$b_m=\{\begin{array}{ll}1,\hfill & m=1,2,3,\hfill \\ 0,\hfill & m=4,5,6,\hfill \end{array}$$

(14)

where $C_{ij}^{sat}$ is the stiffness element of fractured saturated rock, $K_f$ is the bulk modulus of the mixed fluid, $f_i$ is the volume fraction of each fluid, and $K_i$ is the bulk modulus of each fluid. The specific derivation and calculation results of the stiffness element of fractured saturated rock are shown in Appendix B.

2.2.6. Calculating the Anisotropic Parameters of Fractured Saturated Rock

After obtaining all the stiffness elements of the saturated rock, we can calculate the Thomsen-style anisotropic parameters using the following equations [45,46]:

$${\epsilon }^{(v)}=\frac{C_{11}^{sat}-C_{33}^{sat}}{2C_{33}^{sat}}$$

(15)

$${\gamma }^{(v)}=\frac{C_{66}^{sat}-C_{44}^{sat}}{2C_{44}^{sat}}$$

(16)

$${\delta }^{(v)}=\frac{(C_{13}^{sat}+C_{55}^{sat}{)}^2-{(C_{33}^{sat}-C_{55}^{sat})}^2}{2C_{33}^{sat}(C_{33}^{sat}-C_{55}^{sat})}$$

(17)

where ${\epsilon }^{(v)}$ can be used to describe the anisotropy of P-wave velocity in HTI medium, ${\gamma }^{(v)}$ can be used to describe the anisotropy of S-wave velocity in HTI medium, and ${\delta }^{(v)}$ is the change rate of P-wave NMO velocity.

By the relationship between the Thomsen-style anisotropic parameters and the fracture weaknesses, Equations (15)–(17) can be linearized, under the assumption of weak anisotropy [47]:

$${\epsilon }^{(v)}=-2g(1-g){\mathsf{\Delta }}_N$$

(18)

$${\gamma }^{(v)}=-\frac{{\mathsf{\Delta }}_T}{2}$$

(19)

$${\delta }^{(v)}=-2g[(1-2g){\mathsf{\Delta }}_N+{\mathsf{\Delta }}_T]$$

(20)

where ${\mathsf{\Delta }}_N$ and ${\mathsf{\Delta }}_T$ are the normal and tangential weaknesses parameters of fractured saturated rock. The parameter $g=V_s^2/V_p^2$ is the S-wave to P-wave velocity ratio to the square, where $V_p$ and $V_s$ are the P- and S-wave velocities in the background medium.

Since the vertical P-wave and S-wave velocities are parallel to the fracture direction in HTI media, the formula for calculating the P-wave and S-wave velocities of background rock are as follows, under the assumption of weak anisotropy [48]:

$$V_p=\sqrt{\frac{C_{33}^{sat}}{\rho }}$$

(21)

$$V_s=\sqrt{\frac{C_{44}^{sat}}{\rho }}$$

(22)

Additionally, we can get the fracture density $e$ from the tangential weaknesses using the following equation [49]:

$$e=\frac{3\left(3-2g\right)}{16}{\mathsf{\Delta }}_T$$

(23)

Therefore, all the anisotropic parameters and elastic modulus of formation can be derived from the well logging data.

2.3. The Workflow of Rock Physics Model Construction

According to the analysis above, we obtained a workflow for the rock physics modeling of tight sandstone. As shown in Figure 5, it comprises four steps:

Step 1. Calculate the elastic modulus of the rock matrix minerals using the VRH equation.

Step 2. Append the connected dry pores and the isolated pores containing bound water into the rock matrix based on the MT model, and then calculate the elastic modulus of the porous dry rock.

Step 3. Add the fractures into the dry-rock skeleton using the linear slip fracture model, and then calculate the effect of fracture on the dry rock elastic modulus.

Step 4. Mix fluid into the pore system of rock utilizing the anisotropic fluid substitution equation, and finally obtain the elastic modulus of the saturated fractured rock.

3. Results

3.1. Analyzing Pore Structure Characteristics Using the Rock Physics Model

The relationship between the velocity and porosity of the rock samples in the study area becomes very complicated because the tight sandstone reservoirs have undergone a series of diagenetic transformations. Figure 6 shows the cross plot of P-wave velocity (V_P) versus porosity in the dry rock samples.

It can be seen from Figure 6 that although V_p generally decreases with increasing porosity from the overall trend, the data points present a certain degree of scatting distributions, so the porosity is not conclusive to subsequent seismic quantitative interpretation.

The established rock physics model was used to simulate the variation trend of the P-wave velocity, it shows that the velocity depends on the aspect ratio (AR) of the rock pores, itself depending on the feldspar content. As the dashed lines in Figure 6 show, we found that the scatter was mainly ascribed to the variation of pore structure. Samples with high feldspar content are close to the trends of spherical or tubular pores, while primary matrix pore samples were near the crack-like trend. At a given porosity, the greater the number of secondary dissolution pores, the higher the velocity of the sample.

Additionally, we used the established rock physics model to explore the relationship between the pore structure parameters and the rock elastic properties. As shown in Figure 7, the velocity difference of three rock samples A1–A3 with the same lithology and similar porosity reached 800 m/s. By observing the thin sections of the rock samples, it was found that sample A1 only contained dissolved pores, sample A2 had intergranular pores and dissolved pores together, and sample A3 developed micro-fractures.

Due to the development of micro-fractures and dissolved pores increasing the connectivity between pore systems, the interaction between the fluid and the rock matrix is also enhanced, which greatly reduces the velocity. The pore aspect ratio divides the velocity–porosity relationship into three regions with different characteristics of pore structure which can be used to characterize reservoirs with different pore types.

3.2. Analyzing Reservoir Elastic Properties Using the Rock Physics Model

The established rock physics model was used to predict the shear modulus of the dry rock samples, and the results were compared with those predicted by the traditional Xu–White model. The comparison results are shown in Figure 8.

Since the traditional Xu–White model regards rock pores as fixed types and does not consider the connectivity between rock pores, the prediction results of rock samples with high feldspar content in a complex pore structure system have a large deviation, and the mean relative error of the measured and predicted shear modulus reached 13.05%. The improved model showed good prediction results for these rock samples with different pore structure systems on account of its reasonable modeling scheme. The mean relative error of the measured and predicted shear modulus was 5.62%, so the prediction accuracy of the improved model was better than the traditional Xu–White model.

3.3. Predicting Reservoir Characteristic Parameters Using the Rock Physics Model

Well A in the research area was selected to further verify the validity of the established rock physics model. Its reservoir characteristic parameters, such as S-wave velocity (Vs), pore aspect ratio and fracture density, were calculated from conventional logging data, including acoustic velocity, density, porosity, water saturation and mineral composition content curves, using the established rock physics model [50].

As shown in Figure 9, although all the rock physics models provide reasonable estimates of Vs, a better result of Vs prediction was produced by the established rock physics model, and the estimated results were in good agreement with the true Vs Data, with the relative error being less than 10%. Compared with the traditional Xu–White model for the prediction of Vs, the improved rock physics model gives a satisfactory result for velocity prediction.

Next, we analyzed the results of estimated pore aspect ratio and fracture density utilizing the imaging log and the reservoir lithofacies data. As shown in Figure 10, the red rectangular regions with high pore aspect ratio and fracture density correspond withthe dissolution-dominated and fracture developing zone, which is the position of favorable reservoirs. The change trend of the pore aspect ratio and fracture density conformed to the logging response characteristics of favorable lithofacies, which can effectively indicate the location of high quality oil–gas reservoirs.

4. Discussion

Porosity is not the only factor affecting the elastic behavior of tight sandstone with low porosity and permeability. Under the same porosity condition, the change of seismic wave elastic characteristics mainly depends on the structural characteristics of the pores, rather than porosity. If the tight sandstone aspect ratio is assumed to be constant, a large prediction error will be introduced in an area with complex geological conditions. Hence, according to the characteristics of micro-pore structures in tight sandstone, we equated its pore space to the combination of multiple types of pores, based on inclusion theory. Meanwhile, in order to further consider the connectivity between rock pores, the connected pores and isolated pores were added to the rock matrix using logging interpretation data and the MT model. Through comparison with the experimental data, the improved model gave more accurate prediction results, which demonstrates its good applicability to the study area.

Obviously, dissolution and fracturing contribute to increasing the storage space and the pore connectivity of tight sandstone reservoirs, and further tends to lead to the formation of high-quality reservoirs. By comparing the test data of the tight sandstone samples with the calculation results of the theoretical model, the pore aspect ratio explained the heterogeneity of the tight sandstone, and the fracture density effectively indicated the location of fracture development. Therefore, this model can provide an effective basis for the logging evaluation and seismic quantitative interpretation of tight sandstone reservoirs.

Although we can obtain accurate microstructure information from micro-CT scanning images of rock samples, the measurement results are often limited by economic costs. Considering the problem of sample representativeness and upscaling in integrating the core and logging data, it is not suggested to use the equivalent aspect ratio of a single rock sample as the initial constraint condition. If there are sufficient core measurement data from the study area, the aspect ratio mean value and variance should be used to represent the statistical normal distribution of pores with varied aspect ratios in the multiple pore aspect ratio model. Doing so will help the improved model achieve better prediction results.

Since the elastic constants of minerals and fluids used in the modeling process are based on the reference values obtained from previous work, the reference values of different layers and regions vary greatly, which affects the prediction accuracy of the model. Furthermore, the acquisition quality of logging data and the calculation accuracy of reservoir parameters, such as porosity and saturation, also affect the final prediction results of the theoretical model.

Furthermore, considering that the rock pore space is composed of stiff pores (insensitive to pressure) and micro- fractures (sensitive to pressure) distributed in the rock [51], because the effective aspect ratios of the soft pores and open cracks decrease under increasing pressure, the soft pores and micro-fractures of rock samples observed in lab conditions are different than in high pressure reservoir conditions. Hence, the velocity of rock samples measured under laboratory conditions will be biased, which affects the evaluation of the prediction effect of the model. One potentially significant over-simplification was that the pressure dependence of the porosity was ignored in our modeling process, and more work is required to develop a proper effective model in the future.

5. Conclusions

The improved anisotropic rock physics modeling method considers the effects of connectivity between different shapes of pores and the anisotropy caused by fractures. Our approach can estimate shear-wave velocity, pore aspect ratio and fracture density more accurately from the conventional well logging data, as compared with the traditional rock physics model. It provides the key parameters to other workflows, such as logging evaluation, seismic anisotropy inversion and fracture development in tight sandstones. In addition, our modeling workflow achieved good results in the application of the fine characterization of the deep-buried tight sandstone reservoirs with strong heterogeneity.

It should be noted that this model can achieve good application results under the condition that the geological situation conforms to the assumption of effective medium theory. However, it may be that not all the tight sandstones are the same as in our study area. Therefore, the applicability of this rock physics model should be carefully verified before applying it to other regions.

The theoretical model cannot fully explain the elastic properties of real tight sandstones underground because the simplified pore and fracture geometry models using effective theory are still far away from the structure of actual rock. Therefore, we still have more work to do to develop more accurate and practical rock physical models in the future.