A New Evaluation Method of Total Organic Carbon for Shale Source Rock Based on the Effective Medium Conductivity Theory

.


Introduction
In assessing the exploration potential of the shale oil system, the first task is to clarify the abundance, type, and maturity of the source rock organic matter [1][2][3]. The organic matter content is an essential indicator for judging the hydrocarbon generation potential of source rock. The common source rocks are mudstone or limestone rich in organic matter. In contrast, the source rocks are usually shale beds for unconventional oil and gas reservoirs [4]. The abundance of organic matter in these rocks is generally determined with pyrolysis, which is time-consuming and costly. In addition, during the experiment analysis, the selection of samples is often subjectively affected by the analyst. Because these samples are discrete, the analysis results may not accurately reflect the actual situation of underground source rock. In this case, the use of well logs to evaluate source rocks came into being [3,[5][6][7][8][9], because logging data is more available than core samples and continuous recording can eliminate the statistical errors of discrete samples.
Using logging curves to evaluate source rocks has evolved from the simple definition of shale to the combination of porosity and gamma logging to identify source rocks [10,11]. Previous studies classified source rocks based on multivariate analysis and established the petrophysical model of source rocks. They treated the organic matter as a part of the rock, and the logging curve response was a function of organic matter content [12,13]. On this basis, Passey et al. [3] proposed the quantitative calculation method of total organic carbon (TOC) of the source rocks by superimposing the sonic and resistivity curves in a reverse scale. This method is also called the Δ log R method, which realizes the high-precision evaluation of the abundance of organic matter in different lithological sections using conventional logging data. In addition, based on the logging responses to organic matter and hydrocarbons, Liu et al. [14,15] proposed the variable coefficient Δ log R method by introducing the proportional coefficient K and wave impedance, thereby improving the accuracy in predicting TOC. This method can also evaluate the free hydrocarbon content in the pores of shale source rocks. The precondition of using the Δ log R method assumes that the TOC of the formation where the resistivity and sonic curve coincide is zero, i.e., nonsource rock section. However, the natural shale formation generally contains a certain amount of organic matter. Therefore, it is necessary to correct the background value according to the regional characteristics when using the Δ log R method to evaluate the abundance of source rock organic matter.
In this work, a petrophysical model for evaluating the TOC of shale and mudstone was established by introducing the theory of effective medium conductivity. The influencing factor analysis proved that the prediction law of the model is consistent with the actual situation. Moreover, this work also presents the determination and solution methods of parameters for the new model. The newly established model is applicable to different regions by adjusting the parameters.

Geological Setting
The Songliao Basin, located in Northeast China, is the largest Mesozoic and Cenozoic continental petroliferous basin in China [16][17][18]. It is approximately rhomboid in shape, located between 42°25 ′ -49°23 ′ N and 119°40 ′ -128°24 ′ E, adjacent to the Greater Khingan Mountains, Changbai Mountains, and Zhangguangcai Mountains in the west, and adjacent to the Lesser Khingan Mountains in the north ( Figure 1). The central axis of the basin, about 750 km long, 330-370 km wide, and covering an area of 26000 km 2 , is distributed in an NNE direction, spanning parts of Inner Mongolia and the three northeastern provinces (Liu et al., 2019a). The Daqing oilfield, located in the north of the Songliao Basin, is a first-order tectonic unit, including Qijia-Gulong sag, Sanzhao sag, Daqing placanticline, and Chaoyanggou terrace (Figure 1(c)). The shale oil in the Qijia-Gulong sag of the Songliao Basin is representative in the China shale oil exploration. Up to now, a series of industrial oil flow wells have been drilled in the shale of the Upper Cretaceous Qingshankou Formation (Figure 1(d)), showing a good exploration prospect [19,20].

Samples and Methods
3.1. Samples and Experiments. The samples used for building the model are primarily from the Qingshankou Formation of the Qijia-Gulong sag. The lithology of these samples includes shale, mudstone, and siltstone. The TOC of 78 samples (cores and cuttings) was measured at the Heilongjiang Key Laboratory of Unconventional Oil and Gas Accumulation and Development, Northeast Petroleum University, China. Carbonates in the crushed samples have been removed with 10% hydrochloric acid, and all samples are heated to 50°C for 1 hour. After all the carbonates were removed, the excess hydrochloric acid solvent was wiped off with distilled water. Then, the TOC was determined on a LECO CS-230 carbon/sulfur analyzer with an experimental error of ±0.2%. In addition, the corresponding physical property data (porosity, mineral composition, etc.) and well logging data of these samples were provided by Daqing Logging Company.
3.2. The Effective Medium Conductive Model. The effective medium conductive theory (EMCT) is a macroscopic model for inhomogeneous media based on numerical and experimental analysis. It is a common theory for studying, predicting, and designing the electromagnetic response of natural and structural materials. When using the EMCT to solve electromagnetic problems, it generally depends on the electrical and magnetic properties of the constituent materials and the volume fraction of each component. Maxwell-Garnett's theory and Bruggeman's theory are two critical theories in developing the EMCT. Maxwell-Garnett's theory is a classic method to homogenize the medium, in which discrete particles are dispersed in a continuous host medium or matrix [22]. However, for mixtures with two or more components randomly distributed, Bruggeman's theory based on statistical formulas is more appropriate [23].
Maxwell-Garnett' and Bruggeman's theories assume that neither the dispersed and continuous phases conduct electricity. Hanai et al. [24] further extended the theory to the case where both the dispersed and continuous phases conduct electricity. Bussian [25] applied the theory to solve the conductivity of argillaceous sandstone. Berg [26] established an effective medium resistivity model for determining the water saturation of argillaceous sandstones. Koelman and de Kuijper [27] processed the components of the rock in a completely symmetrical form and finally established an effective medium symmetrical conductivity (EMSC) model suitable for n-component rock. The EMSC model has been introduced to describe a shaly sandstone reservoir's conductive mechanism and evaluate the reservoir oil-bearing property [28], which is also the basis for this study to assess shale's conductivity.

Four-Component Conductive
Model of Shale Rock. The EMSC model processed each component as a completely symmetrical form and considered the anisotropy for each element [27,29]. In this study, the effective medium symmetrical conduction theory is introduced to the TOC evaluation of shale rock.
It is well known that the main component of argillaceous sandstone is usually the nonconductive framework (such as quartz and feldspar). The argillaceous content is generally less than 20%, and the pores are composed of conductive formation water and nonconductive hydrocarbons, while, in shale source rocks, clay minerals and nonconductive organic matter are the main components. The content of the nonconductive framework is significantly lower than that of the reservoir. Meanwhile, the pores are mainly filled with formation water in immature-low mature source rocks, 2 Geofluids but oil and gas fill up some of the pore space as maturity increases. Solid organic matter and free hydrocarbon have similar physical properties, such as nonconduction, high acoustic slowness, and low density [30]. Moreover, the total organic carbon content primarily represents solid organic matter and free hydrocarbons. Therefore, we treated the solid organic matter and hydrocarbons as one component and defined them as the organic component in this study. Based on the above hypothesis, a four-component shale EMSC model was established. The model assumes that the shale source rock comprises the nonconductive mineral rock matrix, conductive 3 Geofluids clay minerals, nonconductive organic components, and formation water. The four-component source rock EMSC model is shown in Figure 2, and its material balance equation is shown in where V ma and V cl are the volume content of the rock matrix and clay minerals, respectively. ϕ o and ϕ w are the volume content of organic components and formation water, respectively. ϕ represents the total porosity, and ϕ ker and ϕ h are the volume content of kerogen and free hydrocarbons, respectively.
According to the EMSC theory and the derivation principle of the effective medium SATORI resistivity model proposed by the previous studies [27,29], the conductivity of the four-component source rock can be expressed as follows: where C t , C ma , C cl , C o , and C w are the conductivity of macrorock, nonconductive rock matrix, clay mineral, organic components, and water, respectively (in S/m); C 0g is the conductivity of the virtual medium (in S/m).
where ϕ 1 , ϕ 2 , ϕ 3 , and ϕ 4 are the relative contents of  4 Geofluids nonconductive rock matrix, clay mineral, organic components, and water, respectively (in decimal). It is considered that when the argillaceous content is high, the argillaceous particles are present in a continuous phase; then, λ sh ≠ 0; according to the C 0g parameterization method given by de Kuijper and Koelman [27], the result is as follows: , the calculation formula of C 0g can be expressed as follows: where h ma , h cl , h o , and h w are the geometric parameters of the rock matrix, clay mineral, organic components, and formation water in the rock, respectively (in dimensionless); λ ma , λ sh , λ o , and λ w are the percolation rates of the rock matrix, shale, organic components, and formation water in the rock, respectively (in dimensionless); and γ ma , γ sh , γ o , and γ w are the percolation indexes of the rock matrix, clay mineral, organic components, and formation water, respectively (in dimensionless) [27].

Effects of Clay
Contents on the Model. The clay content is an essential factor influencing the rock resistivity (Rt) and resistivity increase coefficient (I) (Figure 3(a)). It can be found that with the increase in clay content, the resistivity of shale gradually decreases. In contrast, when the content of clay minerals is constant, the resistivity of shale decreases as the water saturation (S w ) increases. It is due to the relative content of organic components which gradually decreases. In addition, the increase coefficient of shale (I shale ) decreases with the increase in water saturation, which is consistent with the electrical properties of porous reservoir rocks (Figure 3(b)). The feature is consistent with the actual observations, indicating that the model can describe the electrical characteristics of shale.
It is worth noting that I shale is nonlinearly related to S w , which is different from conventional reservoir rocks ( Figure 3(b)). Generally, for sandstone reservoirs with a clay mineral content of less than 10%, the resistivity increase coefficient is proportional to the water saturation in logarithmic coordinates, which is also the basis for Archie's formula to predict the water saturation of reservoir rocks [31]. However, when the clay mineral content is greater than 10%, the I shale and water saturation will show a nonlinear correlation. Previous studies believe that this nonlinear relationship is mainly due to the increase in the content of clay minerals, which leads the conductive path of the rock to become more complicated (De Witte, 1957; Poupon et al., 1954 [32, 33]). It is consistent with the simulation result in this study, in which the clay mineral content of shale rock is generally greater than 20%.
The curvature of the I shale -S w correlation line gradually decreases with the increase in clay content line, and the 75% clay content curve is a nearly linear correlation (Figure 3(b)). It may be due to the increase in clay content that has caused a change in the conductive path of the rock. For the sample with a clay content of 35%, the conduction path mainly relies on the pore water and the clay-bound water, and the conduction path is relatively complicated. By contrast, for rock with 75% clay content, its conductive mechanism is more like pure sandstone rock saturated with water.

The Influence of the Geometrical Parameters of
Nonconductive Components. The percolation rate (λ) and index (γ) are geometric parameters to describe the discrete phases in the model. The percolation rate refers to the connectivity of discrete phases, while the percolation index describes the influence of the geometry on its electrical conductivity [27,29]. For the matrix phase of nonconductive minerals, when λ ma and γ ma are constant, the electric conductivity of the nonconductive matrix is fixed, and the contribution of this component to the overall rock system mainly depends on the relative content (Figure 4), while with the increase in λ ma , the rock resistivity increases significantly, indicating that its conductivity to electricity becomes worse (Figure 4(a)). On the country, with the increase in γ ma , the rock resistivity shows a decreasing trend, but the decreasing amplitude is weaker than the change caused by λ ma (Figure 4(b)). And when γ ma is less than 3.0, its influence on the rock resistivity is more prominent. When γ ma is greater than 3.0, the influence of this parameter tends to be constant. Therefore, for the nonconductive matrix phase, the sensitivity of λ ma is significantly higher than that of γ ma , and the effective range for γ ma is between 0 and 3.0.
For the organic components, the influence of the percolation rate and index is similar to that of the mineral matrix ( Figure 5). However, since both organic components and formation water are parts of the total pores of the rock, when the rock contains no organic matter (i.e., S w = 1:0), the macroscopic electrical properties of the rock will not be affected by the organic components. In that case, the resistivity curve will converge at one point ( Figure 5).

The Influence of the Geometrical Parameters of
Conductive Components. The conductive elements considered in the model are mainly clay minerals and pore water. When the rock composition is constant, with the increase 5 Geofluids in λ cl and the decrease in γ cl , the resistivity of the rock appears to decrease ( Figure 6). As the value of λ cl increases, the change rate of rock resistivity is decreasing, approaching a particular constant value. In this case, the influence of λ cl on conductive components becomes insignificant (Figure 6(a)). Moreover, the effects of the percolation rate and index on conductive and nonconductive constituents exhibit opposite properties (Figures 4 and 6), which may be caused by the electrical conductivity of the components. In other words, in the model established in this study, for nonconductive or weakly conductive members, their connectivity (percolation rate) is the main controlling factor affecting their electrical properties. On the contrary, good conductive components' geometric properties (percolation index) are the main controlling factor that affects their electrical properties. As for the pore water part, λ w is not sensitive to the conductivity characteristics of water (Figure 7(a)), which may be caused by the water being a good conductor of electricity and the relatively low water content in shale rock. Then, the main factor that influences the conductivity characteristics of the pore water is the geometric characteristics of water, i.e., the spatial distribution form (Figure 7(b)).
In summary, the four-component shale EMSC model established in this paper can be used to characterize the electrical conductivity of shale. The percolation rate and percolation index have opposite laws for nonconductive and conductive components. For nonconductive constituents (such as mineral  6 Geofluids framework and kerogen) and weakly conductive components (clay minerals), the percolation rate and index are sensitive parameters that affect the model. On the other hand, for good electrical conductors (such as pore water), the percolation index is more sensitive to the results predicted by the model. Table 1 lists the sensitive ranges of geometric parameters for each model component obtained from theoretical analyses.   Figure 7: Effect of changes in the percolation rate (λ w ) and percolation index (γ w ) of formation water on the model.

Determination and Solution Method of
7 Geofluids tion of the mineral and porosity of the rock was based on the XRD and physical property analysis data combined with conventional logging calculation methods of clay content and porosity: (a) Clay content calculation formula is the classic natural gamma method: where V sh is the clay content, GCUR is the geochronological coefficient, SH is the clay content index, GR is the curve value of natural gamma, and GR max and GR min are the maximum and minimum values of the natural gamma curve.
(b) Porosity calculation: Using the three porosity curves (density, acoustic, and neutron) to calculate porosity is a commonly used method in well logging to determine the porosity of rocks. In this study, the density curve is selected to calculate porosity. The model treated the organic matter as a part of the total pores, but the porosity calculated by the density curve represents the porosity filled with fluid (pore water and hydrocarbons). Therefore, the measured rock organic matter content is used to correct the porosity when calculating the total porosity. The calculation formula is as follows: where Port is the total porosity, Pore org is the organic matter in pore volume, DEN is the value of the density curve, DG is the density of the rock matrix, and DF is the fluid density.
(c) Nonconductive mineral matrix content is derived from the equation of material balance of the shale model: where V ma is the volume content of the shale rock matrix, V sh is the clay content, and Port is the total porosity.

(d) Model geometric parameters:
The geometric parameters involved in the EMSC model include λ ma , λ sh , λ o , λ w , γ ma , γ sh , γ o , and γ w . Previous studies suggest defining λ ma = 1:0 and γ ma = 1:0 [27,28]. In that case, we use the measured core data as a constraint to optimize other parameters, i.e., λ sh , λ o , λ w , γ sh , γ o , and γ w . A total of 20 measured core data were used in the optimization process. The model can predict the organic component with an acceptable error (Figure 8). Table 2 shows the EMSC model parameter for the shale oil interval in well X1.

Model Solution Method.
The key to solving the ESMC model is to get the organic component content ϕ o in equation (2). Then, use the equation TOC = ððϕ o × ρ o Þ/1:25Þ/DEN to calculate the TOC, where ρ o is the average density of the organic matter component [34]. Equation (2) is an implicit function about ϕ o , so it should be solved by an iterative method. In order to ensure the convergence of the iterative process, a dichotomy algorithm is selected to solve the problem. The iterative form of the dichotomy is as follows: The following is the simulation verification and solution process. First, calculate C t from the given ϕ o , and then, based on the C t value, use the binary iterative algorithm to reverse calculating the ϕ o value under different parameter conditions. Compare the calculated ϕ o with the real ϕ o to verify the reliability of the algorithm, and examine the convergence of the model. Table 3 shows that the average relative errors are  lithology of the shale oil reservoir interval usually presents an obvious sand-shale interbed. The corresponding logging response is also very complicated in the area, which leads to poor accuracy of conventional logging evaluation TOC methods in shale oil intervals. The thickness of the shale bed in the selected interval is generally less than 5 m, and the TOC is between 0.2% and 3.2%. The porosity of the tight  9 Geofluids sandstone layer is less than 10%, and the permeability is less than 1 mD. Apart from the EMSC model, the classic △log R method [3] and the neural network method [35] were used to evaluate the accuracy. Figure 9 shows the calculated TOC results of well X1 derived from the three methods. It can be seen that these three methods can give a relatively accurate TOC value compared to the measured results. However, the TOC trends predicted by the three methods show apparent differences in the same interval, suggesting that the results cannot be judged only based on the TOC used for scaling.
Finally, to verify the generality of the new model, we selected well X2 with a conventional mudstone section from the Baiyun Sag to calculate the TOC (Figure 10). The Baiyun Sag of the Pearl River Mouth Basin is located in the deepwater slope area of the northern continental margin of the South China Sea. Source rocks of the Baiyun Sag are mainly the lacustrine and marine shales in the Wenchang, Enping, and Zhuhai Formation [36]. In this study, the main mudstone layer encountered by well X2 is the Zhuhai Formation, composed of deltaic, neritic, and bathyal sediments. The measured TOC of the mudstone interval range from 1.28% to 2.14%. Although the TOC predicted by the three methods show an excellent fitting degree with the measured TOC ( Figure 10), differences can be observed among the predicted results.

Discussion
Section 4.4 compares the application effect of the new model and the traditional TOC evaluation method. The calculated results are acceptable compared with the measured data, but the TOC trends predicted by the three methods show apparent differences in the same interval (Figures 9 and  10). We may find that the TOC curve predicted by the ESMC model is similar to that of the △log R method. It is because these two methods are based on the physical response characteristics of the rock. In contrast, the TOC curve derived from the neural network did not show apparent fluctuation except the section with scaling data.
In some tight sandstone reservoir sections (2444~2448 m and 2660~2664 m) (Figure 9), the sonic time curve and the resistivity curve are approximately coincident, which leads to the TOC value predicted by the △log R method being lower than the measured value. Still, the model established in this study has achieved better prediction results. Moreover, at the calcareous shale between 4430 m and 4440 m, the TOC calculated by the new model is slightly higher. In contrast, the value calculated by the △log R method is lower ( Figure 10). It is mainly due to the calcium mineral in this shale section which causes a relatively low acoustic slowness, and the resistivity curve shows a relatively high value.

Geofluids
Finally, these two curves are close to coincidence. Therefore, the TOC value predicted by the △log R method is relatively low. By contrast, the EMSC model is a four-component model and does not consider the calcium content of the rock. Thus, the high resistivity of calcium leads to a false high TOC. Although the neural network method has the highest prediction accuracy, the predicted TOC value is unacceptable in the two cases of this study. It is inconsistent with the results of both the model established in this study and the △log R method. We inferred that the neural network method generally needs a large number of samples to support the modeling. However, the measured TOC value of well X1 has a small range (1.29%~2.14%), which makes the neural network unable to extract the sufficient characteristic values of different lithological conditions effectively.

Conclusion
(1) Based on the effective medium symmetrical conductivity theory, a logging evaluation model of the total organic carbon of the rock was established. The new model is a four-component petrophysical model, which mainly considers the nonconductive rock skeleton, conductive clay, organic components (kerogen and hydrocarbon), and pore water (2) Compared with the classic △log R and neural network methods, the new model considers the influence of rock mineral composition, component structure, and shape on rock conduction. Therefore, the new TOC evaluation model is suitable for intervals with a large section of shale rock, and lithology frequently changes intervals (3) In the actual application process, the EMSC model has achieved high evaluation accuracy, which is comparable with the △log R method and neural network method, confirming the practical application value of this method. However, it is worth noting that the model will cause more parameters to be involved in the application process, which requires systematic analysis and test data as support

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.