Fluid Discrimination in Ultra-Deep Reservoirs Based on a Double Double-Porosity Theory

We develop a methodology, based on rock-physics templates, to effectively identify reservoir fluids in ultra-deep reservoirs, where the poroelasticity model is based on the double double-porosity theory. P-wave attenuation, the ratio of the first Lamé constant to mass density (λ/ρ) and Poisson ratio are used to build the templates at the ultrasonic and seismic frequency bands to quantitatively predict the total and crack (soft) porosities and oil saturation. Attenuation on these frequency bands is estimated with the spectral-ratio and frequency-shift methods. We apply the methodology to fault-controlled karst reservoirs in the Tarim Basin (China), which contain ultra-deep hydrocarbon resources with a diverse pore-crack system, low porosity/permeability and complex oil-water spatial distributions. The results are consistent with well-log data and actual oil recovery. Crack porosity can be used as an indicator to find regions with high oil saturation, since high values implies a good pore connectivity.


INTRODUCTION
With the depletion of relatively shallow conventional reservoirs, ultra-deep resources will become an important target for oil and gas exploration. Carbonate reservoirs are typical examples (e.g., Zou et al., 2014), which, due to their burial depth and high degree of diagenesis, their pay zone is located in areas with secondary pores or cracks (or rift caves and karst caves) (Pang, 2010). The tectonic action, abnormal high pressure, the presence of fractures and thermal fluid activity favor the development of this secondary pore space. Fracture-cavern carbonate reservoirs are among the most important types (Tian et al., 2017;Li et al., 2019). A "fault-controlled karst reservoir" is a largescale fracture-cavern carbonate reservoir formed by a series of geological processes including the dissolution of atmospheric fresh water or hydrothermal fluid that flows along faults and associated complex fracture system under unexposed conditions (Lan et al., 2015).
The fabric of carbonate rocks has at least two dissimilar porosities, leading to a heterogeneous distribution of immiscible fluids or patchy saturation (Ba et al., 2017). Early models consider a single porosity, such as that of White (1975), who assumed patchy saturation with spherical gas inclusions, and a generalization to non-spherical patches by Johnson (2001). Then, Pride et al. (2004) studied attenuation caused by patchy saturation based on a double-porosity theory, while Dvorkin and Nur (1993) combined the Biot and squirt flow mechanisms into a set of poroelasticity equations, explicitly considering attenuation due to global flow at cracks or grain contacts. The model has been applied to predict wave velocity and attenuation (Dvorkin et al., 1994(Dvorkin et al., , 1995. Deng et al. (2015) extended the squirt-flow model and Tang (2011) analyzed the effect of fractures on the elastic waves. Ba et al. (2011) introduced the Biot-Rayleigh theory to model wave propagation in double-porosity media with a single saturating fluid, and Ba et al. (2017) extended this model to the case of two fluids, i.e., a double double-porosity medium model (DDP).
Rock-physics templates (RPT) relate the reservoir and seismic properties (Odegaard and Avseth, 2004;Carcione and Avseth, 2015). For instance, Chi and Han (2009) interpreted the lithology and obtained fluid content by using RPTs and logging data, and RPTs were used to compute reservoir porosity and oil/gas saturation by Michel (2010). Picotti et al. (2018) analyzed the effects of different fluid saturations, porosities, and permeabilities on seismic waves, and introduced an attenuation RPT for sandstone reservoirs. Pang et al. (2020) estimated the microfracture porosity in carbonate reservoirs by using these attenuation templates. Tran et al. (2020) predicted organic matter content of source rocks using RPTs.
Building a RPT implies to choose the proper properties sensitive to lithology or fluid type and saturation (Batzle et al., 2001;Dillon et al., 2003;Qiao and An, 2007). Goodway et al. (1997) proposed to identify and detect reservoir fluids by using the Lamé constants. Hedlin (2000) proposed the pore modulus as a fluid recognition factor through experiments to distinguish shale, calcite, gas sand, oil sand, and wet sand. The Gassmann fluid term and its product with mass density were used as a fluid discriminator by Russell et al. (2003Russell et al. ( , 2006, and Quakenbush et al. (2006) used Poisson's ratio for the same purpose. Quintal (2012) suggested that frequency-dependent attenuation due to waveinduced fluid flow is a potential index of saturation, and Xue et al. (2017) used a combination of the normal crack compliance of the saturated and dry rocks as a fluid factor.
We choose P-wave attenuation, Poisson's ratio and the ratio of the first Lamé constant to mass density (λ/ρ) to build the templates at ultrasonic and seismic frequencies, based on the DDP theory, and associate the total and crack porosities and oil saturation to the seismic properties. The theory considers the effects of fabric heterogeneity and heterogeneous distribution of immiscible fluids. Ultrasonic, well-log and seismic data are used to calibrate the templates. Then, these are applied to quantitatively predict the rock properties and oil saturation.

Overview of the Work Area
The area is located at the Tarim Basin, China, and the reservoir section is constituted by the upper members of the Yijianfang to Yingshan formation (Ordovician), formed during the Caledonian-Hercynian period. It is a marine fault-controlled karst reservoir (Li et al., 2019). The multi-stage fault activity has a controlling effect on the accumulation and distribution of oil and gas in the reservoir. Fractures and fluid transformation FIGURE 1 | (A) Schematic of a fault-controlled karst reservoir (orange represents caverns and vugs, while black refers to cracks); (B) thin section of a carbonate sample; (C) schematic of a double-porosity structure with patchy saturation of fluids (two skeletons and two fluids; orange represents oil and blue is water). have formed a diverse reservoir system including dissolved pores, inter-crystalline pores and cracks, which make the structure quite complex (Lan et al., 2015). The reservoir along the fault zone is oil-bearing and unevenly enriched. It is characterized by an ultra-deep burial (>7 km), high temperature (>150 • C) and high pressure (>160 MPa). Figure 1A shows a diagram of the reservoir, where the orange color represents caverns and vugs, and the black lines represent cracks. Figure 1B shows a thin section of a rock sample collected at a depth of 7443 m. The lithology is limestone, with a dolomite content of 6% and a calcite content of 94%. Dissolution of pores and multiple cracks of various diameters are observed and the surface porosity is extremely low (less than 1%).

Rock-Physics Modeling Flowchart
The reservoir characteristics are simplified with a doubleporosity structure with cracks and pores and saturation with two immiscible fluides, oil and water ( Figure 1C). Figure 2 shows the flowchart of the modeling methodology, which can be summarized as follows.
(1) The minerals are mainly calcite and a small amount of dolomite and clay. The Voigt-Reuss-Hill (VRH) equation (see Supplementary Appendix Eq. A1) is used to calculate the elastic moduli of the mineral mixture (Voigt, 1910;Reuss, 1929;Hill, 1952). The bulk and shear moduli of calcite, dolomite, and clay are 63.7 and 31.7, 76.4 and 49.7, and 21 and 7 GPa, respectively (Mavko et al., 2009). (2) The pores and cracks are assumed spherical and oblate in shape (Shapiro, 2003), with aspect ratios of 1 and 0.00035 (generally in range of 10 −5 ∼10 −4 , according to Tan et al., 2020 andPimienta et al., 2015), respectively. The connectivity of these two kinds of pores affects P-wave attenuation and dispersion, as well as the elastic parameters (Agersborg et al., 2008). In this model, the two pore systems are connected and the composite medium is isotropic. Based on the self-consistent approximation (SCA) (Berryman, 1980(Berryman, , 1992(Berryman, , 1998, the pores and cracks are treated as inclusions and mixed with the minerals to obtain two partial dry-rock elastic moduli and then mixed again with the SCA model to obtain the dry-rock elastic moduli. At this point, the frame (or skeleton) is a double-porosity structure. The relevant equations are given in Supplementary Appendix Eq. A2.
(3) The volume content of the host medium v 1 and the volume content of inclusion v 2 satisfy v 1 +v 2 =1. The total porosity φ = φ 10 · v 1 + φ 20 · v 2 , where φ 10 is the porosity of the host medium and φ 20 is the porosity of the inclusions, and the crack porosity is φ c = φ 20 · v 2 . φ 20 is 0.006 (in range of 0.005∼0.1, by considering that a larger cracked grain is modeled as a cylinder of axial length around 100 µm, containing a transverse crack with an average effective aperture of several microns, according to Pride et al., 2004;Pang et al., 2020, andBa et al., 2017). The stiff porosity is The equations proposed by Batzle and Wang (1992) are used to estimate the bulk moduli, densities, and viscosity of the fluids at 180 • C and 80 MPa. The bulk modulus, density, and viscosity of oil and water in the oil-water mixture are 0.79 and 2.2 GPa, 0.70 FIGURE 4 | P-wave velocity (A) and attenuation (B) as a function of frequency for varying crack (soft) porosity. and 0.93 g/cm3, and 0.002289 and 0.000647 Pa s, respectively. Then, fluid substitution is performed by using the DDP theory to obtain the acoustic properties. The main skeleton is composed of two components, namely, a host one with dissolved pores and an inclusion one with cracks, and each of them is patchy saturated. Each component can be considered as a secondary doubleporosity structure, i.e., two pore phases saturated with different fluids (there are four types of pores globally). By performing a plane-wave analysis (see Supplementary Appendix Eq. A5), we obtain the P-wave velocity and attenuation.

P-Wave Phase Velocity and Attenuation
The stiff porosity, crack porosity, and oil saturation are basic factors related to oil production in reservoirs. First, we study the effects of these properties on P-wave velocity and attenuation in the frequency range of 10 1 ∼10 8 Hz. Attenuation is represented by the dissipation factor, defined as 1000/Q, where Q is the quality factor. Oil saturation is 100%, stiff porosity is 6%, crack porosity is 0.06%. In the analysis, we vary one property and set the others constant. The fluid patch radius of the host medium is 50 µm, the inclusion radius is 50 µm (at pore/grain scale), and the fluid patch radius in the inclusion is 5 µm. Figures 3A,B show the effect of the stiff porosity on the P-wave phase velocity and dissipation factor, respectively, where it can be seen that velocity dispersion and attenuation increase with increasing porosity and the relaxation peaks move slightly to the low frequencies. Figure 4 shows that increasing crack porosity, the velocity dispersion becomes stronger, attenuation increases and the peaks move to the high frequencies, while Figure 5 indicates that dispersion and attenuation increase with increasing oil saturation, and the peaks move to the low frequencies.

Fluid Identification Factors
Poisson's ratio and λ/ρ are used to identify the fluids (Goodway et al., 1997;Quakenbush et al., 2006). Figures 6, 7 show the effects of stiff porosity, crack porosity, and oil saturation at 1 MHz, respectively. As can be seen, increasing stiff porosity and oil saturation, λ/ρ gradually decreases, while increasing crack porosity has the opposite effect. On the other hand (Figure 7),  Poisson's ratio increases with increasing stiff and crack porosities and decreases with increasing saturation.

Rock-Physics Templates
Total and crack porosities and oil saturation are set as variables to build a 3D RPT at 1 MHz, which is displayed in Figure 8A, where the parameters are the same as in section "P-Wave Phase Velocity and Attenuation." The blue and red line refers to water and oil saturations, respectively, and the behavior of the template agrees with the previous results. On the other hand, Figure 8B shows the 3D RPT at 35 Hz, where the radii of the inclusion, fluid patch in the host medium and fluid patch in the inclusion are 50, 50, and 5 mm, respectively. The characteristics of the seismic template are identical to those of the ultrasonic template.

Ultrasonic Band
Ultrasonic experiments were performed on five water-and oilsaturated samples at 1 MHz, a pore pressure of 70 MPa and 140 • C, whose properties are given in Table 1. The carbonate  samples used in the ultrasonic measurements are mainly made of limestone, and parts of cracks are filled with asphalt, calcite, and pyrite. Moreover, there are micro-cracks and a relative amount of dissolved pore cavities. In the experiments, the sample is fully saturated, jacketed, and put into a container according to the experimental setup of Guo et al. (2018a). The first arrivals are picked to obtain the wave velocities. Figure 9A,B show the velocities and the attenuation factor estimated with the spectral-ratio method (Toksöz et al., 1979;Guo and Fu, 2006; see Supplementary Appendix Eq. A3). Figure 9A shows that P-wave velocity at full water saturation is higher than that at full oil saturation, while the S-wave velocity FIGURE 9 | (A) Experimental P-S wave velocity crossplot; (B) P-wave dissipation factor as a function of the total porosity; (C) 3D RPT at 1 MHz compared to the experimental data. The blue and red dots refer to the samples saturated with water and oil, respectively. exhibits the opposite behavior. Figure 9B indicates that the oilsaturated samples have a stronger P-wave attenuation. As shown in the template of Figure 9C, the experimental data of these five samples are used for the calibration at ultrasonic frequencies, where the blue and red circles refer to the water-and oil-saturated samples, respectively. The average error between the predicted total porosities (based on ultrasonic wave attributes) and the measured porosities under oil saturation is 0.00344. The template shows good agreement with the experimental measurements.

Seismic Band
Based on well-log and seismic data, the RPT is calibrated at 35 Hz. The seismic quality factor is estimated with the improved frequency-shift method (Pang et al., 2019(Pang et al., , 2020, as shown in Supplementary Appendix Eq. A4. Figures 10A,B show the post-stack seismic data and attenuation (the dissipation factor 1000/Q) profiles of the target layer, respectively, where the red and blue colors indicate high and low attenuation ( Figure 10B). The results show that the rocks at the two wells exhibit significant attenuation. On the other hand, Poisson's ratio and λ/ρ are obtained from AVO inversion (Guo et al., 2018b;Luo et al., 2020a,b). The pre-stack angle gather data contains useful information about the amplitude change with angle. The threeparameter inversion method of pre-stack seismic data, based on the AVO theory, can directly estimate P-wave impedance, S-wave impedance and density, and from these the Poisson's ratio and λ/ρ are obtained. Figures 10C,D show 2D profiles of λ/ρ and Poisson's ratio obtained with pre-stack seismic inversion. We can see that Poisson's ratio is low in regions with a high λ/ρ. Attenuation, Poisson's ratio and λ/ρ at well A were used to calibrate the RPT at seismic frequencies (see Figure 11) and the log porosity of well A is assumed. Comparing the seismic inversion data with the RPT, we can observe that the seismic data   with porosity is consistent with the template. According to the color bar of Figure 11, the porosity of the seismic data is low (less than 7%), which agrees with the characteristics of the ultra-deep reservoirs. Attenuation, Poisson's ratio, and λ/ρ can be obtained from seismic data, and a quantitative prediction of total porosity, crack porosity, and oil saturation can be achieved by overlapping the data on the template.

PREDICTION OF THE RESERVOIR PROPERTIES
The calibrated RPTs are applied to estimate the reservoir properties. Figure 12 shows the reservoir total and crack porosities, and oil saturation on a 2D survey line crossing wells A and B. In Figure 12A, total porosity is generally less than 10% and areas around the wells exhibit higher porosity. As shown in Figure 13, the predicted results of well A are consistent with the porosity log profile. The average error between the predicted porosities (based on seismic wave attributes) and the measured porosities of well A is 0.0125. By combining Figures 12B,C, high crack porosity enables good connectivity and better oil storage capacity. In Figure 12, we also can see that the upper part of well B has a pore-crack system with low oil saturation, while the lower part has high oil saturation, indicating that the reservoir is not homogeneous. Comparative analysis shows that the predicted results are consistent with the reported oil production of well A (90.26 tons/day) and well B (50.31 tons/day). The prediction results of well A are used to analyze the correlation between crack porosity and oil saturation. The correlation coefficient is 0.6296 between fracture porosity and oil saturation.

Figure 14
gives horizontal sections of total and crack porosities and oil saturation, showing that well A exhibits higher porosity and oil saturation than well B, while the difference in crack porosity is small. The storage space of fault-controlled karst reservoirs is formed by the dissolution of water or hydrothermal fluid along the fault and the distribution of the fluid is extremely uneven. Combining Figures 14A,B, the distributions of total and crack porosities are complementary and the oil saturation distribution in Figure 14C is irregular. These are typical characteristics of fault-controlled karst reservoirs.
It can be seen in Figures 12, 14 that fault-controlled karst reservoirs do not have a defined oil-water interface, but a patchy fluid saturation due to dissimilar crack-pore characteristics. Crack porosity can be used as an indicator to find regions with high oil saturation. Stiff pores mainly store oil and soft pores act as flow channels and also play an important storage role. These findings are consistent with the actual characteristics of the reservoirs, indicating that the proposed method can be successfully applied.

CONCLUSION
We have implemented a double double-porosity theory to build RPT in order to estimate the properties of deep carbonate reservoirs. The dry-rock moduli were obtained with a combination of the VRH equation and the self-consistent model. The templates axes involve attenuation, Poisson's ratio and the ratio Lamé constant to density. Ultrasonic experiments were performed on rock samples of the study area and used to calibrate the templates, as well as information from well logs. Then, we have predicted total porosity, crack porosity and oil saturation of the carbonate reservoirs. According to the results, the reservoir has low porosity with developed cracks, which is consistent with the local geological features. Moreover, the results are basically consistent with the actual oil production in the area.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, and further inquiries can be directed to the corresponding author/s.