Method of Calculating Secondary Porosity of Reef Limestone Reservoir by Casting Thin Section Calibrating Nuclear Magnetic T2 Spectrum

Key Laboratory of Exploration Technologies for Oil and Gas Resources, Ministry of Education, Yangtze University, Wuhan Hubei, China College of Geophysics and Oil Resources, Yangtze University, Wuhan Hubei, China Research Institute of Petroleum Exploration and Development, Jilin Oil & Gas Field Branch, PetroChina, Songyuan, China Zhanjiang Branch, CNOOC Energy Technology & Services Co., Zhanjiang, China


Introduction
Reef limestone reservoir is one of the most important types of carbonate reservoirs [1][2][3][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][7][8][9][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][12][13][14][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][17][18][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: cut-off value of porosity spectrum = a × δ + median porosity, where δ means the porosity variance between the minimum porosity and the median porosity and a 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 a, 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 T 2 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 T 2 spectrum, the content of different T 2 components can be regarded as the number of different pore size components. Therefore, the multicomponent pore can be determined by different T 2 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-T 2 value of NMR is found by Monte Carlo method. Secondary porosity cutoff value of NMR T 2 is determined by this relationship. Finally, reservoir secondary porosity is calculated according to the cutoff value of T 2 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 T 2 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. 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 highprecision 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.

Geofluids
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 T 2 of NMR is proportional to the pore throat radius r c .
where C means the conversion coefficient between T 2 and r c . As long as the C value is determined, the pore radius distribution curve of rock can be obtained by using nuclear magnetic resonance T 2 spectrum, which is the theoretical basis of quantitative evaluation of rock pore structure by nuclear magnetic resonance technology. The magnitude of C value depends on the lateral surface relaxation rate ρ 2 and the pore shape factor F S . Because these two parameters of rock cannot be measured directly in practical research, the conversion coefficient C value can be obtained by calibrating the NMR T 2 spectrum with more accurate mercury injection data based on the correlation between the NMR T 2 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 T 2 spectrum and pore throat radius distribution in reef limestone reservoirs has not been solved well.

The Principle of Calculating Secondary Porosity of Reef
Limestone Reservoir by NMR T 2 Spectrum. Different relaxation time in NMR T 2 spectrum corresponds to different pore radius sizes. Based on the assumption of normal distribution model, for homogeneous reservoirs, rock NMR T 2 spectrum shows unimodal characteristics. For highly heterogeneous reservoirs, such as reef limestone reservoirs, primary and secondary pores are developed simultaneously. Rock NMR T 2 spectrum shows bimodal or triple-modal characteristics. Figure 1 is a sketch map of NMR T 2 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 T 2 spectrum of core before centrifugation shows three peaks, while the T 2 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 T 2 relaxation time to distinguish primary pore from secondary pore. The latter is referred to as the NMR T 2 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 T 2 spectrum of core in well DF51 at 1402.15 m before and after centrifugation. Similar to the NMR T 2 spectrum in Figure 1, the porosity component 3 Geofluids 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 T 2 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 T 2 cutoff (points A and B in Figure 2 correspond to φ1 and φ2) 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 T 2 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 T 2 spectrum in some cores is not as typical as that in Figure 1 such as the distribution of NMR T 2 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 T 2 cutoff value of the secondary pore.
Therefore, it is a new problem how to obtain the general rule of the T 2 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 T 2 cutoff value of the secondary pore and then use this T 2 cutoff value of the secondary pore calculate the secondary porosity of reef limestone reservoir.

Calculating Secondary Porosity by the T 2 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 T 2 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 sam-ples 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

Geofluids
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 T 2 cutoff value of secondary pore, and then, use this T 2 cutoff value of secondary pore to calculate the secondary porosity of reef limestone reservoir.

Calculating the T 2 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 T. The configurations begin with an initial guess and computing corresponding energy of the system E. 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 T.
Our paper replaces the energy function in simulated annealing method using the cost function shown in Then, we maximize the energy function EðQÞ with Q with restricted to Ω, a subset of R M , where M is the number    Initially, a starting quality factor model, Q 0 , is chosen arbitrarily on Ω, and the initial value of the temperature, T 0 , is determined as described below. We denote A k at the iteration number of simulation at temperature T k . Then, the Q optimization process essentially comprises the following steps.
(1) Employ Equation (4) to compute the T 2 relaxation time at each depth point using the initial quality fac- Then, use Equation (5) to compute the similarity values for T 2 relaxation time. At the same time set k = 0, u = 0, L k = 0, and A k = 0.
(2) The quality factor model is perturbed randomly, and each Q 0,i value may remain unchanged changed with equal probability. We may decrease or increase the changed Q 0,i value by a predefined step, ΔQ, with equal probability. We then obtain a new quality factor model, Q 1 , by applying a three-dimension smoothing the altered quality factor model. We calculate ΔE = EðQ 1 Þ − EðQ 0 Þ if Q 1 ∈ Ω, increase L k by 1. This step is named as simulation step (3) We accept the simulation step unconditionally if ΔE ≫ 0, set Q 0 = Q 1 and increase A k by 1 (4) We accept the simulation step with probability, P = e ΔE/T k if ΔE < 0. We first generate a uniformly distributed number between 0 and 1. We accept the simulation step if ρ < P, set Q 0 = Q 1 and increase A k 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 A k < A min and L k < L min (6) Increase U by 1 if A k = 0; otherwise, set U = 0 (7) Update the temperature, T k+1 < 0:99T k if U < U max , and increase k by 1. Then, go to step 2 (8) The simulation procedure stops if there is no accepted system configuration, U < U max . Then, the best solution is Q 0 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 T 2 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 T 2 relaxation time corresponding to the best correlation with the secondary porosity of thin films. The relaxation time is the T 2 cutoff value of the secondary pore. Figure 8 shows the distribution of nuclear magnetic resonance T 2 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 T 2 cutoff value described above, it is considered that the secondary pore T 2 cutoff value in the target formation of the well is distributed in the range of relaxation time of T 2 from 10000 to 1 000 000 000 μm, and its corresponding log 10 T 2 is between 4 and 6.
The range of log 10 T 2 value from 4 to 6 is set as the search interval, in which the search step is 0.5. The secondary pore NMR T 2 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 8 Geofluids (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 T 2 ; based on the NMR data, the Monte Carlo method is used in the search interval of log 10 T 2 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

Field Example
Equation (6) for calculating the secondary porosity of NMR is used when the T 2 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 log 10 T 2 = 5:1 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 (log 10 T 2 = 5:1) 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.    Surface porosity by porosity spectrum (%) Fluid production index per meter (bbl/psi·m·day) y = 0.0106x + 0.0111 R 2 = 0.902 Figure 11: Intersection diagram of surface porosity by porosity spectrum and fluid production index per meter in production layer X3.

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 T 2 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 T 2 cutoff value. And then, we analyze the function relationship between T 2 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.