# Method of Calculating Secondary Porosity of Reef Limestone Reservoir by Casting Thin Section Calibrating Nuclear MagneticT2 Spectrum

Authors: Bin Zhao; Zhaoping Li; Chuqiao Gao; Guoyi Zhang; Jinbo Wu; Gang Gao; Yang Tang
Journal: Geofluids
 (2021)

---

## 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 resonanceT2 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 T2 relaxation time corresponding to the best correlation between the secondary porosity and the secondary surface porosity of casting thin section, that is, the T2 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.


---

## Body

## 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:
(1)cut‐offvalueofporosityspectrum=a×δ+medianporosity,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 NMRT2 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 T2 spectrum, the content of different T2 components can be regarded as the number of different pore size components. Therefore, the multicomponent pore can be determined by different T2 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-T2 value of NMR is found by Monte Carlo method. Secondary porosity cutoff value of NMR T2 is determined by this relationship. Finally, reservoir secondary porosity is calculated according to the cutoff value of T2 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 T2 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 NMRT2 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 timeT2 of NMR is proportional to the pore throat radius rc.
(2)T2=C×rc,where C means the conversion coefficient between T2 and rc.As long as theC value is determined, the pore radius distribution curve of rock can be obtained by using nuclear magnetic resonance T2 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 FS. 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 T2 spectrum with more accurate mercury injection data based on the correlation between the NMR T2 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 T2 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 NMRT2 Spectrum

Different relaxation time in NMRT2 spectrum corresponds to different pore radius sizes. Based on the assumption of normal distribution model, for homogeneous reservoirs, rock NMR T2 spectrum shows unimodal characteristics. For highly heterogeneous reservoirs, such as reef limestone reservoirs, primary and secondary pores are developed simultaneously. Rock NMR T2 spectrum shows bimodal or triple-modal characteristics. Figure 1 is a sketch map of NMR T2 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 T2 spectrum of core before centrifugation shows three peaks, while the T2 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 T2 relaxation time to distinguish primary pore from secondary pore. The latter is referred to as the NMR T2 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 1
NMRT2 spectrum of reef limestone reservoir core before and after centrifugation.Figure2 shows the NMR T2 spectrum of core in well DF51 at 1402.15 m before and after centrifugation. Similar to the NMR T2 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 T2 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 T2 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%.Figure 2
NMRT2 spectrum of core in well DF51 at 1402.15 m.In fact, the above method of judging theT2 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 T2 spectrum in some cores is not as typical as that in Figure 1 such as the distribution of NMR T2 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 T2 cutoff value of the secondary pore.Figure 3
NMRT2 spectrum of core in well DF51 at 1428.52 m.Figure 4
NMRT2 spectrum of core in well DF51 at 1442.54 m.Therefore, it is a new problem how to obtain the general rule of theT2 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 T2 cutoff value of the secondary pore and then use this T2 cutoff value of the secondary pore calculate the secondary porosity of reef limestone reservoir.


### 2.2. Calculating Secondary Porosity by theT2 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. Figure5 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.Figure 5
The casting thin-section photograph of the corresponding depth of NMR core in well DF51.Table1 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%.Table 1
Identification report of casting thin section (well DF51 1402.15 m).
Well: DF51Depth: 1402.15 mType of specimen: coreReservoir space typePoreRatio of pore’s area (%)FissureRatio of fissure’s area (%)Primary poreSecondary poreTypeGenesisQuantityWidthFilling compositionFilling conditionIntergranular poreIntragranular poreOrganic visceral poreIntergranular dissolved poreIntragranular dissolved poreCemented dissolved poreMatrix dissolved poreKarst cave7.5Open fractureStructural dissolutionSeveral0.01-0.1 mmNoNo1.5//331/0.5///////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 theT2 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 T2 cutoff value of secondary pore, and then, use this T2 cutoff value of secondary pore to calculate the secondary porosity of reef limestone reservoir.Table 2
Statistical tables of surface porosity, secondary surface porosity, and total porosity by logging interpretation of 25 casting thin sections in well DF51.
No.Depth (m)Surface porosity (%)Secondary surface porosity (%)Total porosity by logging interpretation (%)11389.88.84.513.4221389.9312.35.514.1431391.27.5617.1541399.5211230.0851400.3201028.6761400.5517726.8971402.159611.5681404.5585.510.491409.49161020.38101410.5111.59.515.15111412.66161121.55121421.07188.521.46131422.0115.5522.27141428.5210217.86151430.84154.524.45161434.3912120.46171435.0521128.66181437.39251632.29191441.2915520.21201442.546.5514.11211446.94181224.94221449.4517823.46231453.2414618.48241459.99241430.23251465.92121013.71Figure 6
Trend charts of surface porosity, secondary surface porosity, and total porosity by logging interpretation of casting thin section with depth.Figure 7
Relation between surface porosity and total porosity by logging interpretation.


#### 2.2.2. Calculating theT2 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 temperatureT. 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(3)Q=argmaxEQQ,EQ=∑m=1MEmQ,EmQ=∑f=fminfmax∑h=2HB¯f,tm,xh+Bf,t1,x1H∑f=fminfmax∑h=2HB¯2f,tm,xh+B2f,t1,x1.Then, we maximize the energy functionEQ with Q with restricted to Ω, a subset of RM, where M is the number of the unknown Q values. Initially, a starting quality factor model, Q0, is chosen arbitrarily on Ω, and the initial value of the temperature, T0, is determined as described below. We denote Ak at the iteration number of simulation at temperature Tk. Then, the Q optimization process essentially comprises the following steps.
(1)
Employ Equation (4) to compute the T2 relaxation time at each depth point using the initial quality factor model Q0(4)B¯f,tm,xh=Bf,tm,x1⋅e∑k=1mπfΔtxh,k/qk.Then, use Equation (5) to compute the similarity values for T2 relaxation time. At the same time set k=0, u=0, Lk=0, and Ak=0.
(5)t2=t02+x2vrms2t0.(2)
The quality factor model is perturbed randomly, and eachQ0,i value may remain unchanged changed with equal probability. We may decrease or increase the changed Q0,i value by a predefined step, ∆Q, with equal probability. We then obtain a new quality factor model, Q1, by applying a three-dimension smoothing the altered quality factor model. We calculate ∆E=EQ1−EQ0 if Q1∈Ω, increase Lk by 1. This step is named as simulation step(3)
We accept the simulation step unconditionally if∆E≫0, set Q0=Q1 and increase Ak by 1(4)
We accept the simulation step with probability,P=e∆E/Tk if ∆E<0. We first generate a uniformly distributed number between 0 and 1. We accept the simulation step if ρ<P, set Q0=Q1 and increase Ak 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 ifAk<Amin andLk<Lmin(6)
IncreaseU by 1 if Ak=0; otherwise, set U=0(7)
Update the temperature,Tk+1<0.99Tk if U<Umax, and increase k by 1. Then, go to step 2(8)
The simulation procedure stops if there is no accepted system configuration,U<Umax. Then, the best solution is Q0Because 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 forT2 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 T2 relaxation time corresponding to the best correlation with the secondary porosity of thin films. The relaxation time is the T2 cutoff value of the secondary pore.Figure8 shows the distribution of nuclear magnetic resonance T2 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 T2 cutoff value described above, it is considered that the secondary pore T2 cutoff value in the target formation of the well is distributed in the range of relaxation time of T2 from 10000 to 1 000 000 000 μm, and its corresponding log10T2 is between 4 and 6.Figure 8
NMRT2 spectrum of cores in well DF51 at different depths.The range oflog10T2 value from 4 to 6 is set as the search interval, in which the search step is 0.5. The secondary pore NMR T2 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.Table 3
Statistical table of secondary porosity and slice secondary face porosity in 16 cores of well 3DF51.
No.12456710121415161720232425Depth (m)1389.81389.931399.51400.31400.551402.151410.511421.071428.521430.841434.391435.051442.541453.241459.991465.92SSP (%)4.55.51210769.58.524.511561410log10T2=410.3110.9526.8525.9923.019.2422.4820.414.7012.525.676.0815.0015.5723.6113.69log10T2=4.059.9810.6826.8125.8723.018.9622.4220.384.6712.305.666.0115.0015.5723.4913.64log10T2=4.19.6310.3926.7625.7223.018.6722.3520.324.6412.075.645.9415.0015.5623.3413.62log10T2=4.159.2610.0726.6725.5423.018.3522.2420.204.6011.815.625.8515.0015.5223.1613.62log10T2=4.28.899.7426.5525.3323.018.0322.1020.024.5511.565.595.7415.0015.4322.9713.62log10T2=4.258.519.3826.3925.0723.017.7021.9219.784.5011.285.545.6215.0015.3022.7913.62log10T2=4.38.129.0026.1824.7822.987.3621.6919.474.4411.005.495.4715.0015.1122.6313.62log10T2=4.357.738.6125.9324.4422.887.0121.4119.104.3810.715.415.3115.0014.8822.5013.62log10T2=4.47.338.2125.6324.0622.696.6721.0618.654.3110.415.325.1414.9714.5922.4213.62log10T2=4.456.917.7725.2423.5822.376.2920.6218.084.2410.095.214.9314.8914.2322.4213.62log10T2=4.56.577.4024.8623.1422.015.9820.1717.544.179.835.094.7414.7813.8722.4213.62log10T2=4.556.216.9924.4022.6121.535.6619.6216.904.109.554.954.5314.6213.4622.4213.62log10T2=4.65.866.5923.8922.0420.995.3519.0216.224.049.294.804.3114.4213.0322.4213.62log10T2=4.655.526.2023.3521.4320.385.0618.3715.513.999.054.634.0914.1612.5922.4213.62log10T2=4.75.225.8522.8120.8419.784.8017.7314.833.968.844.473.8913.8812.1822.4213.62log10T2=4.754.895.4722.1720.1419.064.5516.9814.063.948.644.283.6713.5111.7422.4213.62log10T2=4.84.615.1421.5619.4618.394.3416.2813.373.938.494.103.4813.1111.3822.4213.62log10T2=4.854.344.8220.9618.7817.734.1615.6012.733.938.363.933.3112.6711.0722.4213.62log10T2=4.94.094.5320.3618.1117.114.0014.9512.143.938.263.783.1612.1910.8222.4213.62log10T2=4.953.854.2519.7817.4416.543.8714.3611.613.938.183.643.0411.6910.6522.4213.62log10T2=53.624.0019.2116.7916.023.7713.8411.173.938.123.522.9411.1710.5622.4213.62log10T2=5.053.413.7618.6716.1615.583.6913.4110.823.938.083.432.8710.6310.5622.4213.62log10T2=5.13.223.5418.1315.5515.203.6313.0710.543.938.033.362.8310.0810.5622.4213.62log10T2=5.153.033.3317.6014.9314.853.5712.8010.323.937.973.312.809.5210.5622.4213.62log10T2=5.22.853.1317.0414.3214.533.5112.5810.143.937.883.272.798.9510.5622.4213.62log10T2=5.252.682.9316.4513.6914.203.4412.409.993.937.743.242.788.3810.5622.4213.62log10T2=5.32.512.7115.7612.9813.823.3512.219.823.937.533.222.787.7710.5622.4213.62log10T2=5.352.362.5215.0312.2813.383.2412.009.623.937.253.182.757.2010.5622.4213.62log10T2=5.42.212.3014.1411.4712.823.0911.699.333.896.863.122.706.5910.4422.4213.62log10T2=5.452.062.0613.1210.5712.112.8911.278.933.796.383.032.615.9510.1822.4213.62log10T2=5.51.911.8211.969.5811.252.6610.718.423.605.812.892.485.299.7422.3113.62log10T2=5.551.781.5910.808.6110.342.4110.087.853.385.222.742.334.689.2021.9813.62log10T2=5.61.631.309.247.339.052.069.096.993.014.412.482.093.928.2920.5713.04log10T2=5.651.491.047.726.117.721.728.016.072.603.622.201.833.217.2918.5611.96log10T2=5.71.360.796.174.876.331.366.835.072.152.831.891.542.536.1716.0510.49log10T2=5.751.220.554.653.664.921.015.574.041.692.061.551.241.884.9813.188.73log10T2=5.81.090.353.222.533.550.684.293.001.231.351.210.931.283.7710.146.80log10T2=5.850.950.181.981.552.320.403.072.050.800.750.880.640.772.637.174.88log10T2=5.90.820.060.960.751.260.181.931.190.440.300.560.380.361.584.363.02log10T2=5.950.690.010.310.240.520.041.030.550.180.070.310.180.110.782.171.54log10T2=60.560.000.040.030.130.000.390.160.040.000.130.050.010.240.670.50Through Monte Carlo simulation calculation, the secondary porosity of NMR is calculated point by point in the set relaxation time search interval ofT2; based on the NMR data, the Monte Carlo method is used in the search interval of log10T2 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 log10T2 equals 5.1. The corresponding secondary porosity is 125892.54 μm, and the calculation equation of secondary porosity of NMR is as follows:
(6)ϕNMRsecondaryporosity=1.4477×ϕsecondarysurfaceporosity+0.2196.Figure 9
Linear relationship between secondary porosity and secondary porosity (16 cores) in NMR experiments atlog10T2=5.1.


## 2.1. Calculating Secondary Porosity of Reef Limestone Reservoir by NMRT2 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 timeT2 of NMR is proportional to the pore throat radius rc.
(2)T2=C×rc,where C means the conversion coefficient between T2 and rc.As long as theC value is determined, the pore radius distribution curve of rock can be obtained by using nuclear magnetic resonance T2 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 FS. 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 T2 spectrum with more accurate mercury injection data based on the correlation between the NMR T2 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 T2 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 NMRT2 Spectrum

Different relaxation time in NMRT2 spectrum corresponds to different pore radius sizes. Based on the assumption of normal distribution model, for homogeneous reservoirs, rock NMR T2 spectrum shows unimodal characteristics. For highly heterogeneous reservoirs, such as reef limestone reservoirs, primary and secondary pores are developed simultaneously. Rock NMR T2 spectrum shows bimodal or triple-modal characteristics. Figure 1 is a sketch map of NMR T2 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 T2 spectrum of core before centrifugation shows three peaks, while the T2 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 T2 relaxation time to distinguish primary pore from secondary pore. The latter is referred to as the NMR T2 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 1
NMRT2 spectrum of reef limestone reservoir core before and after centrifugation.Figure2 shows the NMR T2 spectrum of core in well DF51 at 1402.15 m before and after centrifugation. Similar to the NMR T2 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 T2 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 T2 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%.Figure 2
NMRT2 spectrum of core in well DF51 at 1402.15 m.In fact, the above method of judging theT2 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 T2 spectrum in some cores is not as typical as that in Figure 1 such as the distribution of NMR T2 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 T2 cutoff value of the secondary pore.Figure 3
NMRT2 spectrum of core in well DF51 at 1428.52 m.Figure 4
NMRT2 spectrum of core in well DF51 at 1442.54 m.Therefore, it is a new problem how to obtain the general rule of theT2 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 T2 cutoff value of the secondary pore and then use this T2 cutoff value of the secondary pore calculate the secondary porosity of reef limestone reservoir.


## 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 timeT2 of NMR is proportional to the pore throat radius rc.
(2)T2=C×rc,where C means the conversion coefficient between T2 and rc.As long as theC value is determined, the pore radius distribution curve of rock can be obtained by using nuclear magnetic resonance T2 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 FS. 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 T2 spectrum with more accurate mercury injection data based on the correlation between the NMR T2 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 T2 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 NMRT2 Spectrum

Different relaxation time in NMRT2 spectrum corresponds to different pore radius sizes. Based on the assumption of normal distribution model, for homogeneous reservoirs, rock NMR T2 spectrum shows unimodal characteristics. For highly heterogeneous reservoirs, such as reef limestone reservoirs, primary and secondary pores are developed simultaneously. Rock NMR T2 spectrum shows bimodal or triple-modal characteristics. Figure 1 is a sketch map of NMR T2 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 T2 spectrum of core before centrifugation shows three peaks, while the T2 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 T2 relaxation time to distinguish primary pore from secondary pore. The latter is referred to as the NMR T2 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 1
NMRT2 spectrum of reef limestone reservoir core before and after centrifugation.Figure2 shows the NMR T2 spectrum of core in well DF51 at 1402.15 m before and after centrifugation. Similar to the NMR T2 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 T2 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 T2 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%.Figure 2
NMRT2 spectrum of core in well DF51 at 1402.15 m.In fact, the above method of judging theT2 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 T2 spectrum in some cores is not as typical as that in Figure 1 such as the distribution of NMR T2 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 T2 cutoff value of the secondary pore.Figure 3
NMRT2 spectrum of core in well DF51 at 1428.52 m.Figure 4
NMRT2 spectrum of core in well DF51 at 1442.54 m.Therefore, it is a new problem how to obtain the general rule of theT2 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 T2 cutoff value of the secondary pore and then use this T2 cutoff value of the secondary pore calculate the secondary porosity of reef limestone reservoir.


## 2.2. Calculating Secondary Porosity by theT2 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. Figure5 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.Figure 5
The casting thin-section photograph of the corresponding depth of NMR core in well DF51.Table1 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%.Table 1
Identification report of casting thin section (well DF51 1402.15 m).
Well: DF51Depth: 1402.15 mType of specimen: coreReservoir space typePoreRatio of pore’s area (%)FissureRatio of fissure’s area (%)Primary poreSecondary poreTypeGenesisQuantityWidthFilling compositionFilling conditionIntergranular poreIntragranular poreOrganic visceral poreIntergranular dissolved poreIntragranular dissolved poreCemented dissolved poreMatrix dissolved poreKarst cave7.5Open fractureStructural dissolutionSeveral0.01-0.1 mmNoNo1.5//331/0.5///////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 theT2 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 T2 cutoff value of secondary pore, and then, use this T2 cutoff value of secondary pore to calculate the secondary porosity of reef limestone reservoir.Table 2
Statistical tables of surface porosity, secondary surface porosity, and total porosity by logging interpretation of 25 casting thin sections in well DF51.
No.Depth (m)Surface porosity (%)Secondary surface porosity (%)Total porosity by logging interpretation (%)11389.88.84.513.4221389.9312.35.514.1431391.27.5617.1541399.5211230.0851400.3201028.6761400.5517726.8971402.159611.5681404.5585.510.491409.49161020.38101410.5111.59.515.15111412.66161121.55121421.07188.521.46131422.0115.5522.27141428.5210217.86151430.84154.524.45161434.3912120.46171435.0521128.66181437.39251632.29191441.2915520.21201442.546.5514.11211446.94181224.94221449.4517823.46231453.2414618.48241459.99241430.23251465.92121013.71Figure 6
Trend charts of surface porosity, secondary surface porosity, and total porosity by logging interpretation of casting thin section with depth.Figure 7
Relation between surface porosity and total porosity by logging interpretation.


### 2.2.2. Calculating theT2 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 temperatureT. 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(3)Q=argmaxEQQ,EQ=∑m=1MEmQ,EmQ=∑f=fminfmax∑h=2HB¯f,tm,xh+Bf,t1,x1H∑f=fminfmax∑h=2HB¯2f,tm,xh+B2f,t1,x1.Then, we maximize the energy functionEQ with Q with restricted to Ω, a subset of RM, where M is the number of the unknown Q values. Initially, a starting quality factor model, Q0, is chosen arbitrarily on Ω, and the initial value of the temperature, T0, is determined as described below. We denote Ak at the iteration number of simulation at temperature Tk. Then, the Q optimization process essentially comprises the following steps.
(1)
Employ Equation (4) to compute the T2 relaxation time at each depth point using the initial quality factor model Q0(4)B¯f,tm,xh=Bf,tm,x1⋅e∑k=1mπfΔtxh,k/qk.Then, use Equation (5) to compute the similarity values for T2 relaxation time. At the same time set k=0, u=0, Lk=0, and Ak=0.
(5)t2=t02+x2vrms2t0.(2)
The quality factor model is perturbed randomly, and eachQ0,i value may remain unchanged changed with equal probability. We may decrease or increase the changed Q0,i value by a predefined step, ∆Q, with equal probability. We then obtain a new quality factor model, Q1, by applying a three-dimension smoothing the altered quality factor model. We calculate ∆E=EQ1−EQ0 if Q1∈Ω, increase Lk by 1. This step is named as simulation step(3)
We accept the simulation step unconditionally if∆E≫0, set Q0=Q1 and increase Ak by 1(4)
We accept the simulation step with probability,P=e∆E/Tk if ∆E<0. We first generate a uniformly distributed number between 0 and 1. We accept the simulation step if ρ<P, set Q0=Q1 and increase Ak 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 ifAk<Amin andLk<Lmin(6)
IncreaseU by 1 if Ak=0; otherwise, set U=0(7)
Update the temperature,Tk+1<0.99Tk if U<Umax, and increase k by 1. Then, go to step 2(8)
The simulation procedure stops if there is no accepted system configuration,U<Umax. Then, the best solution is Q0Because 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 forT2 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 T2 relaxation time corresponding to the best correlation with the secondary porosity of thin films. The relaxation time is the T2 cutoff value of the secondary pore.Figure8 shows the distribution of nuclear magnetic resonance T2 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 T2 cutoff value described above, it is considered that the secondary pore T2 cutoff value in the target formation of the well is distributed in the range of relaxation time of T2 from 10000 to 1 000 000 000 μm, and its corresponding log10T2 is between 4 and 6.Figure 8
NMRT2 spectrum of cores in well DF51 at different depths.The range oflog10T2 value from 4 to 6 is set as the search interval, in which the search step is 0.5. The secondary pore NMR T2 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.Table 3
Statistical table of secondary porosity and slice secondary face porosity in 16 cores of well 3DF51.
No.12456710121415161720232425Depth (m)1389.81389.931399.51400.31400.551402.151410.511421.071428.521430.841434.391435.051442.541453.241459.991465.92SSP (%)4.55.51210769.58.524.511561410log10T2=410.3110.9526.8525.9923.019.2422.4820.414.7012.525.676.0815.0015.5723.6113.69log10T2=4.059.9810.6826.8125.8723.018.9622.4220.384.6712.305.666.0115.0015.5723.4913.64log10T2=4.19.6310.3926.7625.7223.018.6722.3520.324.6412.075.645.9415.0015.5623.3413.62log10T2=4.159.2610.0726.6725.5423.018.3522.2420.204.6011.815.625.8515.0015.5223.1613.62log10T2=4.28.899.7426.5525.3323.018.0322.1020.024.5511.565.595.7415.0015.4322.9713.62log10T2=4.258.519.3826.3925.0723.017.7021.9219.784.5011.285.545.6215.0015.3022.7913.62log10T2=4.38.129.0026.1824.7822.987.3621.6919.474.4411.005.495.4715.0015.1122.6313.62log10T2=4.357.738.6125.9324.4422.887.0121.4119.104.3810.715.415.3115.0014.8822.5013.62log10T2=4.47.338.2125.6324.0622.696.6721.0618.654.3110.415.325.1414.9714.5922.4213.62log10T2=4.456.917.7725.2423.5822.376.2920.6218.084.2410.095.214.9314.8914.2322.4213.62log10T2=4.56.577.4024.8623.1422.015.9820.1717.544.179.835.094.7414.7813.8722.4213.62log10T2=4.556.216.9924.4022.6121.535.6619.6216.904.109.554.954.5314.6213.4622.4213.62log10T2=4.65.866.5923.8922.0420.995.3519.0216.224.049.294.804.3114.4213.0322.4213.62log10T2=4.655.526.2023.3521.4320.385.0618.3715.513.999.054.634.0914.1612.5922.4213.62log10T2=4.75.225.8522.8120.8419.784.8017.7314.833.968.844.473.8913.8812.1822.4213.62log10T2=4.754.895.4722.1720.1419.064.5516.9814.063.948.644.283.6713.5111.7422.4213.62log10T2=4.84.615.1421.5619.4618.394.3416.2813.373.938.494.103.4813.1111.3822.4213.62log10T2=4.854.344.8220.9618.7817.734.1615.6012.733.938.363.933.3112.6711.0722.4213.62log10T2=4.94.094.5320.3618.1117.114.0014.9512.143.938.263.783.1612.1910.8222.4213.62log10T2=4.953.854.2519.7817.4416.543.8714.3611.613.938.183.643.0411.6910.6522.4213.62log10T2=53.624.0019.2116.7916.023.7713.8411.173.938.123.522.9411.1710.5622.4213.62log10T2=5.053.413.7618.6716.1615.583.6913.4110.823.938.083.432.8710.6310.5622.4213.62log10T2=5.13.223.5418.1315.5515.203.6313.0710.543.938.033.362.8310.0810.5622.4213.62log10T2=5.153.033.3317.6014.9314.853.5712.8010.323.937.973.312.809.5210.5622.4213.62log10T2=5.22.853.1317.0414.3214.533.5112.5810.143.937.883.272.798.9510.5622.4213.62log10T2=5.252.682.9316.4513.6914.203.4412.409.993.937.743.242.788.3810.5622.4213.62log10T2=5.32.512.7115.7612.9813.823.3512.219.823.937.533.222.787.7710.5622.4213.62log10T2=5.352.362.5215.0312.2813.383.2412.009.623.937.253.182.757.2010.5622.4213.62log10T2=5.42.212.3014.1411.4712.823.0911.699.333.896.863.122.706.5910.4422.4213.62log10T2=5.452.062.0613.1210.5712.112.8911.278.933.796.383.032.615.9510.1822.4213.62log10T2=5.51.911.8211.969.5811.252.6610.718.423.605.812.892.485.299.7422.3113.62log10T2=5.551.781.5910.808.6110.342.4110.087.853.385.222.742.334.689.2021.9813.62log10T2=5.61.631.309.247.339.052.069.096.993.014.412.482.093.928.2920.5713.04log10T2=5.651.491.047.726.117.721.728.016.072.603.622.201.833.217.2918.5611.96log10T2=5.71.360.796.174.876.331.366.835.072.152.831.891.542.536.1716.0510.49log10T2=5.751.220.554.653.664.921.015.574.041.692.061.551.241.884.9813.188.73log10T2=5.81.090.353.222.533.550.684.293.001.231.351.210.931.283.7710.146.80log10T2=5.850.950.181.981.552.320.403.072.050.800.750.880.640.772.637.174.88log10T2=5.90.820.060.960.751.260.181.931.190.440.300.560.380.361.584.363.02log10T2=5.950.690.010.310.240.520.041.030.550.180.070.310.180.110.782.171.54log10T2=60.560.000.040.030.130.000.390.160.040.000.130.050.010.240.670.50Through Monte Carlo simulation calculation, the secondary porosity of NMR is calculated point by point in the set relaxation time search interval ofT2; based on the NMR data, the Monte Carlo method is used in the search interval of log10T2 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 log10T2 equals 5.1. The corresponding secondary porosity is 125892.54 μm, and the calculation equation of secondary porosity of NMR is as follows:
(6)ϕNMRsecondaryporosity=1.4477×ϕsecondarysurfaceporosity+0.2196.Figure 9
Linear relationship between secondary porosity and secondary porosity (16 cores) in NMR experiments atlog10T2=5.1.


## 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. Figure5 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.Figure 5
The casting thin-section photograph of the corresponding depth of NMR core in well DF51.Table1 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%.Table 1
Identification report of casting thin section (well DF51 1402.15 m).
Well: DF51Depth: 1402.15 mType of specimen: coreReservoir space typePoreRatio of pore’s area (%)FissureRatio of fissure’s area (%)Primary poreSecondary poreTypeGenesisQuantityWidthFilling compositionFilling conditionIntergranular poreIntragranular poreOrganic visceral poreIntergranular dissolved poreIntragranular dissolved poreCemented dissolved poreMatrix dissolved poreKarst cave7.5Open fractureStructural dissolutionSeveral0.01-0.1 mmNoNo1.5//331/0.5///////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 theT2 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 T2 cutoff value of secondary pore, and then, use this T2 cutoff value of secondary pore to calculate the secondary porosity of reef limestone reservoir.Table 2
Statistical tables of surface porosity, secondary surface porosity, and total porosity by logging interpretation of 25 casting thin sections in well DF51.
No.Depth (m)Surface porosity (%)Secondary surface porosity (%)Total porosity by logging interpretation (%)11389.88.84.513.4221389.9312.35.514.1431391.27.5617.1541399.5211230.0851400.3201028.6761400.5517726.8971402.159611.5681404.5585.510.491409.49161020.38101410.5111.59.515.15111412.66161121.55121421.07188.521.46131422.0115.5522.27141428.5210217.86151430.84154.524.45161434.3912120.46171435.0521128.66181437.39251632.29191441.2915520.21201442.546.5514.11211446.94181224.94221449.4517823.46231453.2414618.48241459.99241430.23251465.92121013.71Figure 6
Trend charts of surface porosity, secondary surface porosity, and total porosity by logging interpretation of casting thin section with depth.Figure 7
Relation between surface porosity and total porosity by logging interpretation.


## 2.2.2. Calculating theT2 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 temperatureT. 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(3)Q=argmaxEQQ,EQ=∑m=1MEmQ,EmQ=∑f=fminfmax∑h=2HB¯f,tm,xh+Bf,t1,x1H∑f=fminfmax∑h=2HB¯2f,tm,xh+B2f,t1,x1.Then, we maximize the energy functionEQ with Q with restricted to Ω, a subset of RM, where M is the number of the unknown Q values. Initially, a starting quality factor model, Q0, is chosen arbitrarily on Ω, and the initial value of the temperature, T0, is determined as described below. We denote Ak at the iteration number of simulation at temperature Tk. Then, the Q optimization process essentially comprises the following steps.
(1)
Employ Equation (4) to compute the T2 relaxation time at each depth point using the initial quality factor model Q0(4)B¯f,tm,xh=Bf,tm,x1⋅e∑k=1mπfΔtxh,k/qk.Then, use Equation (5) to compute the similarity values for T2 relaxation time. At the same time set k=0, u=0, Lk=0, and Ak=0.
(5)t2=t02+x2vrms2t0.(2)
The quality factor model is perturbed randomly, and eachQ0,i value may remain unchanged changed with equal probability. We may decrease or increase the changed Q0,i value by a predefined step, ∆Q, with equal probability. We then obtain a new quality factor model, Q1, by applying a three-dimension smoothing the altered quality factor model. We calculate ∆E=EQ1−EQ0 if Q1∈Ω, increase Lk by 1. This step is named as simulation step(3)
We accept the simulation step unconditionally if∆E≫0, set Q0=Q1 and increase Ak by 1(4)
We accept the simulation step with probability,P=e∆E/Tk if ∆E<0. We first generate a uniformly distributed number between 0 and 1. We accept the simulation step if ρ<P, set Q0=Q1 and increase Ak 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 ifAk<Amin andLk<Lmin(6)
IncreaseU by 1 if Ak=0; otherwise, set U=0(7)
Update the temperature,Tk+1<0.99Tk if U<Umax, and increase k by 1. Then, go to step 2(8)
The simulation procedure stops if there is no accepted system configuration,U<Umax. Then, the best solution is Q0Because 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 forT2 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 T2 relaxation time corresponding to the best correlation with the secondary porosity of thin films. The relaxation time is the T2 cutoff value of the secondary pore.Figure8 shows the distribution of nuclear magnetic resonance T2 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 T2 cutoff value described above, it is considered that the secondary pore T2 cutoff value in the target formation of the well is distributed in the range of relaxation time of T2 from 10000 to 1 000 000 000 μm, and its corresponding log10T2 is between 4 and 6.Figure 8
NMRT2 spectrum of cores in well DF51 at different depths.The range oflog10T2 value from 4 to 6 is set as the search interval, in which the search step is 0.5. The secondary pore NMR T2 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.Table 3
Statistical table of secondary porosity and slice secondary face porosity in 16 cores of well 3DF51.
No.12456710121415161720232425Depth (m)1389.81389.931399.51400.31400.551402.151410.511421.071428.521430.841434.391435.051442.541453.241459.991465.92SSP (%)4.55.51210769.58.524.511561410log10T2=410.3110.9526.8525.9923.019.2422.4820.414.7012.525.676.0815.0015.5723.6113.69log10T2=4.059.9810.6826.8125.8723.018.9622.4220.384.6712.305.666.0115.0015.5723.4913.64log10T2=4.19.6310.3926.7625.7223.018.6722.3520.324.6412.075.645.9415.0015.5623.3413.62log10T2=4.159.2610.0726.6725.5423.018.3522.2420.204.6011.815.625.8515.0015.5223.1613.62log10T2=4.28.899.7426.5525.3323.018.0322.1020.024.5511.565.595.7415.0015.4322.9713.62log10T2=4.258.519.3826.3925.0723.017.7021.9219.784.5011.285.545.6215.0015.3022.7913.62log10T2=4.38.129.0026.1824.7822.987.3621.6919.474.4411.005.495.4715.0015.1122.6313.62log10T2=4.357.738.6125.9324.4422.887.0121.4119.104.3810.715.415.3115.0014.8822.5013.62log10T2=4.47.338.2125.6324.0622.696.6721.0618.654.3110.415.325.1414.9714.5922.4213.62log10T2=4.456.917.7725.2423.5822.376.2920.6218.084.2410.095.214.9314.8914.2322.4213.62log10T2=4.56.577.4024.8623.1422.015.9820.1717.544.179.835.094.7414.7813.8722.4213.62log10T2=4.556.216.9924.4022.6121.535.6619.6216.904.109.554.954.5314.6213.4622.4213.62log10T2=4.65.866.5923.8922.0420.995.3519.0216.224.049.294.804.3114.4213.0322.4213.62log10T2=4.655.526.2023.3521.4320.385.0618.3715.513.999.054.634.0914.1612.5922.4213.62log10T2=4.75.225.8522.8120.8419.784.8017.7314.833.968.844.473.8913.8812.1822.4213.62log10T2=4.754.895.4722.1720.1419.064.5516.9814.063.948.644.283.6713.5111.7422.4213.62log10T2=4.84.615.1421.5619.4618.394.3416.2813.373.938.494.103.4813.1111.3822.4213.62log10T2=4.854.344.8220.9618.7817.734.1615.6012.733.938.363.933.3112.6711.0722.4213.62log10T2=4.94.094.5320.3618.1117.114.0014.9512.143.938.263.783.1612.1910.8222.4213.62log10T2=4.953.854.2519.7817.4416.543.8714.3611.613.938.183.643.0411.6910.6522.4213.62log10T2=53.624.0019.2116.7916.023.7713.8411.173.938.123.522.9411.1710.5622.4213.62log10T2=5.053.413.7618.6716.1615.583.6913.4110.823.938.083.432.8710.6310.5622.4213.62log10T2=5.13.223.5418.1315.5515.203.6313.0710.543.938.033.362.8310.0810.5622.4213.62log10T2=5.153.033.3317.6014.9314.853.5712.8010.323.937.973.312.809.5210.5622.4213.62log10T2=5.22.853.1317.0414.3214.533.5112.5810.143.937.883.272.798.9510.5622.4213.62log10T2=5.252.682.9316.4513.6914.203.4412.409.993.937.743.242.788.3810.5622.4213.62log10T2=5.32.512.7115.7612.9813.823.3512.219.823.937.533.222.787.7710.5622.4213.62log10T2=5.352.362.5215.0312.2813.383.2412.009.623.937.253.182.757.2010.5622.4213.62log10T2=5.42.212.3014.1411.4712.823.0911.699.333.896.863.122.706.5910.4422.4213.62log10T2=5.452.062.0613.1210.5712.112.8911.278.933.796.383.032.615.9510.1822.4213.62log10T2=5.51.911.8211.969.5811.252.6610.718.423.605.812.892.485.299.7422.3113.62log10T2=5.551.781.5910.808.6110.342.4110.087.853.385.222.742.334.689.2021.9813.62log10T2=5.61.631.309.247.339.052.069.096.993.014.412.482.093.928.2920.5713.04log10T2=5.651.491.047.726.117.721.728.016.072.603.622.201.833.217.2918.5611.96log10T2=5.71.360.796.174.876.331.366.835.072.152.831.891.542.536.1716.0510.49log10T2=5.751.220.554.653.664.921.015.574.041.692.061.551.241.884.9813.188.73log10T2=5.81.090.353.222.533.550.684.293.001.231.351.210.931.283.7710.146.80log10T2=5.850.950.181.981.552.320.403.072.050.800.750.880.640.772.637.174.88log10T2=5.90.820.060.960.751.260.181.931.190.440.300.560.380.361.584.363.02log10T2=5.950.690.010.310.240.520.041.030.550.180.070.310.180.110.782.171.54log10T2=60.560.000.040.030.130.000.390.160.040.000.130.050.010.240.670.50Through Monte Carlo simulation calculation, the secondary porosity of NMR is calculated point by point in the set relaxation time search interval ofT2; based on the NMR data, the Monte Carlo method is used in the search interval of log10T2 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 log10T2 equals 5.1. The corresponding secondary porosity is 125892.54 μm, and the calculation equation of secondary porosity of NMR is as follows:
(6)ϕNMRsecondaryporosity=1.4477×ϕsecondarysurfaceporosity+0.2196.Figure 9
Linear relationship between secondary porosity and secondary porosity (16 cores) in NMR experiments atlog10T2=5.1.


## 3. Field Example

Equation (6) for calculating the secondary porosity of NMR is used when the T2 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 log10T2=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 (log10T2=5.1) is 0.391, and the average absolute error is 3.13, which indicates that the calculated secondary porosity has high accuracy.Table 4
Statistical table of errors in the calculation of secondary porosity by NMR.
No.Depth (m)Secondary surface porosity (%)NMR surface porosity (%)Surface porosity by NMR experiment (log10T2=5.1)Absolute error (%)Relative error (%)31391.268.9068.610.2963.32181404.555.58.1827.940.2422.95791409.491014.69714.780.0830.567111412.661116.14415.350.7944.920131422.0157.4587.740.2823.780181437.391623.38323.640.2571.100191441.2957.4587.220.2383.193211446.941217.59218.641.0485.957221449.45811.80111.520.2812.383Average value0.3913.131The accuracy of secondary porosity calculation can be further verified by using production performance data of well DF51. Table5 shows the productivity statistics of each sand layer in the production horizon (X3 horizon) of well DF51.Table 5
The productivity statistics of each sand layer in the production horizon (X3 horizon) of well DF51.
Sand layerEffective thicknessCores of NMR experimentNMR surface porositySurface porosity by porosity spectrumProducing pressure dropProduction splitFluid production index per meter(m-m)(m)No.Depth (m)(%)(%)(psi)(bbl/day)(bbl/psi·m· day)1412.5-14207.5111412.6616.1418.8597767.40.0796421420-14244121421.078.999.2597244.90.04563131422.011424-1428.64.6141428.523.933.1597124.60.0531071428.6-1431.83.2151430.848.037.5597231.10.0482071431.8-1435.63.8161434.393.13.559775.30.040687171435.051435.6-14404.4181437.3923.3821.1597696.70.049915According to bubble diagram (bubble size corresponds to fluid production per meter) Figure10 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.Figure 10
Intersection diagram of NMR secondary porosity and fluid production index per meter in production layer X3.Figure 11
Intersection diagram of surface porosity by porosity spectrum and fluid production index per meter in production layer X3.


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



---

*Source: 6647072-2021-04-14.xml*