TOC estimation of shale oil reservoir by combining nuclear magnetic resonance logging and nuclear physics logging

The evaluation of source rock properties has become a vital step in logging interpretation. Total organic carbon (TOC) content is the key to estimating the quality and hydrocarbon generation potential of source rocks. In the shale oilfield of the Junggar Basin, the conventional method of calculating the TOC of hydrocarbon source rocks cannot satisfy logging evaluation requirements. This paper predominantly deals with a method for the quantitative estimation of TOC in source rocks via nuclear physics and nuclear magnetic resonance (NMR) logs. According to this method, the total hydrogen index of the source rock is the sum of the response of kerogen, clay minerals and fluid, expressed by corrected neutron porosity. The hydrogen index of fluid and clay minerals is indicated by the effective porosity of NMR and the estimated clay content, respectively. To eliminate the hydrogen index of fluid, the effective NMR porosity is subtracted from the corrected neutron porosity. On this basis, a new and overlapping method suitable for clay-rich rocks and oil reservoirs is proposed. This method was developed by overlaying the scaled clay content curve on the hydrogen index curve. In non-source rocks, the two curves regularly overlap. However, in organic-rich rocks the two curves will separate. The separation distance between the two curves was used to estimate TOC continuously. Possessing sound application and benefiting from the measured results of sweet spots, this method provides new insights for TOC quantitative prediction in shale oil reservoirs.


Introduction
Shale oil reservoirs, a significant field of unconventional oil and gas exploration, are an important replacement resource of oil and gas for the world of tomorrow (Stevens et al. 2013;Madi et al. 2015). However, shale oil reservoirs possess strong heterogeneous rock physical properties as well as source-reservoir integration characteristics (Zou et al. 2013). Previous studies show that brittleness, in situ stress anisotropy, quality of source rock combined with the 'four properties' of conventional reservoirs (comprising lithology, physical property, electrical property and oiliness) form the 'seven properties' for shale oil reservoirs (Kuang et al. 2021).
Evaluation of source rock properties includes the abundance, type and maturity of organic matter. These jointly control and influence the hydrocarbon generation potential of source rocks. Total organic carbon (TOC) plays an important role in estimating the abundance of organic matter. Organic matter can be obtained through organic geochemical analysis of core samples. The traditional experimental method, however, is expensive and time-consuming, and produces results that do not explain TOC's vertical and horizontal variation trend in the reservoir (Xu & Zhu 2010).
To overcome the problem of limited core samples and to further obtain continuous TOC evaluation results, scholars have carried out numerous studies on TOC estimation based on logging data. These mainly involved the single curve method, natural gamma-ray logging and density logging curves, as well as the multiple regression, carbon/oxygen logging, curve overlap, ΔlogR and matrix density overlap methods. The advantages of the single curve method are its basic principles and convenient operation. However, it is compromised by environmental factors such as well diameter, fluid and lithology (Sehmoker 1981;Fertl & Chilingar 1988). Usually, TOC estimation results calculated in this way are inaccurate. Compared with the single curve method, the multiple regression method considers more conventional logging information and is more accurate. However, more regression parameters mean more complex weight allocation (Tan et al. 2021). Carbon/oxygen logging can exclude the influence of formation water salinity (Zhao et al. 1994) but is not suitable for complex mineral types of reservoirs. The ΔlogR method establishes the relationship among logging curve, organic matter abundance and kerogen maturity (Passey et al. 1990). When the resistivity curve changes frequently, however, it is difficult to determine the baseline position and maturity parameters. To enhance the applicability of the ΔlogR method and improve prediction accuracy, some scholars have refined the model form and parameters (Wang et al. 2020a;Zhu et al. 2019). The matrix density overlap method proposed by Jacobi et al. (2008) can effectively eliminate the influence of inorganic carbon. However, it needs elemental capture spectroscopy logging data. Ramirez et al. (2011) proposed a method to estimate TOC by subtracting the density porosity from nuclear magnetic total porosity based on the difference between the measurement principles of nuclear magnetic and density logging in shale reservoirs. However, complex mineral types can lead to instability of skeleton density parameters. Zhao et al. (2016) proposed a new method based on the overlap of natural gamma ray and clay indicator parameters, making up for the deficiency of the ΔlogR method in overmature source rocks. It is not suitable, however, for continental shale reservoirs with high potassium content. In addition to the previous mentioned evaluation methods, artificial intelligent algorithms such as artificial neural networks (Mahmoud et al. 2017), the support vector machine (De et al. 2019), extreme learning machine (Shi et al. 2016), Gaussian sequence simulation (Nezhad et al. 2018) and the comprehensive application of a variety of intelligent algorithms (Zhu et al. 2020) have been developed to evaluate TOC. While highly precise, these intelligent algorithms maintain high requirements on data sets but lack the theoretical guidance of physical models.
Based on the physical volume model of source rocks and the typical logging response characteristics of kerogen, this paper proposes a TOC estimation method combining nuclear physics logging (Density logging (DEN), compensated neutron logging (CNL) and gamma logging) and nuclear magnetic resonance (NMR) logging. The two sweet spots interval core data of the Permian Lucaogou Formation shale oil reservoir in the Jimsar Sag are applied, showing a comparison between estimated results and experimental data. The comparison demonstrates that the method is reliable in the TOC quantitative prediction of a shale oil reservoir.

Geological background
As shown in figure 1a and b, the Jimsar Sag, located in the southwest of the eastern uplift belt of Junggar Basin, is bounded by thrust faults in the south, west and north of the sag ( Jiang et al. 2015;Cao et al. 2016). The Jimsar Sag has experienced Hercynian, Indosinian, Yanshan, Himalayan and other multi-stage orogenies, and closed in the late Middle Permian, forming a relatively independent lacustrine basin environment (Guo et al. 2019). As the main production layer, the Lucaogou Formation experienced multiple periods of lacustrine regression and transgression during the sedimentary period (Qiu et al. 2016a). As shown in figure 1c, the lithology of the Lucaogou Formation is diverse, characterized by the interbedded distribution of thin bedded clastic and chemical rocks, including arenaceous dolomite, feldspathic litharenite, and mudstone, as well as a small amount of limestone (Kuang et al. 2015;Qiu et al. 2016b;Su et al. 2019).
According to lithologic characteristics and logging curves, Lucaogou Formations can be divided into the upper and lower Lucaogou Formation (P 2 l 2 , P 2 l 1 ). There are two 'sweet spots' (P 2 l 2 2 , P 2 l 1 2 ) composed of the interbedding of fine-silty siltstone, carbonate rock and organic-rich mudstone. Source rocks of Lucaogou Formation are high quality, with mature source rocks widely distributed in the Jimsar Sag possessing strong hydrocarbon generation potential (Kuang et al. 2014;Qiu et al. 2016a). Organic geochemical data show that the kerogen type of shale reservoir in the Lucaogou Formation is mainly type II; siliceous mudstone and dolomitic mudstone are predominantly type I and type II, and silky mudstone is type II and type III (Cao et al. 2016;Qiu et al. 2016b). The source rocks have high organic matter abundance: the floor of TOC is >1% with the average TOC reaching >4%. The organic matter abundance of mudstone and dolomite is higher than that of silty mudstone (Xiang et al. 2013). At present, the source rocks of 834 Figure 1. Location of the study area (Cao et al. 2016). Lucaogou Formation are in the mature stage of the 'oil generation window'-vitrinite maturity is between 0.7 and 1%-showing high production potential ( Jiang et al. 2015;Su et al. 2019).

Rock physical volume model
For shale reservoirs, scholars maintain that organic matter is enriched only in source rock such as mud shale or calcareous mudstone, while only a small or negligible amount exists in other lithologic reservoirs (Wang & Guo 2000;Liu 2009). Figure 2 shows the differences between non-source and source rocks in rock physical volume models (Passy et al. 1990;Lei et al. 2010). Non-source rocks mainly consist of non-clay minerals, clay minerals, pore water and hydrocarbon. In addition to these four elements, the source rocks also contain some kerogen. The evaluation of TOC by logging curve is based on the typical logging response characteristics of kerogen (figure 3), summarized as 'four high and one low' . 'Four high' refers to high natural gamma ray, neutron porosity, interval transit time and resistivity, while 'one low' refers to low density.
Previous studies on the physical and chemical properties of kerogen have concluded that kerogen is rich in hydrogen and has high radioactivity and low conductivity. Its den- sity and acoustic wave propagation speed is lower than common minerals (Wang & Guo 2000;Lewis et al. 2004;Ward 2010). According to the physical volume models of source rocks, the resistivity curve is affected by many factors such as the mineral type of rock matrix, pore fluid type, pore structure, wettability and oil-bearing properties. Thus, the logging response mechanism is complicated. The three porosity curves (density, neutron porosity and acoustic time) are the comprehensive response of rock matrix, kerogen and pore fluid. The natural gamma-ray logging signal mainly comes from the radioactivity of the source rock matrix. The organic matter of source rock can absorb uranium ions in the 835  formation and increase the radioactivity, the result of which is the natural gamma-ray curve reflecting the influence of organic matter. However, previous studies suggest that volcanic hydrothermal activity is an important geological origin of Permian Lucaogou Formation source rocks (Zhang 2019;Ding et al. 2019;Wang et al. 2020b). Therefore, there are many tuffaceous components composed of alkaline feldspar and quartz derived from alkaline magma in the Lucaogou Formation ( Jiang et al. 2015;Zhang 2019), resulting in poor response sensitivity of natural gamma-ray logging to organic matter (figure 4).

Theory and methodology
According to logging data collected so far, there is a lack of data on natural gamma spectrometry. However, high-quality NMR logging data is sufficient. Consequently, this study pro-poses an estimation method of TOC in source rocks based on nuclear physics logging (density, neutron and gamma logging) combined with NMR logging.
The compensated neutron logging principle dictates that the formation deceleration for neutrons is largely dependent on the hydrogen content (Huang et al. 2008). Compensated neutron logging converts the thermal neutron counting rate into hydrogen index (HI) or neutron porosity by calibrating in water-bearing limestone calibration test pits. In the shale oil reservoir of the Lucaogou Formation, the main fluid types include free oil, absorbed oil, free water and bound water. The hydrogen index of liquid hydrocarbons, however, is not at all dissimilar to that of fresh water (HI = 1), resulting in a negligible influence of fluid property change on compensated neutron logging in shale oil reservoirs. According to the XRD analysis data in Table 1, the common minerals of the Lucaogou Formation are primacy dolomite, calcite, quartz and plagioclase, followed by potash feldspar and clay minerals (mainly montmorillonite and illite). The compensated neutron logging responses of these common minerals are shown in Table 2 (Schlumberger 1998;Yong & Zhang 2002). Except for clay minerals, the compensated neutron logging response values of other minerals (dolomite, calcite, quartz, plagioclase and potash feldspar) are very small. In summary, the compensated neutron logging response of source rocks can be simplified into a volume model composed of clay minerals, kerogen and pore fluid, as where Φ N is the neutron log response (total hydrogen index) of shale reservoir, percentages of Φ k , Φ cl and Φ f are the hydrogen indices of kerogen, clay mineral and pore fluid, respectively; V k , V cl and Φ are the volume specific gravity of 836 kerogen, clay minerals and porosity, respectively, as a percentage.

The total hydrogen index of shale
In Table 2, the neutron deceleration ability of quartz and feldspar minerals is weaker than that of calcite, resulting in an equivalent hydrogen index of <0. On the contrary, the neutron deceleration ability of dolomite is stronger than that of calcite, resulting in the equivalent hydrogen index >0. To further improve the sensitivity of the compensated neutron log-ging to the kerogen (V k Φ k ) in equation (1), and to reduce the influence of other minerals in the rock matrix on the measurement results, the same calibration standard of CNL and DEN logging is used to correct the neutron log response of the shale reservoir. It is well established that CNL and DEN logs are calibrated in water-saturated limestone standard wells. When the same calibration is used for both, the CNL and DEN log curves of limestone overlap, as would be the case if neutron porosity ranged from 45 to −15% and density ranged from 2.95 to 1.95 g/cm 3 (Yong & Zhang 2002). Dolomite has strong neutron deceleration ability, leading to high neutron porosity. However, its matrix density is higher than that of limestone, resulting in low density porosity calculated via the volume model. Similarly, sandstone (quartz and feldspar minerals) has a weak neutron deceleration ability, leading to low neutron porosity. Its matrix density is lower than that of limestone, resulting in higher density porosity calculated via the volume model (Huang et al. 2008).
Based on this, the average value of neutron porosity and apparent density porosity is taken to reduce the influence of minerals on the total hydrogen index of shale. In addition, because density logging instruments consider kerogen response as part of the pore response, density porosity contains both kerogen and pore fluid effects . Therefore, the aforementioned hydrogen index correction method does not greatly influence the hydrogen index of kerogen and pore fluid. Equations (2) and (3) are used to correct the total hydrogen index of shale (the Φ N of equation (1)): where Φ is the neutron log response (total hydrogen index) of shale after correction (the correction of Φ N ), as a percentage; Φ D is apparent density porosity (%); Φ N is neutron porosity (%); b is the density logging value (g/cm 3 ); ma is the matrix density of limestone (2.71 g/cm 3 ) and f is the density of the pore fluid (1.0 g/cm 3 ).

The hydrogen index of pore fluid
NMR logging can also measure the hydrogen index but the measurement objects are different from compensated neutron logging. As shown in Table 1, paramagnetic minerals such as pyrite are rarely developed in the source rocks. We can regard the NMR response as unaffected by rock skeleton (Deng & Xie 2010;Feng et al. 2017Feng et al. , 2021. On the other hand, as hydrogen nuclei of organic matter are often attached to pore surfaces in solid form, their transverse relaxation time is very quick in detecting their signals via NMR logging (Ramirez et al. 2011). Thus, the response of NMR logging is almost unaffected by the solid skeleton and organic matter of source rocks.
Previous studies have shown that the NMR response characteristics of clay minerals depend on the form and type of clay (Prammer et al. 1996;Jia 2017). Hydrogen in clay minerals mainly takes three forms: OH − of the clay structure (structural water), OH − at the edges of clay platelets (interlayer water) and OH − on the surface of clay particles (adsorbed water). As the T 2 time of the first two forms is <0.1 ms, NMR logging failed to obtain relevant information (Fleury & Romero 2016). The third form can be calculated by ensuring the cutoff value of the T 2 distribution. Hence, the T 2 spectrum contains a portion of hydrogen nuclei from the contribution of clay minerals. To whittle down the efforts of clay minerals on the hydrogen index of pore fluid, the clay-bound water portion needs to be separated from the T 2 spectrum. This means that the calculation of NMR effective porosity (ɸ e in figure 5) is required.
There are two methods to obtain the value of ɸ e . The first is experimentally measuring the effective porosity of the core directly and the corresponding T 2 spectrum. The second is to compare the porosity calculated by NMR logging at different T 2 times with the effective porosity analysis data of the core. Wang et al. (2019) achieved an accurate calculation of the effective porosity of shale oil reservoirs in the Lucaogou Formation of the Jimsar Sag by using an iterative method and determined that the start time to obtain the effective porosity was 1.7 ms. Effective porosity can be obtained by summing up the porosity components of the T 2 spectrum after 1.7 ms, as where ɸ e is the the NMR effective porosity (%) and ɸ i is the porosity component corresponding to different T 2 times (%). Due to the hydrogen index of liquid hydrocarbons being very close to that of fresh water (Φ f = 1), the meaning of ɸ e can be regarded as the hydrogen index of pore fluid (ɸΦ f ).

The hydrogen index of clay mineral
Combining the forms equations (3) and (4)   see the sum of hydrogen index contributed to by clay minerals and kerogen, as To further discern the contribution of kerogen to the hydrogen index (V k Φ k in equation (5)), it is necessary to use other log curves and parameters to imply the clay minerals content and separate it from the residual hydrogen index. Previous studies have summarized the distribution range of neutron porosity and density of four common clay mineral types (Schlumberger 1998;Mao 2001;Yong & Zhang 2002). These studies showed that the subtracted result between neutron porosity and density porosity (Φ ND ) of clay minerals generally shows a significant positive difference and is much larger than other mineral types. The average density of kero-gen was taken as 1.175 g/cm 3 (Yao 2019), with the neutron porosity of type I and type II kerogen usually at 86 porosity unit (Zhao et al. 2015). The estimated Φ ND of kerogen is ∼3.54%, differing considerably from that of clay minerals in Table 2. In conclusion, Φ ND is an important sensitive parameter to characterize clay minerals content: where Φ ND is the subtracted result between neutron porosity and density porosity (%) that will be positively correlated with clay content. In addition, clay-bound water is related to clay content and type. Thus, ɸ e is negatively correlated with clay content (Deng & Xie 2010). Due to the lack of a natural gamma-ray spectral log, a multivariate statistical regression 839 Figure 7. TOC was estimated by overlapping Φ k+cl curve and V cl curve. model for quantitative prediction of clay minerals content was established by using natural gamma ray, neutron porosity, density, acoustic time difference, ɸ e and Φ ND (He et al. 2020). As shown in figure 6, the single factor correlation analysis above these parameters and the clay minerals content obtained by the XRD experiment was carried out. From figure 6, the correlations between the clay content and Φ ND , ɸ e and GR are relatively high, while the correlations between the clay content and AC, CNL and DEN are relatively low. Ultimately, we established a multivariate statistical regression model for clay content using Φ ND , ɸ e and GR, as shown in equation (7), with correlation coefficient is 0.8388: where V cl is the estimated clay content (%).

The new overlap method
For the sake of description, Φ-ɸ e in equation (5) is defined as parameter Φ k+cl . This includes the contribution of kerogen and clay minerals to the hydrogen index. V cl is positively correlated with the contribution of clay minerals to the hydrogen index. In the non-source rock segment both Φ k+cl and V cl mainly reflect the contribution of clay minerals, rendering their curves very similar. However, in the source rock segment Φ k+cl not only reflects the contribution of clay minerals but the additional contribution of kerogen is also seen. The two curves will present different variations. Therefore, the curve scale range is adjusted appropriately until the Φ k+cl curve and V cl curve overlap in the non-source rock segment, bringing about the separation of the two in the source rock segment ( figure 7).
At this point, the quantitative relationship can be established between TOC and the separation distance of Φ k+cl curve and V cl curve in source rock, as shown in equations (8) and (9) as follows: where Δs is the distance between the Φ k+cl curve and the V cl curve, indicated by the yellow area in figure 7 (%); (Φ k+cl ) min and (Φ k+cl ) max are the minimum scale and maximum scale of the Φ k+cl curve, respectively; (V cl ) min and (V cl ) max are the minimum and maximum scales of V cl curve, respectively, and k and b are the slope and intercept of TOC quantitative evaluation model, respectively. The details of the calibration of the values of k and b in equation (9) are as follows: first, as shown in equation (8), (Φ k+cl ) min , (Φ k+cl ) max , (V cl ) min and (V cl ) max were used to normalize the Φ k+cl curve and V cl curve respectively; next, based on the normalized Φ k+cl curve and normalized V cl curve, the distance Δs was obtained by deducing the difference; finally, the values of k and b in equation (9) are calibrated using linear regression between Δs and the actual measured TOC.

Application and verification
The aforementioned methods are applied to the actual well data, with application results shown in figure  V cl curve, while the eighth track includes the TOC predicted by the ΔlogR method (the dotted green line) and the actual measured TOC (red dot). The final track includes the TOC estimated by the proposed method (the solid black line) and the actual measured TOC.
The experimental measurement process of TOC in cores conforms to the corresponding Chinese national mandatory standard (GB/T 19145-2003). The general experimental process is as follows: first, the sample was crushed into particles <80 mesh before diluted hydrochloric acid was added to remove the inorganic carbon of carbonate minerals; next, the particles were combusted in a high-temperature oxygen stream after drying, releasing CO 2 ; finally, the LECO carbon analyzer was used to measure the TOC.
As shown in figure 8b, the section below 'XX60m' is the non-source rock zone. This can be distinguished from the source rock zone in two respects. First, the TOC measured in the experiment is low (<1%), indicating that the abundance of organic matter in this section is very low, with the thickness of this section stable. Second, the log curves of source rocks above 'XX60m' are dentiform, with logs characterized by higher natural gamma, resistivity, neutron porosity and lower density than those of non-source rocks.
According to the non-source rock section of the lower member of Well JI A, the scale range of Φ k+cl curve and V cl curve were adjusted to ensure they basically overlap. The confirmed values of (Φ k+cl ) min and (Φ k+cl ) max are 0 and 50%, with (V cl ) min and (V cl ) max at 0 and 100%. As seen from the final track in figure 8, the estimated results are in line with the actual measured TOC. In section 'XX42∼XX52m' in figure 8a, the measured TOC of the cores is very low. However, the TOC predicted by the ΔlogR method is high. In contrast, section 'XX20∼XX28m' in figure 8b shows, measured TOC of the cores is very high with low TOC predicted by the ΔlogR method. TOC predicted by the ΔlogR method is close to the overall trend of measured TOC in cores. However, the new method is much more accurate. Figure 9 parts a and b are the cross plots of actual measured TOC and the distance between the Φ k+cl curve and V cl curve (Δs) in the two sweet spots of Well JI A, with correlation coefficients of 0.797 and 0.846, indicating that TOC and Δs maintain sound linear correlation. As shown in figures 8a and 9a, the equation to estimate the TOC of the upper sweet spot of Well JI A can be expressed as TOC(%) = 0.311 × Δs + 0.0312, R 2 = 0.797.
In the same way, as shown in figures 8b and 9b, the equation to evaluate the TOC of the lower sweet spot of Well JI A can be expressed as Table 3 sums up the forecasting statistics from the two sweet spot formations for contrasting the evaluated outcome with the experimental data. The average error and root mean squared error (RMSE) of the estimated TOC and that of sealed coring measurement are ∼10 and 3%, respectively, showing sound application effect. Using figure 9 and Table 3, the feasibility and accuracy of this novel way of shale oil reservoirs in the Lucaogou Formation were confirmed. As shown in figure 10, to further verify the extensibility of equations (10) and (11), the measured TOC of another well (Well JI B) was selected to compare with the 842 Figure 10. Cross plots of measured TOC in cores and TOC estimated by models of the upper sweet spot (a) and the lower sweet spot (b) in Well JI B. model estimation results. The result shows that the slopes of the upper and lower sweet spots are 0.894 and 0.904 (close to a diagonal distribution), with the correlation coefficient at 0.780 and 0.862, respectively, proving that the model estimation results correspond with the experimental results.

Conclusion
Based on accumulated knowledge on the rock physical volume model of hydrocarbon source rocks and neutron response characteristics of common minerals and fluids in sedimentary rocks, this paper proposes a new overlapping method combining NMR and nuclear physics logs to estimate the TOC.
First, considering the diversity of shale minerals, neutron porosity is corrected by density porosity to obtain the final total hydrogen content index (Φ); second, the T 2 spectrum after the optimal T 2 cutoff value accumulated to obtain the effective porosity of NMR (ɸ e ), characterizing the contribution of pore fluid to the hydrogen index. The contribution of pore fluid to the total hydrogen index of shale is eliminated by the difference between Φ and ɸ e , with the remaining hydrogen index (Φ k+cl ) forming the sum contribution of clay minerals and kerogen. This study developed a suitable estimated model of clay minerals content (V cl ) via the multivariate statistical regression method, the correlation coefficient of which is 0.8388. Except in the case of clay minerals, in source rock layers, the difference between the Φ k+cl curve and the V cl curve will expand as Φ k+cl reflects the contribution of kerogen. Following calibration with core data, the separation distance between the two curves was used to evaluate TOC sequentially.
This new method was successfully used in the sweet spots of a single well. The comparison of the new method with the TOC measured by the actual core shows that the new method achieves sound results. The coefficients of correlation between the two sweet spots estimated TOC and the actual measured data were 0.797 and 0.846. The average errors and RMSE between the two sweet spots estimated TOC and the experimental data were <10 and 3%, respectively. In addition, the introduction of another well and coring data validated the extensibility of the new method. This method corrected the influence of some matrix minerals on TOC evaluation results, overcoming the common failure to estimate TOC with abnormal RT and GR. This method is of great significance for the continuous evaluation of the hydrocarbon generation potential of shale oil reservoirs.