Estimation of the Transverse Wave Velocity in Siliceous Carbonate Reservoirs of the Dengying Formation in the Gaoshiti–Moxi Area, Sichuan Basin, China

Abstract

Siliceous minerals of the Dengying Formation in the Gaoshiti–Moxi area in the central Sichuan Basin exhibit four types of quartz crystals (cryptocrystalline quartz, chalcedony, microcrystalline quartz, and megacrystalline quartz) and three structural types: cryptocrystalline, microcrystalline, and mosaic (laminated mosaic, window-hole interrupted mosaic, and arc-laminated mosaic). Siliceous minerals have a great influence on the storage performance of the reservoirs in the Dengying Formation. According to the petrophysical parameters of the Dengying Formation and porosity intersection diagrams, the siliceous dolomite and the reservoirs have low impedance characteristics, which makes it difficult to distinguish between them and leads to difficulties in the characterization and prediction of the reservoirs. The transverse wave velocity is favorable for reservoir characterization. Currently, the main method used to estimate the transverse wave velocity is petrophysical modeling, which establishes a relationship between the elastic and physical parameters of the reservoir. In this paper, the siliceous minerals in the dolomite in the study area are regarded as solid inclusions, and the calculation method of the rock matrix modulus is improved by using solid replacement. Then, an improved petrophysical model is constructed by combining the KT (Kuster–Toksöz) model, the DEM (Discrete Element Method) model, the Gassmann equation, and the Wood equation. The transverse wave velocity is estimated using the improved model under the constraint of the longitudinal wave velocity. The shapes of the transverse wave velocity curves obtained by the improved model and the deviations from the measured velocities are significantly better than those of the Xu–Payne model and other models. The results show that the improved model can effectively estimate the transverse wave velocity of the reservoir in this area, which provides a basis for future reservoir predictions in this area.

1. Introduction

Carbonate oil and gas reservoirs occupy an important position in the distribution of oil and gas around the world. The distribution area of carbonate rocks accounts for about 20% of the total area of sedimentary rocks around the world; carbonate oil and gas reserves account for about 50% of the world’s total oil and gas reserves; and carbonate oil and gas production accounts for about 60% of the world’s total oil and gas production. Many important oil and gas areas around the world are dominated by carbonate reservoirs. Oil and gas fields composed of carbonate reservoirs have the characteristics of large reserves and high yields of single wells, which make it easy to form large-scale oil and gas fields. There are a total of nine high-yield wells around the world that produce more than 10,000 tons per day, of which eight are carbonate oil and gas fields. At present, conflicts between international energy supply and demand are prominent, and energy security is increasingly becoming the focus of attention of all countries. The exploration and development of increasingly complex and hidden targets, carbonate rocks, and other special lithological reservoirs have gradually become a hot spot in oil and gas exploration and development.

The area of carbonate rocks in China accounts for about 30% of the country’s land area, of which the Tarim Basin, Sichuan Basin, Ordos Basin, and North China are widely developed and are potential areas for oil and gas exploration. Starting from the breakthrough of the GS1 well in the Dengying Formation in the Gaoshiti–Moxi area of the central Sichuan Basin, several gas wells with capacities of millions of cubic feet have been found, and carbonate reservoirs of the Dengying Formation have become the key areas for deep natural gas exploration and development in the Sichuan Basin [1]. Siliceous minerals of the Dengying Formation in the Gaoshiti–Moxi area of the central Sichuan Basin exhibit four types of quartz crystals (cryptocrystalline quartz, chalcedony, microcrystalline quartz, and megacrystalline quartz) and have three structural types: cryptocrystalline, microcrystalline, and mosaic structures (layered inlays, interrupted window-shaped hole inlays, and curved layered inlay structures) [2]. Different types of siliceous minerals have different structures, which leads to non-homogeneity of the reservoirs in terms of their lithology and storage properties, especially in deep carbonate karst formations which are often difficult to characterize because of their greater burial depth and longer evolution time, resulting in unusually complex reservoir mineral components, pore structures, and cementation types [1]. With the increasing demand for exploration and development, more accurate quantitative interpretation of seismic forward and inverse results and more accurate identification of reservoir lithology and fluids are required [3]. The transverse wave velocity is suitable for describing reservoir characteristics and has a very important role in the process of pre-stack seismic inversion, AVO attribute analysis, etc., which can help realize the above objectives. However, in actual production, many wells lack transverse wave data due to the high cost of transverse wave logging, immature acquisition technology, and difficult data interpretation [3]. Therefore, transverse wave velocity estimation methods have always been a popular but difficult research topic.

Previous studies have shown that petrophysical modeling is a bridge to establish a connection between the elastic and physical parameters of reservoirs, which implies the quantitative characterization of the “physical parameter → elastic parameter” relationship, and it has become the main transverse wave velocity estimation method. However, in the dolomite reservoirs of the Dengying Formation, siliceous minerals are widely present [4] and their genesis is complex, so their mineral structure characteristics are different [5]. Therefore, the influence of siliceous mineral structural features on the petrophysical parameters needs to be considered when establishing a petrophysical model to ensure its accuracy [6].

Currently, the most commonly used carbonate rock model is the extension of the Xu–White model [7] by Xu and Payne [8], which was further improved by Sun et al. [9] and Zhang et al. [10]. These models mainly focus on describing the pore structure, and the pore structure parameters (porosity, pore throat radius, etc.) are usually calculated using iterative methods [11]. However, when using the above models, when the structural characteristics of siliceous minerals in carbonate reservoirs change, the usual model used to obtain the matrix mineral modulus cannot respond, which causes the changes due to the mineral structural characteristics to be transferred to the pore structure characteristic parameters; ultimately, accurate pore characteristic parameters cannot be obtained, and model calculations cannot be performed. Li Jingye et al. [12] divided the rock using a small grid according to the pore features to calculate the modulus of elasticity, but if the mineral structure features are considered in this method, the calculation process becomes more complicated. Cao Xiaochu et al. used a combination of SCA (Self-Consistent Approximation) and DEM (Discrete Element Method) models, but also did not state the meaning of the geometric factors used in constructing the rock matrix [13]. Regarding the effect of diagenesis (dissolution, dolomitization, and recrystallization) on the model, Pride et al. introduced consolidation parameters to establish the relationship between the dry rock modulus and matrix modulus [14]. However, this model only considers the factor of diagenesis and neglects the effect of pore structure on the longitudinal and transverse wave velocity.

In general, the existing physical models of carbonate rocks cannot effectively characterize the complex siliceous minerals contained in the strata of the study area. Therefore, it is difficult to estimate the transverse wave velocity of the siliceous dolomite strata of the Dengying Formation in the Gaoshiti–Moxi region of central Sichuan using the existing models. To address the above problems, this study improves the calculation method for the rock matrix modulus by adding a solid substitution method and estimates the transverse wave velocity of the dolomite reservoirs in this area by using the improved model. By comparing the fitting degree of its curves with the measured transverse wave velocity curves, the validity of the model is verified, laying a foundation for reservoir prediction in the study area in the later stage.

2. Carbonate Reservoir Characteristics

Siliciclastic development is a common occurrence in carbonate reservoirs within the Sichuan Basin. Based on their distinctive output characteristics, these reservoirs can be categorized into four main types, as shown in

Table 1. Relationship between siliceous characteristics and reservoirs in carbonate rocks in the Sichuan Basin.

	Source of SiO2 and Deposition Mechanism	Characteristics	Relationship with the Reservoir
Interlayer silica	Volcanism brings SiO2 into the ocean through various pathways, and the deposition mechanism is dominated by purely chemical interactions	Has a shell-like fracture, a cryptocrystalline–microcrystalline structure, fewer radial aggregates of fibrous chalcedony, little lateral extension, and large variation in the thickness of the single layer	Large sediment thickness and wide distribution, which has a hindering effect on carbonate reservoir development
Mixed silica	Terrestrial quartz grains are transported to shallow marine environments to co-deposit with carbonate rocks to form mixed silica, and the deposition mechanism is mixed	The sorting and rounding are excellent, and the quartz crystals are floating between the carbonate grains, and the two are often filled with clay and other materials, which is slightly disordered	Hard nature, strong resistance to weathering, not easily dissolved, difficult to form storage space in, resulting in the formation of mixed rock lithology; dense
Silicone- filled	The SiO2-rich volcanic debris material and water vapor enter the seawater and chemically precipitate in the dissolution pore crevices to form authigenic quartz	Authigenic quartzes can be in direct contact with each other or with authigenic dolomite, or filled with bitumen simultaneously in the hole	Hard nature, strong resistance to weathering, not easily re-dissolved after filling in the dissolved holes
Interchangeable silica	SiO2-rich seawater chemically changes with and converts minerals in carbonate rocks completely or partially to silica, forming accountable silica	Silica is present in unaccounted-for minerals in the form of raw minerals or components	The original non-dolomite grains will be accounted for as dolomite grains, resulting in intergranular holes, while dolomite is prone to dissolution, resulting in dissolution holes

.

The depositional environment of the Dengying Formation in the study area is dominated by carbonatite terrace phases, with the development of algal dolomite, crystalline dolomite, sand (granular) dolomite, and sand (granular) debris dolomite, with a few thinly bedded sandstones, mudstones, siliciclastics, and paste rocks [15].

Siliceous minerals in the study area have the following characteristics: ① They are dense and hard, often with streaks and bands; radiolarian chalcedony and cryptocrystalline siliceous minerals can be observed under the microscope (Figure 1). ② Carbonatite reservoirs have high SiO₂ content and low TiO₂, Al₂O₃, and MgO contents (

Table 2. The main element content of siliceous carbonate rock samples from the fourth member of the Dengying Formation in the Gaoshiti-Moxi area.

Main Elements	Sample 1	Sample 2	Sample 3	Sample 4
	GS18 Well	GS18 Well	GS20 Well	GS20 Well
	5184 m	5207 m	5229 m	5230 m
SiO2	94.11%	96.65%	96.80%	95.85%
TiO2	0.05%	0.03%	0.04%	0.04%
Al2O3	0.09%	0.18%	0.19%	0.21%
MgO	1.32%	0.15%	0.08%	0.55%

) [16].

Siliceous minerals of the Dengying Formation in the Gaoshiti–Moxi area of the central Sichuan Basin show four types of quartz crystals (cryptocrystalline quartz, chalcedony, microcrystalline quartz, and megacrystalline quartz) and three structural types: cryptocrystalline, microcrystalline, and mosaic (laminar mosaic, window-hole mosaic, and interruptive and arcuate laminar mosaic) [2].

Siliceous minerals of different genesis in the dolomite reservoirs in the study area have different structural characteristics, and the contact relationship with the minerals and surrounding rocks is also different (Figure 1); this, in turn, affects the pore characteristics of the dolomite reservoirs. Siliceous minerals formed under secondary action generally have coarse grains and are mainly filled or accounted for in pore seams [16]; siliceous minerals selectively account for carbonatite grains formed by deposition, which are mainly manifested as fine-crystalline or microcrystalline, radial chalcedony, and cryptocrystalline siliceous minerals [17].

The storage performance of the reservoirs in the study area is negatively correlated with the siliceous mineral content in the formation. Formations with high siliceous mineral content tend to have relatively low porosity and poor reservoir performance (green box area in Figure 2) [18,19]. According to the diagram showing both the longitudinal wave velocity and the porosity of the formations in the study area, the reservoirs with high porosity and the formations with low porosity (high siliceous mineral content) are both characterized by low longitudinal wave velocity (green ellipse area in Figure 3) [19]. This indicates that it is not possible to distinguish between siliciclastic layers and reservoirs based on the zone of low longitudinal wave velocity in the study area; therefore, reservoir prediction in this area is difficult.

3. Petrophysical Modeling

Based on the characteristics of the siliciclastic dolomite reservoirs in the study area, they can be categorized into four parts in petrophysical modeling, namely, rock matrix, rock skeleton, mixed fluids in the rock, and saturated rock [3]. According to the aforementioned, the dolomite reservoir in the study area contains a large number of siliceous minerals and a small number of clay minerals. The siliceous minerals are of different origins, their structural characteristics are different, and there are differences in the pore spaces of the dolomite reservoir and its petrophysical characteristics, which, in turn, affects the accuracy of the subsequent transverse wave velocity estimation in the reservoir. To address this problem, this study, on the basis of the existing petrophysical model, makes some improvements to the method of obtaining the elastic modulus of the rock matrix to better describe the structural characteristics of siliceous minerals in dolomite reservoirs [20].

3.1. Calculation of the Rock Matrix Modulus of Elasticity

3.1.1. Commonly Used Models

When the wavelength of seismic waves is much larger than the non-uniform scale (particle scale) of rocks, the rocks can be treated as statistically uniform objects, and the concept of an equivalent medium is used to describe and characterize their properties [3,20]. Estimating the equivalent modulus for rock mixtures typically relies on the elastic modulus, volume fraction, and specific combination details of each mineral. In existing petrophysical models, siliceous minerals are usually considered part of the rock matrix, which is a mixture of various minerals, without taking into account their structural characteristics. However, as mentioned earlier, the siliceous minerals in the study area have different origins and exhibit diverse structural characteristics. Therefore, in the petrophysical modeling of silica-bearing dolomite reservoirs, special attention needs to be given to the minerals’ combination details. Typically, siliceous minerals with different structural characteristics are treated as inclusions with unique shapes within the dolomite. The equivalent aspect ratio of siliceous minerals is used to describe their structural characteristics, and the KT, DEM, and SCA models are employed to calculate the matrix elastic modulus of silica-bearing dolomite rocks [3]. Nevertheless, when applying the KT, DEM, and SCA models to determine the elastic modulus of the rock matrix, there are certain limitations associated with the equivalent aspect ratio of the pores (

Table 3. Comparison table of theoretical basis and applicable conditions of commonly used inclusion models.

	Theoretical Basis	Applicable Situation	Calculation Formula
KT Model	Based on the long-wavelength first-order scattering principlededuced for a two-phase medium full of fluid. Equivalent modulus of a fluid-saturated rock in a two-phase medium	Porosity to pore aspect ratio is much less than 1. Suitable for low porosity and low fracture density of rocks	$\begin{array}{c}\left({\mathrm{K}}_{\mathrm{K}\mathrm{T}}^{\ast }-{\mathrm{K}}_{\mathrm{m}}\right)\frac{3{\mathrm{K}}_{\mathrm{m}}+4{\mathrm{G}}_{\mathrm{m}}}{3{\mathrm{K}}_{\mathrm{K}\mathrm{T}}^{\ast }+4{\mathrm{G}}_{\mathrm{m}}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{v}}_{\mathrm{i}}\left({\mathrm{K}}_{\mathrm{i}}-{\mathrm{K}}_{\mathrm{m}}\right){\mathrm{P}}^{\left(\mathrm{m},\mathrm{i}\right)}\\ \left({\mathrm{G}}_{\mathrm{K}\mathrm{T}}^{\ast }-{\mathrm{G}}_{\mathrm{m}}\right)\frac{{\mathrm{G}}_{\mathrm{m}}+{\mathsf{\zeta }}_{\mathrm{m}}}{{\mathrm{G}}_{\mathrm{K}\mathrm{T}}^{\ast }+{\mathsf{\zeta }}_{\mathrm{m}}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{v}}_{\mathrm{i}}\left({\mathrm{G}}_{\mathrm{i}}-{\mathrm{G}}_{\mathrm{m}}\right){\mathrm{Q}}^{\left(\mathrm{m},\mathrm{i}\right)}\end{array}$ $\mathsf{\zeta }=\frac{\mathrm{G}}{6}\times \frac{9\mathrm{K}+8\mathrm{G}}{\mathrm{K}+2\mathrm{G}}$
DEM Model	By gradually adding the inclusion phase to the main mineral phase to simulate a biphasic mixture to obtain the equivalent modulus of the biphasic mineral equivalent modulus	Only one wrapper can be added at a time, and the result depends on the order of addition (the order of addition has no physical meaning)	$\begin{array}{c}\left(1-\mathrm{y}\right)\frac{\mathrm{d}\left({\mathrm{K}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast }\right)}{\mathrm{d}\mathrm{y}}=({\mathrm{K}}_{\mathrm{i}}-{\mathrm{K}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast }){\mathrm{P}}^{(\ast ,\mathrm{i})}\\ \left(1-\mathrm{y}\right)\frac{\mathrm{d}\left({\mathrm{G}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast }\right)}{\mathrm{d}\mathrm{y}}=({\mathrm{G}}_{\mathrm{i}}-{\mathrm{G}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast }){\mathrm{Q}}^{\left(\ast ,\mathrm{i}\right)}\end{array}$
SCA Model	Continuously adjust the substrate elastic parameters of the matrix until the plane wave incident on the porous medium no longer induces scattering by the incident plane wave. At this point, the elastic modulus of the porous medium can be equated to the elastic modulus of the medium. The modulus of elasticity of the matrix can be equated to the effective elastic modulus of the matrix	The pores are not interconnected, and the wavelength is much larger than the size of the inclusion. Suitable for rocks with large porosity	$\begin{array}{c}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{v}}_{\mathrm{i}}({\mathrm{K}}_{\mathrm{i}}-{\mathrm{K}}_{\mathrm{S}\mathrm{C}}^{\ast }){\mathrm{P}}^{(\ast ,\mathrm{i})}=0\\ {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{v}}_{\mathrm{i}}({\mathrm{G}}_{\mathrm{i}}-{\mathrm{G}}_{\mathrm{S}\mathrm{C}}^{\ast }){\mathrm{Q}}^{(\ast ,\mathrm{i})}=0\end{array}$

Note: In the calculation formula in the table, K_m and G_m are the moduli of the matrix mineral (K and G represent the bulk modulus and shear modulus, respectively); vi and y are the volume fraction of each inclusion;${\mathrm{K}}_{\mathrm{K}\mathrm{T}}^{\ast }$, ${\mathrm{G}}_{\mathrm{K}\mathrm{T}}^{\ast },{\mathrm{K}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast },{\mathrm{G}}_{\mathrm{D}\mathrm{E}\mathrm{M}}^{\ast },{\mathrm{K}}_{\mathrm{S}\mathrm{C}}^{\ast },{\mathrm{G}}_{\mathrm{S}\mathrm{C}}^{\ast }$ are each the equivalent modulus to be solved for each model. The DEM model and the SCA model are generally coupled because the P and Q coefficients contain the parameters to be solved, so they are all difficult to solve directly. In this paper, the fourth-order Runge–Kutta method is used for the DEM model to obtain its approximate solution, and for the SCA model, the approximate solution is obtained by an iterative algorithm with the initial value of VRH (Voigt–Reuss–Hill) averaged over the iterations. ${\mathrm{P}}^{\left(\mathrm{m},\mathrm{i}\right)}$, ${\mathrm{Q}}^{\left(\mathrm{m},\mathrm{i}\right)}$, ${\mathrm{P}}^{\left(\ast ,\mathrm{i}\right)}$, ${\mathrm{Q}}^{\left(\ast ,\mathrm{i}\right)}$ are the geometric factors related to the pore aspect ratio, reflecting the influence of the pore on the rock matrix.

). These limitations prevent an accurate depiction of the structural characteristics of siliceous minerals and hinder these models’ practical applicability [3].

3.1.2. Matrix Modulus of Silica-Bearing Dolomite

To address the above problems, in this study, we borrow the solid replacement method to better characterize the structure of siliceous minerals when finding the rock matrix modulus. The siliceous minerals are considered as inclusions into which empty pore space (the volume percentage occupied by siliceous minerals) is first added using a model, then silica is added using the solid replacement equation. The solid replacement equation given by Ciz and Shapiro [21] is a generalization of the anisotropic Gassmann equation to the case where elastic solids fill the pore space. This equation contains newly defined parameters related to the solids in the pore space [22], so, in this paper, a simplification is used, and the final solid replacement equation can be written as follows:

$${\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}^{-1}={\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-\frac{({\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-{\mathrm{K}}_{0\_\mathrm{m}}^{-1})}{{\mathsf{\phi }}_{\mathrm{m}}\left({\mathrm{K}}_{\mathrm{f}\_\mathrm{m}}^{-1}-{\mathrm{K}}_{0\_\mathrm{m}}^{-1}\right)+({\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-{\mathrm{K}}_{0\_\mathrm{m}}^{-1})}$$

(1)

$${\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}^{-1}={\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-\frac{({\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-{\mathrm{G}}_{0\_\mathrm{m}}^{-1})}{{\mathsf{\phi }}_{\mathrm{m}}\left({\mathrm{G}}_{\mathrm{f}\_\mathrm{m}}^{-1}-{\mathrm{G}}_{0\_\mathrm{m}}^{-1}\right)+({\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}^{-1}-{\mathrm{G}}_{0\_\mathrm{m}}^{-1})}$$

(2)

In the above, ${\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}$, ${\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}$ are the elastic modulus values of solid saturated rock (where K and G represent the bulk modulus and shear modulus, respectively); ${\mathrm{K}}_{0\_\mathrm{m}},{\mathrm{G}}_{0\_\mathrm{m}}$ are the modulus of elasticity values of the skeletal minerals that make up the solid rock; ${\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}},{\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}$ are the effective modulus of elasticity values of the solid rock skeleton; ${\mathrm{K}}_{\mathrm{f}\_\mathrm{m}}$, ${\mathrm{G}}_{\mathrm{f}\_\mathrm{m}}$ are the modulus of elasticity values of the solid inclusion; and ${\mathsf{\phi }}_{\mathrm{m}}$ is the volume fraction of solid inclusions. (The rock physical elastic parameters used in this article are shown in

Table 4. Mineral composition and petrophysical parameters of fluids in the dolomite reservoir of the Dengying Formation in the Gaoshitai area of central Sichuan.

Composition	Bulk Modulus (GPa)	Shear Modulus (GPa)	Density (g/cm3)
Dolomite	94.9	45	2.87
Quartz	37.9	43.7	2.65
Muddy	21	9	2.54
Water	2.25	0	1
Gas	0.13 × 10−3	0	0.65 × 10−3

.)

The siliceous dolomite in the study area is mainly composed of dolomite, siliceous minerals, and a few clay minerals. In the calculation using the solid replacement equation, dolomite is used as the rock skeleton mineral, and both siliceous minerals and clay minerals are added into the dolomite as solid inclusions. The structural morphology of siliceous minerals in the study area is closer to that of a coin-shaped seam, so in this paper, we use coefficients P and Q as shape factors of the inclusions in the coin-shaped seam to describe the structural characteristics of siliceous minerals (Equation (3)).

$$\left\{\begin{array}{c}{\mathrm{P}}^{\left(\mathrm{m},\mathrm{i}\right)}=\frac{{\mathrm{K}}_{\mathrm{m}}+\frac{4}{3}{\mathrm{G}}_{\mathrm{i}}}{{\mathrm{K}}_{\mathrm{i}}+\frac{4}{3}{\mathrm{G}}_{\mathrm{i}}+\mathsf{\pi }\mathsf{\alpha }{\mathsf{\beta }}_{\mathrm{m}}}\\ {\mathrm{Q}}^{\left(\mathrm{m},\mathrm{i}\right)}=\frac{1}{5}(1+\frac{8{\mathrm{K}}_{\mathrm{m}}}{4{\mathrm{G}}_{\mathrm{i}}+\mathsf{\pi }\mathsf{\alpha }\left({\mathrm{G}}_{\mathrm{m}}+2{\mathsf{\beta }}_{\mathrm{m}}\right)}+2\frac{{\mathrm{K}}_{\mathrm{i}}+\frac{2}{3}\left({\mathrm{G}}_{\mathrm{i}}+{\mathrm{G}}_{\mathrm{m}}\right)}{{\mathrm{K}}_{\mathrm{i}}+\frac{4}{3}{\mathrm{G}}_{\mathrm{i}}+\mathsf{\pi }\mathsf{\alpha }{\mathsf{\beta }}_{\mathrm{m}}})\end{array}\right.$$

(3)

In the above, $\mathsf{\beta }=\mathrm{G}\frac{3\mathrm{K}+\mathrm{G}}{3\mathrm{K}+4\mathrm{G}}$; $\mathsf{\alpha }$ is the aspect ratio of the coin inclusion and takes a value within the interval [0,1]; ${\mathrm{K}}_{\mathrm{m}}$ and ${\mathrm{G}}_{\mathrm{m}}$ are the modulus of elasticity values of the matrix mineral; and ${\mathrm{K}}_{\mathrm{i}}$ and ${\mathrm{G}}_{\mathrm{i}}$ are the modulus of elasticity values of the inclusion.

It has also been mentioned that all three models, KT, DEM, and SCA, can be used to incorporate empty porosity into the rock matrix. However, the KT and SCA models have certain requirements on the range when calculating the porosity; thus, the DEM model, which has better applicability, was chosen for the incorporation of solid inclusions (empty porosity) in the dolomite. To avoid the occurrence of error in the final results due to different sequences of the addition of solid inclusions, in this paper, we treat the siliceous minerals and small number of clay minerals as a whole as a kind of hybrid mineral with consistent equivalent aspect ratios and identical structural features.

The specific steps are as follows:

① The DEM model is used to add empty pore space with the same volume fraction (${\mathsf{\phi }}_{\mathrm{m}}={\mathrm{v}}_{\mathrm{s}\mathrm{i}}+{\mathrm{v}}_{\mathrm{s}h}$) as the siliceous minerals and clay and then obtain ${\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}$, ${\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}\_\mathrm{m}}$; the pore shape factors (P, Q) used in this step are for coin-shaped fractures (Equation (3)).

② The VRH model is used to calculate the mixed elastic modulus values of the siliceous minerals and clay and then obtain ${\mathrm{K}}_{\mathrm{f}\_\mathrm{m}}$, ${\mathrm{G}}_{\mathrm{f}\_\mathrm{m}}$.

③ Equations (1) and (2) are used to calculate the modulus of elasticity values of solid saturated rocks, ${\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}},{\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}$, which form the “new matrix modulus of elasticity of rocks characterized by siliceous minerals” described in this paper.

3.2. Calculating the Elastic Modulus of the Dry Rock Skeleton

Carbonate reservoir pores are mostly of the suture type, and the influence of pore type on the physical characteristics of the reservoir rock should be considered in the calculation. The DEM model also needs to consider the order in which pores are added, so the KT model is used to add empty pores to the rock matrix and then calculate the elastic modulus of the rock skeleton (${\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$,${\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$).

$$\left\{\begin{array}{c}\left({\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}-{\mathrm{K}}_{\mathrm{m}\_\mathrm{s}}\right)\frac{3{\mathrm{K}}_{\mathrm{m}\_\mathrm{s}}+4{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}}{3{\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}+4{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{w}}_{\mathrm{i}}\left({\mathrm{K}}_{\mathrm{i}}-{\mathrm{K}}_{\mathrm{m}\_\mathrm{s}}\right)\mathsf{P}\\ \left({\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}}-{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}\right)\frac{{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}+{\mathsf{\zeta }}_{\mathrm{m}\_\mathrm{s}}}{{\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}}+{\mathsf{\zeta }}_{\mathrm{m}\_\mathrm{s}}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}}{\mathrm{w}}_{\mathrm{i}}\left({\mathrm{G}}_{\mathrm{i}}-{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}\right)\mathsf{Q}\end{array}\right.$$

(4)

In the above, ${\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$ and ${\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$ are the modulus of elasticity values of the dry rock skeleton to be obtained; M is the number of inclusions; ${\mathsf{\zeta }}_{\mathrm{m}\_\mathrm{s}}=\left({\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}/6\right)\times \left(9{\mathrm{K}}_{\mathrm{m}\_\mathrm{s}}+8{\mathrm{G}}_{\mathrm{m}\_\mathrm{s}}\right)/\left({\mathrm{K}}_{\mathrm{m}\_\mathrm{s}}+2\mathrm{G}\mathrm{m}\_\mathrm{s}\right)$; and ${\mathrm{w}}_{\mathrm{i}}$ is the volume percentage of each pore.

In previous studies on petrophysical modeling, more emphasis was placed on describing the shape characteristics (aspect ratio) of pores. Combined with the pore characteristics of sutured carbonate reservoirs, the pores in this paper are broadly classified into hard pores (pore aspect ratio of 0.8 and volume fraction of ${\mathrm{w}}_{\mathrm{s}}$) and soft pores (pore aspect ratio of 0.02 and volume fraction of ${\mathrm{w}}_{\mathrm{c}}=1-{\mathrm{w}}_{\mathrm{s}}$).

3.3. Calculating the Modulus of Elasticity of Mixed Fluids

The bulk modulus of the mixed fluid is calculated using Wood’s equation.

$${\mathrm{K}}_{\mathrm{f}}={\left({\mathrm{S}}_{\mathrm{w}}/{\mathrm{K}}_{\mathrm{w}}+{\mathrm{S}}_{\mathrm{g}}/{\mathrm{K}}_{\mathrm{g}}\right)}^{-1}$$

(5)

In the above, ${\mathrm{K}}_{\mathrm{f}}$ is the mixed-fluid bulk modulus; ${\mathrm{S}}_{\mathrm{w}}$ and ${\mathrm{S}}_{\mathrm{g}}$ are the saturation of water and gas, where ${{\mathrm{S}}_{\mathrm{g}}=1-\mathrm{S}}_{\mathrm{w}}$; and ${\mathrm{K}}_{\mathrm{w}}$ and ${\mathrm{K}}_{\mathrm{g}}$ are the bulk modulus values of water and gas.

3.4. Calculating the Modulus of Elasticity of Saturated Rocks

The bulk modulus of the saturated rock is calculated using the Gassmann equation.

$$\frac{{\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{{\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}-{\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}}=\frac{{\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}}{{\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}-{\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}}+\frac{{\mathrm{K}}_{\mathrm{f}}}{\mathsf{\phi }({\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}\_\mathrm{m}}-{\mathrm{K}}_{\mathrm{f}})}$$

(6)

$${\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}}={\mathrm{G}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$$

(7)

In the above, ${\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ and ${\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ are the modulus of elasticity values of fluid-saturated rock; ${\mathrm{K}}_{\mathrm{d}\mathrm{r}\mathrm{y}}$ is the modulus of elasticity value of the rock skeleton; and $\mathsf{\phi }$ is the porosity of the rock.

4. Transverse Wave Velocity Estimation

The elastic parameters derived from the previous rock physics model were used to calculate the longitudinal and transverse wave velocities for saturated rocks.

$${\mathrm{V}}_{\mathrm{p}\mathrm{c}}=\sqrt{\frac{3{\mathrm{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}+4{\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{{\mathsf{\rho }}_{\mathrm{s}\mathrm{a}\mathrm{t}}}}$$

(8)

$${\mathrm{V}}_{\mathrm{s}\mathrm{c}}=\sqrt{\frac{{\mathrm{G}}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{{\mathsf{\rho }}_{\mathrm{s}\mathrm{a}\mathrm{t}}}}$$

(9)

In the above, ${\mathrm{V}}_{\mathrm{p}\mathrm{c}}$ is the calculated longitudinal wave velocity; ${\mathrm{V}}_{\mathrm{s}\mathrm{c}}$ is the calculated transverse wave velocity; and ${\mathsf{\rho }}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ is the density of the saturated rock, which can be obtained directly by density logging.

In the method described in this paper, two types of parameters—the equivalent aspect ratio $\mathsf{\alpha }$, which characterizes the mineral structure in the process of finding the matrix modulus of silica-bearing dolomite, and two volume fractions, ${\mathrm{w}}_{\mathrm{s}}$ and ${\mathrm{w}}_{\mathrm{c}}$, which characterize the pore structure of the rock in the process of finding the elastic modulus of the dry rock skeleton—cannot be obtained directly. At present, the conventional method is to obtain these parameters from the existing longitudinal velocity inversion calculation, and then substitute these parameters into the model to calculate the transverse wave velocity. In this paper, a relatively simple inversion method of preferentially iterating the pore volume fraction is used, followed by iterating the equivalent aspect ratio parameters of siliceous minerals to obtain a more accurate solution when the porosity is low or when other changes in the porosity volume fraction still do not yield an accurate solution.

In this paper, the simulated annealing iterative inversion method is used to obtain each parameter, and the specific process is shown in Figure 4, in which f = ∣v_pc − v_{p_measured}∣ and ${\mathrm{v}}_{\mathrm{p}\_\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{d}}$ denotes the longitudinal wave velocity values obtained from logging [23].

(1) The initial value ${\mathsf{\alpha }}_{\_0}$ (where subscript _t indicates the number of iterations of siliceous minerals and _0 indicates the initial value) can be estimated according to the actual situation of the reservoir, or by fixing the pore parameters and then quickly retrieving a value with the smallest error (in the retrieval, it can be set to 0.01, 0.05, 0.10, 0.20, 0.50, 0.75, or 0.99, and then substituted into the percentage of pore space at the last sampling point).

The initial value ${\mathrm{w}}_{\mathrm{c}\_0}$ (where the subscript _n indicates the number of pore iterations and _0 indicates the initial value, which can be considered as the same iteration parameter due to ${\mathrm{w}}_{\mathrm{s}}=1-{\mathrm{w}}_{\mathrm{c}}$) can be set to 50% or another empirical value. The initial value ${\mathrm{f}}_{\_0}$ is then calculated so that the current solution $\mathrm{f}={\mathrm{f}}_{\_0}$, and we continue to step (2).

(2) The perturbation generates a new solution ${\mathrm{w}}_{\mathrm{c}\_\mathrm{n}}$ (at this point, the number of pore iterations is incremented, $\mathrm{n}=\mathrm{n}+1$), calculates the new objective function ${\mathrm{f}}_{\mathrm{n}+\mathrm{t}}$, and continues with step (3). The parameter perturbation needs to be optimized. For example, when ${\mathrm{f}}_{\_\mathrm{n}+\mathrm{t}}$ is greater than $\mathsf{\alpha }$, the pore volume fraction is varied by 0.5%; otherwise, it is varied by 0.01%. The value of $\mathsf{\alpha }$ should be taken to ensure that when the volume fraction is varied by 0.5%, the longitudinal wave velocity variation is less than the value of $\mathsf{\alpha }$.

(3) Determine whether ${\mathrm{f}}_{\_\mathrm{n}+\mathrm{t}}$ is smaller than the current solution $\mathrm{f}$. If yes, continue to step (4); otherwise, continue to step (5).

(4) Accept ${\mathrm{w}}_{\mathrm{c}\_\mathrm{n}}$ and ${\mathsf{\alpha }}_{\_\mathrm{t}}$ as the new solutions for ${\mathrm{w}}_{\mathrm{c}}$ and $\mathsf{\alpha }$, and then continue to step (6). When $\mathsf{\alpha }$ is unperturbed, ${\mathsf{\alpha }}_{\_0}$ is the solution.

(5) Accept the new solution according to the Metropolis criterion when $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}<\mathrm{e}\mathrm{x}\mathrm{p}\left[-\left({\mathrm{f}}_{\_\mathrm{n}+\mathrm{t}}-\mathrm{f}\right)/{\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{t}\right]$ (where rand is a random number) and accept ${\mathrm{w}}_{\mathrm{c}\_\mathrm{n}}$ and ${\mathsf{\alpha }}_{\_\mathrm{t}}$ as the solutions for ${\mathrm{d}\mathrm{w}}_{\mathrm{c}}$ and $\mathsf{\alpha }$, respectively; otherwise, maintain the original solution and continue to step (6).

(6) Determine whether n is less than ${\mathrm{n}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (the upper limit of the number of pore iterations). If yes, go back to step (2); otherwise, continue to step (10).

(7) Generate a new perturbation solution ${\mathsf{\alpha }}_{\_\mathrm{t}}$ (increment the number of pore iterations, $\mathrm{t}=\mathrm{t}+1$), calculate the new objective function ${\mathrm{f}}_{\_\mathrm{n}+\mathrm{t}}$, and continue to step (8).

(8) Determine whether t is less than ${\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (the upper limit of the number of mineral iterations). If so, go back to step (3); otherwise, continue to step (9).

(9) Make ${\mathsf{\alpha }}_{\_\mathrm{t}}=\mathsf{\alpha }$ (i.e., make the current perturbed ${\mathsf{\alpha }}_{\_\mathrm{t}}$ reduce to the current optimal solution), set the number of pore iterations to $\mathrm{n}=0$ (i.e., reset the number of pore iterations), and then go back to step (2).

(10) Determine whether the current solution is less than the given minimal value $\mathrm{e}$. (Here, $\mathrm{e}$ is set according to the accuracy requirement; it is recommended to set it to 10~50 m/s. The higher the accuracy, the smaller the perturbation interval of ${\mathrm{w}}_{\mathrm{c}}$ that needs to be set.) If yes, continue to step (11); otherwise, go back to step (7).

(11) End of the operation. Return the solutions for ${\mathrm{w}}_{\mathrm{c}}$ and $\mathsf{\alpha }$, calculate ${\mathrm{v}}_{\mathrm{s}\mathrm{c}}$, and output the result.

5. Results

The transverse velocity of the GS1 well in the Gaoshiti–Moxi area in the central Sichuan Basin was estimated using the improved model proposed in this paper. The estimated results are compared with the measured transverse wave velocity of the GS1 well and the predicted results from the Xu–Payne and Zhang Bingming models (Figure 5).

In Figure 5, within the whole range of the target layer section, formations with higher siliceous mineral content have lower porosity and poorer storage performance, which is consistent with the relationship between siliceous minerals and reservoirs summarized in Table 1 and the results of changes in siliceous mineral content and porosity in Figure 2. According to the measured longitudinal wave velocity curve, the phenomenon of low longitudinal wave velocity values is present not only in formations with high porosity but also in formations with low porosity; this result is consistent with the results presented in Figure 3.

As far as the constrained longitudinal velocities are concerned, the results of the Xu–Payne model are not sufficiently precise when compared to the measured results, especially in the formations with high silica content. On the other hand, the Zhang Bingming model and the model proposed in this paper obtained velocities consistent with the measured velocities. Only when the porosity of the formation is 0, both of them show a slight error.

Regarding the estimation of the transverse wave velocity, the errors between the measured velocities and those predicted by the Xu–Payne model are small in the strata with fewer siliceous minerals at the top and bottom of the target section, but larger in the strata with more siliceous minerals in the middle part of the section, while the errors between the measured velocities and those predicted by the Zhang Bingming model are small in the strata with higher siliceous mineral contents, and larger in the other strata. However, the improved model predicts the transverse wave velocity in the whole target section, and the error between the prediction and the measured velocity is smaller across the whole target section. The transverse wave velocity curves of this model better fit the measured transverse wave velocity curves, and the prediction results are closer to the transverse wave velocity measurements, making this model better than the previous two methods.

The prediction deviation is the absolute value of the difference between the actual measured transverse wave velocity values and the transverse wave velocity values predicted by the model. By statistically analyzing the number of data points in each interval, we can assess the degree of deviation between the predicted and measured values. A comparison of the statistical plots and the fitted plots reveals that the deviation intervals of the improved model mainly occur in a smaller range of values as compared to those of the Xu–Payne model and the Zhang Bingming model. In addition, the deviation between the predicted transverse wave velocity and the measured transverse wave velocity is the smallest for the improved model, among the three models compared (Figure 6).

6. Discussion

In this paper, the calculation method to obtain the matrix modulus of siliceous dolomite was improved by solid replacement, and the improved method is better able to describe the structural morphology of siliceous minerals in a reservoir. On this basis, the transverse wave velocity of the dolomite reservoir was estimated using the improved model in combination with other petrophysical models. According to the curve fit and the degree of deviation, the results of the transverse wave velocity estimation are good and obviously better than those from petrophysical models of carbonate rocks proposed by previous authors, which provides a certain basis for subsequent reservoir prediction work in this area.

The siliceous minerals in this area have three structural types: cryptocrystalline, microcrystalline, and mosaic structure (laminated mosaic structure, discontinuous window-shaped pore mosaic structure, and arc-shaped laminated mosaic structure). According to a combined plot of the siliceous mineral content and porosity in the study area, the storage performance of the reservoirs is negatively correlated with their siliceous mineral content, which echoes the results of many previous studies [24,25,26,27,28,29]. Reservoirs with higher porosity and siliceous layers with lower porosity (higher siliceous content) in the study area both exhibit lower longitudinal wave velocities, suggesting that velocity variations in carbonate reservoirs are sensitive to the siliceous mineral content [30,31]. Therefore, the effect of siliceous minerals on the porosity and wave velocity cannot be ignored when estimating the transverse wave velocity of reservoirs in this area using petrophysical models. However, the complex siliceous minerals in the reservoir cannot be accurately described when predicting the transverse wave velocity in the region using existing petrophysical models [7,8,9,10,11,12]. Therefore, in this paper, the rock matrix modulus was calculated by adding solid replacement to an existing model describing siliceous minerals. The improved model was utilized to predict the transverse wave velocity, and the prediction results were better than those from petrophysical models proposed by previous authors, such as the Xu–Payne model.

In this study, we focused on siliceous carbonate formations. However, when the local lithology changes, the mineral compositions in the formation change and the mineral and pore structures of the formation become more complex; some changes in these factors can cause changes in the elastic modulus of the formation. Therefore, in the face of the above situation, we should clarify the influence of minerals on the pore structure and physical parameters of stratigraphic rocks and seek suitable methods to characterize them. This will be a key step in constructing the whole petrophysical model.

7. Conclusions

The main purpose of this study was to solve the problem of estimating the transverse wave velocity of siliceous dolomite reservoirs of the Dengying Formation in the Gaoshiti–Moxi region of central Sichuan. In this study, the siliceous minerals were regarded as solid inclusions, and an improved model was proposed by using methods such as the solid replacement equation. The transverse wave velocity of the reservoir was estimated using the improved model. According to the shape and deviation of the velocity curves, it was shown that the transverse wave velocity curves estimated by the improved model fit well with the measured transverse wave velocity curves, with minimum deviation, and the prediction results were obviously better than those of petrophysical models of complex carbonate rocks proposed by previous authors. This method provides a reference for petrophysical models with complex mineral structure characteristics and provides a basis for the future identification of reservoirs in low-impedance areas. The applicable object of this study is siliceous dolomite reservoirs. In cases where the reservoir lithology changes and the reservoir minerals become more complex, the relationship between the reservoir mineral characteristics and the reservoir rock’s physical parameters should be carefully analyzed.