Abstract
Secondary pores are the main reservoir space and transportation channel of oil and gas in reef limestone reservoir. At present, the main method of calculating secondary porosity is restricted by the morphological characteristics of porosity spectrum, regional artificial influence, and accuracy of calculation. We present a new method for calculating secondary porosity of reef limestone reservoir by the nuclear magnetic resonance spectrum which is calibrated by casting thin section. We begin with analyzing and determining the high correlation between the surface porosity of casting thin section and the total porosity. The objective is confirming the feasibility of the method of calculating secondary porosity by using thin-section information calibrate. Then, we use the surface porosity of thin section as the calibrating data and find the relaxation time corresponding to the best correlation between the secondary porosity and the secondary surface porosity of casting thin section, that is, the cutoff value of secondary porosity, through the Monte Carlo method. Finally, we calculate the secondary porosity by using the functional relationship between the secondary surface porosity and the surface porosity. The statistical analysis shows that the method of calculating secondary porosity effectively improves the calculation accuracy of secondary porosity. The secondary porosity calculation results have a high correlation with the reservoir productivity.
1. Introduction
Reef limestone reservoir is one of the most important types of carbonate reservoirs [1–4] and plays an extremely important role in world oil reserves [5]. Reef limestone reservoir has both primary and secondary pore, which is characterized by diverse reservoir space, complex pore structure, and fast phase transition [6–10]. Primary pore and secondary pore in reef limestone reservoir have good reservoir and permeability performance, but because of diagenesis and other factors, the primary pore is often difficult to effectively preserve, and the secondary pore becomes the main reservoir space and migration channel of oil and gas [11–15]. Therefore, there is an increasing demand for secondary porosity evaluation in reef limestone reservoir in both geological exploration and reservoir development. In particular, the accurate secondary porosity evaluation conclusion can provide basic data and theoretical basis for the formulation and adjustment of oil field development plan and the effective improvement of oil recovery.
Geophysical logging, as an important technical means and data source for evaluating reservoir porosity, undoubtedly becomes the preferred technology for quantitative evaluation of secondary porosity of reef limestone reservoir. However, because the formation of secondary porosity is controlled by a series of factors such as denudation, cementation, diagenesis, and biological perturbation [16–19], the response relationship between single conventional logging parameters and secondary porosity is not ideal. Using the conventional logging data cannot effectively evaluate secondary porosity [20]. Compared with conventional logging data, electrical imaging logging data has prominent advantages in intuition, directivity, and high resolution. Therefore, a group of geophysicists and workers represented by Newberry [21, 22] have begun to use microresistivity imaging data to evaluate secondary porosity of reef limestone reservoir. At present, the main method of evaluating reservoir secondary porosity by logging is to calculate secondary porosity by means of electrical imaging logging porosity spectrum. In essence, the electrical imaging porosity spectrum is essentially a frequency histogram formed by frequency statistics in different intervals of porosity values for a certain depth of reservoir. Generally, the evaluation process is to calculate the reservoir porosity spectrum by using electrical imaging logging data, then calculate the cutoff value of the porosity spectrum, and then calculate the secondary porosity of the reservoir.
The cutoff value of porosity spectrum is a key parameter to distinguish primary and secondary porosity in the porosity spectrum, which directly affects the accuracy of the calculation results of secondary porosity. At present, the methods to determine the cutoff value of porosity spectrum in electrical imaging logging include Newberry’s cutoff value calculation method based on median porosity proposed by Newberry of Schlumberger Company [21], the method for calculating the cutoff value of discriminant porosity spectrum based on standard linear discriminant analysis proposed by Ramakrishnan et al. [23], and the method of calculating the cutoff value of porosity spectrum fitting by Gauss function based on the shape of porosity spectrum proposed by Shi et al. of China University of Geosciences (Beijing) [24]. Among them, the accuracy of Gauss function fitting method is closely related to the shape of porosity spectrum, which is not suitable for the case of large bimodal spacing in porosity spectrum, and its applicability in reef limestone reservoir is limited. Newberry method and Ramakrishnan method are both mathematical methods based on porosity. Taking Newberry method as an example, the cutoff value of porosity spectrum is obtained by calculating the variance of porosity in the section of minimum porosity and median porosity. The calculation formula is as follows: where means the porosity variance between the minimum porosity and the median porosity and means the fixed coefficient of Newberry method.
It can be seen from Equation (1) that Newberry’s cutoff value method of porosity spectrum can reflect reservoir’s porosity characteristics, and this method is not affected by the morphological characteristics of porosity spectrum. Compared with the calculation method of Gauss function fitting, it has stronger applicability. However, the accuracy of calculating the cutoff value of porosity spectrum depends greatly on the value of fixed coefficient , but its value is usually based on regional empirical value, which lacks a theoretical basis; therefore, the accuracy of secondary porosity calculation is restricted.
In order to improve the calculation accuracy of the secondary porosity, we propose a method of calculating secondary porosity of reef limestone reservoir by casting thin section to calibrate NMR spectrum. At present, nuclear magnetic resonance logging and nuclear magnetic resonance petro physics experiment have been widely used in reservoir geophysical logging evaluation. Because the pore size distribution of reservoir is directly related to the shape of nuclear magnetic resonance spectrum, the content of different components can be regarded as the number of different pore size components. Therefore, the multicomponent pore can be determined by different values, and the pore content of corresponding pore can be obtained by accumulative value of pore size components. The method we propose takes the advantage of nuclear magnetic technology. Firstly, the secondary surface porosity of using casting thin-section recognition as the calibration data and the corresponding relationship between secondary surface porosity of casting thin sections and the pre- value of NMR is found by Monte Carlo method. Secondary porosity cutoff value of NMR is determined by this relationship. Finally, reservoir secondary porosity is calculated according to the cutoff value of determined by casting thin-section data calibration. This method not only makes full use of the advantages of NMR data in pore structure analysis but also takes into account the calibration of the cutoff value of NMR by casting thin-section surface porosity. Compared with the determination of the cutoff value of secondary porosity spectrum by electrical imaging, it has more theoretical basis and further ensures the accuracy of secondary porosity calculation.
2. Materials and Methods
2.1. Calculating Secondary Porosity of Reef Limestone Reservoir by NMR Spectrum
2.1.1. The Principle of Discriminating Different Types of Pores by Nuclear Magnetic Resonance Method
At present, there are many methods to distinguish pore types by analyzing reservoir pore structure. In addition to traditional casting thin sections, scanning electron microscopy (SEM) and mercury injection method, advanced testing techniques such as field emission scanning electron microscopy (FSEM), focused ion beam microscopy (FIB), nuclear magnetic resonance (NMR), and micron-nanometer CT scanning have gradually begun to be used in pore structure research. Among these methods, nuclear magnetic resonance core technology is a new method developed rapidly in recent years. It has the advantages of fast, nondestructive, and repeatable [25]. It can realize high-precision measurement of micron-nanometer pore in different types of reservoir rocks. It can also meet the needs of large-scale and continuous testing in reservoir evaluation [26]. These advantages make nuclear magnetic resonance technology have great application prospects in pore type discrimination research.
According to the relaxation mechanism of NMR and the principle of rock physical property measurement, when the saturated water rock is in a uniform magnetic field, the lateral relaxation time of NMR is proportional to the pore throat radius . where means the conversion coefficient between and .
As long as the value is determined, the pore radius distribution curve of rock can be obtained by using nuclear magnetic resonance spectrum, which is the theoretical basis of quantitative evaluation of rock pore structure by nuclear magnetic resonance technology. The magnitude of value depends on the lateral surface relaxation rate and the pore shape factor . Because these two parameters of rock cannot be measured directly in practical research, the conversion coefficient value can be obtained by calibrating the NMR spectrum with more accurate mercury injection data based on the correlation between the NMR spectrum and the mercury injection curve. However, the transformation method proposed by predecessors can achieve high-precision NMR measurement of pore throat radius distribution in medium-high porosity and permeability sandstone reservoirs, while the problem of transformation between NMR spectrum and pore throat radius distribution in reef limestone reservoirs has not been solved well.
2.1.2. The Principle of Calculating Secondary Porosity of Reef Limestone Reservoir by NMR Spectrum
Different relaxation time in NMR spectrum corresponds to different pore radius sizes. Based on the assumption of normal distribution model, for homogeneous reservoirs, rock NMR spectrum shows unimodal characteristics. For highly heterogeneous reservoirs, such as reef limestone reservoirs, primary and secondary pores are developed simultaneously. Rock NMR spectrum shows bimodal or triple-modal characteristics. Figure 1 is a sketch map of NMR spectrum of secondary porous core before and after centrifugation. The blue curve is the porosity component of core before centrifugation, the red curve is the porosity component of core after centrifugation, and the green curve and the brown curve are the accumulative porosity values calculated according to the porosity components of core before centrifugation and after centrifugation, respectively. It can be seen from the figure that the spectrum of core before centrifugation shows three peaks, while the spectrum after centrifugation shows two peaks. Its shape is close to the corresponding porosity component of different relaxation times. The third peak with the longest relaxation time before centrifugation disappears after centrifugation. The reason is that the fluid in the secondary large pore is discharged by centrifugation. The relaxation time corresponding to the black vertical line in the figure is the cutoff value of relaxation time to distinguish primary pore from secondary pore. The latter is referred to as the NMR cutoff value of secondary pore. The difference between the cumulative total porosity before centrifugation and the cumulative porosity before centrifugation (points A and B in Figure 1, points A and B correspond to the cumulative porosity values 1 and 2) is the secondary porosity of the core.
Figure 2 shows the NMR spectrum of core in well DF51 at 1402.15 m before and after centrifugation. Similar to the NMR spectrum in Figure 1, the porosity component of core before centrifugation presents obvious three-peak characteristics, and the porosity component after centrifugation presents obvious double-peak characteristics. Comparing the porosity component curves of the core before and after centrifugation, the relaxation time corresponding to the disappearance of the third peak with the longest relaxation time before centrifugation is 115000 μm, which is the cutoff value of secondary pore at this depth point. The cumulative total porosity before centrifugation and the cumulative total porosity before centrifugation corresponding to the secondary pore NMR cutoff (points A and B in Figure 2 correspond to and ) are 15.33% and 11.66%, respectively. Therefore, the secondary porosity of core is judged to be 3.67%.
In fact, the above method of judging the cutoff value of secondary pore is influenced by experimental factors such as centrifugal force of NMR experiment and artificial factors such as reading value, which cannot guarantee its accuracy. In addition, the distribution of NMR spectrum in some cores is not as typical as that in Figure 1 such as the distribution of NMR spectrum of core in well DF51 at 1428.52 m and 1442.54 m (as shown in Figures 3 and 4), which makes it difficult to distinguish secondary pore from primary pore clearly, and it is even more difficult to accurately determine the cutoff value of the secondary pore.
Therefore, it is a new problem how to obtain the general rule of the cutoff value of the secondary pore in reservoir with the existing data. In this paper, a new idea is proposed to use surface porosity and secondary surface porosity calculated by casting thin section to calibrate and calculate the cutoff value of the secondary pore and then use this cutoff value of the secondary pore calculate the secondary porosity of reef limestone reservoir.
2.2. Calculating Secondary Porosity by the Cutoff Value of Secondary Pore Calibrated by Casting Thin Section
2.2.1. The Relation between Surface Porosity, Secondary Surface Porosity, and Total Porosity
Using a polarizing microscope to observe the casting thin section of reef limestone reservoir, different pore types (intergranular pore, intragranular pore, organic visceral pore, intergranular dissolved pore, intragranular dissolved pore, cemented dissolved pore, matrix dissolved pore, karst cave, structural fracture, dissolution fracture, pressure solution fracture, etc.) can be identified and the relative content information of various pore types can be counted. Figure 5 is the photograph of some casting thin section in well DF51. From casting thin-section photographs, it can be seen that reservoir pore development and distribution are uneven, mainly secondary pore such as intergranular dissolved pore, intragranular dissolved pore, structural dissolved fracture, and karst cave.
Table 1 is the identification report of 1402.15 m casting thin section in well DF51. The type and relative content of the secondary pore counted by observing the casting thin section are recorded in the report. The identification report shows that the ratio of pore’s area and ratio of fissure’s area of the casting thin section are 7.5% and 1.5%, respectively. The total surface porosity is 9%. The secondary pore types include intragranular dissolved pore, intragranular dissolved pore, matrix dissolved pore, and structural dissolution fissure, and their relative content is 3%, 1%, 0.5%, and 1.5%, respectively, that is, the secondary surface porosity is 6%.
Secondary surface porosity and secondary porosity are two-dimensional and three-dimensional physical quantities, respectively. Therefore, in order to use secondary surface porosity to calibrate the cutoff value of the secondary pore to calculate reservoir secondary porosity, it is necessary to ensure that secondary surface porosity has a good correlation with secondary porosity. Table 2 shows the surface porosity, secondary surface porosity, and total porosity by logging interpretation of 25 casting thin sections at different depths in reef limestone reservoir section of well DF51 (the 25 samples have the similar mineral composition, each sample is a plunger sample, and NMR experiments are completed). Based on this, the corresponding surface porosity, secondary surface porosity, and total porosity by logging interpretation break-line chart are made (Figure 6). From Figure 6, it can be seen that surface porosity, secondary surface porosity, and total porosity by logging interpretation of the casting thin section at the same depth point are obviously different in numerical value, but the variation trend of porosity along the depth is the same, which shows that the three porosities have obvious correlation. Further analysis of the correlation between surface porosity and total porosity by logging interpretation shows that the linear correlation between them is more than 0.85 (Figure 7). This provides a theoretical basis for the idea that using surface porosity and secondary surface porosity calculated by casting thin section calibrate and calculate the cutoff value of secondary pore, and then, use this cutoff value of secondary pore to calculate the secondary porosity of reef limestone reservoir.
2.2.2. Calculating the Cutoff Value of Secondary Pore by Using Monte Carlo Method
Simulated annealing method models a physical process in which the atoms form a crystal (Metropolis et al., 1953). Simulated annealing process gradually decreases the temperature from “melting” to “freezing” until there are no further changes for the system. The physical process consists of thermal equilibriums at a set of temperatures. The simulation process at each temperature must have enough time (configuration) for the system to reach a steady state. Thus, we need a set of configurations to reach the thermal equilibrium at temperature . The configurations begin with an initial guess and computing corresponding energy of the system . We then generate a small random perturbation and compute the energy change for the system. The displacement is accepted if the energy change is smaller than 0 or conditionally accepted based on a user-defined probability function if the energy change is greater than 0. Thus, the main parameters in the annealing process include the starting temperature, the sequence of temperatures, the number of the simulation steps at each temperature, and the stopping criterion. We repeat the perturbation and acceptance steps many times until the system reaches thermal equilibrium at temperature .
Our paper replaces the energy function in simulated annealing method using the cost function shown in
Then, we maximize the energy function with with restricted to , a subset of , where is the number of the unknown values. Initially, a starting quality factor model, , is chosen arbitrarily on , and the initial value of the temperature, , is determined as described below. We denote at the iteration number of simulation at temperature . Then, the optimization process essentially comprises the following steps. (1)Employ Equation (4) to compute the relaxation time at each depth point using the initial quality factor model Then, use Equation (5) to compute the similarity values for relaxation time. At the same time set , , , and . (2)The quality factor model is perturbed randomly, and each value may remain unchanged changed with equal probability. We may decrease or increase the changed value by a predefined step, , with equal probability. We then obtain a new quality factor model, , by applying a three-dimension smoothing the altered quality factor model. We calculate if , increase by 1. This step is named as simulation step(3)We accept the simulation step unconditionally if , set and increase by 1(4)We accept the simulation step with probability, if . We first generate a uniformly distributed number between 0 and 1. We accept the simulation step if , set and increase by 1. This step aims to avoid the simulation step to escape out of local minimization maxima in its search for the global maxima(5)Go to step 2 if and(6)Increase by 1 if ; otherwise, set (7)Update the temperature, if , and increase by 1. Then, go to step 2(8)The simulation procedure stops if there is no accepted system configuration, . Then, the best solution is
Because of the good correlation between the surface porosity of casting thin section and the total porosity by logging interpretation, it can be considered that the secondary surface ratio and secondary porosity have the same good correlation. The Monte Carlo method is used to search for relaxation time and core NMR experimental data based on corresponding depth point B calculates the secondary porosity of NMR point by point according to the set step size and finds the relaxation time corresponding to the best correlation with the secondary porosity of thin films. The relaxation time is the cutoff value of the secondary pore.
Figure 8 shows the distribution of nuclear magnetic resonance spectra of multiple cores at different depths in well DF51, combined with other depths of cores in Figures 2–4. According to the analysis of the determination method of the secondary pore cutoff value described above, it is considered that the secondary pore cutoff value in the target formation of the well is distributed in the range of relaxation time of from 10000 to 1 000 000 000 μm, and its corresponding is between 4 and 6.
The range of value from 4 to 6 is set as the search interval, in which the search step is 0.5. The secondary pore NMR cutoff value is determined by the Monte Carlo method. In this process, the secondary porosity of 16 cores in Table 2 and the secondary porosity of thin section are used to participate in Monte Carlo simulation calculation (Table 3). The remaining 9 cores are used for method verification and accuracy analysis.
Through Monte Carlo simulation calculation, the secondary porosity of NMR is calculated point by point in the set relaxation time search interval of ; based on the NMR data, the Monte Carlo method is used in the search interval of from 4 to 6 according to the set step size and the experimental data of 16 cores, according to the set step size of 0.5. As shown in Figure 9, the calculation results show that the correlation between secondary porosity and secondary porosity of NMR experiment is the highest when equals 5.1. The corresponding secondary porosity is 125892.54 μm, and the calculation equation of secondary porosity of NMR is as follows:
3. Field Example
Equation (6) for calculating the secondary porosity of NMR is used when the cutoff value of secondary porosity is 125892.54 μm. According to the secondary porosity data of thin section, the secondary porosity of other 9 cores in DF51 well is calculated. The secondary porosity of 9 cores in is the true value, and the error analysis of the secondary porosity calculation results is made. Table 4 is the statistical error table of secondary porosity calculated by NMR. The average absolute error between the calculated secondary porosity of 9 cores and the true secondary porosity () is 0.391, and the average absolute error is 3.13, which indicates that the calculated secondary porosity has high accuracy.
The accuracy of secondary porosity calculation can be further verified by using production performance data of well DF51. Table 5 shows the productivity statistics of each sand layer in the production horizon (X3 horizon) of well DF51.
According to bubble diagram (bubble size corresponds to fluid production per meter) Figure 10 made with NMR surface porosity and fluid production index per meter and Figure 11 made with surface porosity by porosity spectrum and fluid production index per meter, it can be seen that the relationship between the NMR secondary porosity and the fluid production index per meter of the production layer X3 is better. It shows that the accuracy of secondary porosity calculation by NMR is higher than that by the popular method of secondary porosity calculation by porosity spectrum.
4. Conclusion
In this paper, we use the secondary porosity from casting thin section to scale the core NMR experimental data at different depths and find and determine the secondary pore NMR cutoff value by the Monte Carlo optimization algorithm. The method avoids the drawbacks of artificial identification method and greatly improves the accuracy of the secondary pore NMR cutoff value. And then, we analyze the function relationship between cutoff value and secondary porosity calculated by NMR experiment and calculate secondary porosity by this function relationship. Through the analysis of the correlation between the calculation results and reservoir productivity, the calculation accuracy of the secondary porosity method proposed in this paper is higher than that of the current mainstream porosity spectrum secondary porosity calculation method.
Data Availability
The analysis of test data used to support the findings of this study have not been made available because of the need for commercial confidentiality.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
We appreciate the financial support provided for the research in this paper from the National Science Foundation of China “Physical Parameters Modeling on Carbonate Fracture-Cavity Reservoirs” (Grant No. 41402113), the National Science and Technology Major Project of China “Logging Identification and Comprehensive Evaluation Technology of Low Permeability Tight Reservoir” (Grant No. 2016ZX05027-002-002), the Scientific Research Project of Hubei Provincial Department of Education “Property Modeling on Carbonate Cavity Reservoirs” (Grant No. B2013283), and the Productive Scientific Research Project of Zhanjiang Branch of CNOOC (China) Co., Ltd. “Logging comprehensive evaluation technology for fracture pore reservoir in Yongle 8 area” (Grant No. CCL2020ZJFN0278).