The Saturation Calculation of NMR Logging Based on Constructing Water Spectrum Function

Abstract

Tight sandstone oil reservoirs are characterized by complex structures, poor pore connectivity, and strong heterogeneity, with features such as low porosity and ultra-low permeability. Conventional methods for calculating saturation cannot accurately evaluate the hydrocarbon saturation of these reservoirs. To address this, a study was conducted from the perspective of non-electrical logging methods, focusing on the inherent nuclear magnetic resonance (NMR) characteristics of different fluids to develop a saturation calculation method that avoids the influence of the rock matrix, thus enabling precise saturation measurement in tight sandstone oil reservoirs. The traditional NMR porosity model was modified by segmenting it using the clay-bound water cutoff value, aiming to identify the distribution pattern of fluids in pores outside the clay-bound water zone. Through theoretical derivation and water spectrum function simulation, a water spectrum function and its parameter range suitable for the NMR T2 distribution in tight sandstone reservoirs were determined. Using core-sealed core saturation as a reference, the particle swarm optimization (PSO) algorithm was applied to optimize the parameter range and construct the final water spectrum function tailored to tight sandstone oil reservoirs. The accuracy and practicality of this method were validated by applying the derived water spectrum function to NMR logging in the exploration block, allowing for precise saturation calculations and the accurate evaluation of tight reservoir saturation.

1. Introduction

Tight oil, as an unconventional oil and gas resource, has garnered significant attention globally. Compared to conventional reservoirs, tight oil reservoirs are characterized by complex structures, poor pore connectivity, and strong heterogeneity, with features such as low porosity and ultra-low permeability. These reservoirs require large-scale fracturing techniques or horizontal multi-branch well technologies to achieve commercial oil and gas flow [1,2,3,4]. Traditional evaluation models and calculation methods are insufficient for effectively assessing these reservoirs, necessitating modifications to traditional models and the adoption of new technologies and methods to obtain accurate and reliable evaluation parameters. Due to the characteristics of hydrocarbon migration and accumulation, the pore structure of the rock plays a key role in controlling fluid flow capacity and accumulation [5,6,7]. The quality of the pore structure also influences the distribution of different fluids within the pores. By considering the relationship between pore structure and saturation, it is possible to evaluate the hydrocarbon potential of tight reservoirs.

To date, the most commonly used method for quantitatively evaluating reservoir fluid saturation through logging is by establishing a resistivity saturation model [8]. The principle of this method is based on the fact that when the reservoir pores are saturated with oil, natural gas, or formation water, the reservoir exhibits different conductivities, resulting in variations in resistivity. By using Archie’s formula, oil saturation in the reservoir can be accurately calculated, thus enabling the identification of oil, gas, and water layers. By considering factors such as irreducible water saturation, pore structure, and clay content, a series of resistivity-based saturation calculation methods have been derived by combining Archie’s formula with the parallel conductivity model [9,10,11]. Some researchers have attempted to calculate saturation using NMR porosity combined with resistivity [12], but the results are consistent with traditional resistivity saturation models, remaining within the scope of resistivity-based models.

Non-resistivity-based saturation models have also seen some development. Methods such as the NMR T2 spectrum differential method, shift method, TDA technique using time-domain difference, and water spectrum construction from echo train signals have shown qualitative identification capabilities for certain reservoirs [13,14], but these methods lack universality and cannot perform quantitative saturation calculations [15]. Additionally, some researchers have attempted to observe and calculate saturation decay through NMR experiments [16], with the feasibility of the method resting on the study of saturation decay patterns in different pore structures. However, discrepancies between the experimentally obtained saturation decay and the actual pore structure saturation changes in reservoirs limit its current application, and it remains in the experimental stage.

Furthermore, big data, artificial intelligence, and machine learning methods have been widely applied in the classification and regression of reservoir parameters using geophysical logging curve information [17,18,19,20,21]. Some researchers have used these methods to predict saturation, but due to challenges such as the difficulty in eliminating well data errors and the lack of standardized data, the accuracy of saturation identification remains low. Although these methods can identify fluid properties under specific geological conditions to some extent, they still lack universality and cannot replace the role of resistivity-saturation models in reservoir fluid identification.

The tight oil reservoirs of the Permian Wuliji Formation in the Jimusaer Depression are characterized by low porosity and low permeability [22]. The effect of the reservoir framework on resistivity is greater than that of the reservoir fluids, making it difficult for conventional resistivity saturation models to accurately reflect the true situation when calculating reservoir saturation. To effectively evaluate the oil-bearing potential of tight reservoirs, it is necessary to accurately calculate oil saturation. However, due to the inherent limitations of resistivity-saturation models, such errors are difficult to eliminate, necessitating the use of non-electrical methods for saturation calculation to achieve precise results. NMR logging has demonstrated some effectiveness in identifying reservoir oil content. Through the study of NMR porosity distribution models, it was found that fluid distribution in the pores exhibits a certain regularity, with oil and water occupying different pore locations, i.e., different intervals of the T2 spectrum [23]. Due to the coexisting states of fluids within the pores (e.g., oil and water), the T2 spectrum only shows changes in the spectrum shape rather than distinct peaks representing different fluids. Conventional NMR porosity models only consider the effect of the pores and do not account for fluid changes, failing to reflect the fluid characteristics in the reservoir.

Therefore, the NMR porosity model was first improved, and a saturation calculation method based on the constructed water spectrum function was proposed, allowing fluid separation in the reservoir through NMR logging. This method is applicable for analyzing the oil-bearing potential of tight sandstone reservoirs. The improved NMR porosity model and saturation calculation method not only enhance the accuracy of oil-bearing evaluations but also provide a more effective and precise tool for oilfield development and assessment under complex geological conditions. The combination of measured NMR logging data has shown excellent application results in evaluating the tight oil reservoirs of the Wuliji Formation in the Jimusaer Depression.

2. Improvement of NMR Logging Porosity Model

Nuclear magnetic resonance (NMR) logging porosity is typically divided using cutoff values, generally by employing the clay-bound water cutoff value (T_2B) and the bound water cutoff value (T_2cutoff). The bound water cutoff value is the general cutoff value, which determines the cutoff value for movable fluid porosity, while the clay-bound water cutoff value determines the cutoff value for effective porosity [14,24]. As shown in Figure 1, from left to right, it includes rock matrix, dry clay, clay-bound water porosity, capillary-bound water porosity, and movable fluid porosity (movable water porosity and movable hydrocarbon porosity), segmented using cutoff values:

$$MBVI={\varphi }_b={\displaystyle {\int }_{T_{2\min }}^{T_{2cutoff}}por(T_2)d{T}_2}$$

(1)

$$MFFI={\varphi }_m={\displaystyle {\int }_{T_{2cutoff}}^{T_{2\max }}por(T_2)d{T}_2}$$

(2)

$$MPHI={\varphi }_e={\displaystyle {\int }_{T_{2B}}^{T_{2\max }}por(T_2)d{T}_2}$$

(3)

$$MSIG={\varphi }_t={\displaystyle {\int }_{T_{2\min }}^{T_{2\max }}por(T_2)d{T}_2}$$

(4)

MBVI represents bound fluid porosity; MFFI represents movable fluid porosity; MPHI represents effective porosity; MSIG represents total porosity.

The traditional nuclear magnetic resonance (NMR) porosity model divides the pore fluids into different components using cutoff values. However, the actual distribution of fluids in geological formations cannot be entirely determined by cutoff values. It is assumed that besides clay-bound water, fluids in the effective pore space exist in an overlapping manner between pores. Using empirical parameters from the Satuo block (clay-bound water cutoff value of 2.5 ms, movable fluid cutoff value of 23.53 ms), the clay-bound water cutoff value T2B = 2.5 ms is used for porosity model segmentation. As shown in Figure 2, this improved model schematic can be understood as follows: in the clay-bound water pore portion, water accounts for 100, i.e., weight equals 1, whereas outside the clay-bound water pores, the distribution of fluids follows a certain regularity of weight distribution. That is, the weight distribution of oil–water distribution can satisfy a certain functional relationship, thus proposing the method of coefficient construction to build the water spectrum.

As shown in Figure 3, the basic process for calculating saturation using nuclear magnetic resonance (NMR) water spectrum is illustrated. Figure 3c represents the T2 distribution in the NMR logging (Figure 3a), and it is weighted with the water spectrum function shown in Figure 3b. The distribution of weights and the weighting process are depicted in Figure 3d. Finally, the water spectrum and oil spectrum are obtained, as shown in Figure 3e.

According to the improved porosity model, the saturation calculation formula for the water spectrum is as follows:

$$S_w={\displaystyle \frac{{\displaystyle {\int }_{T_{2\min }}^{T_{2B}}por(T_2)d{T}_2}+{\displaystyle {\int }_{T_{2B}}^{T_{2\max }}S\left(T_2\right)por\left(T_2\right)d{T}_2}}{{\displaystyle {\int }_{T_{2\min }}^{T_{2\max }}por\left(T_2\right)d{T}_2}}}$$

(5)

In which S(T₂) represents the water spectrum function, which is defined as follows:

$$\{\begin{array}{l}S(T_2)=1,T_2\le T_{2B}\\ S(T_2)={\displaystyle \frac{a_1-a_2}{1+{\left(\frac{T_2}{T_{2CW}}\right)}^m}}+a_2\end{array},T_2>T_{2B}$$

(6)

In the equation, T_2CW represents the oil–water line, which can be understood as the T₂ value when oil and water each occupy 50% of the pore fluid, and a₁, a₂, and m are function parameters.

Calculating reservoir petrophysical parameters based on the improved nuclear magnetic resonance porosity model:

$$MBVI={\varphi }_b={\displaystyle {\int }_{T_{2\min }}^{T_{2B}}por(T_2)d{T}_2}$$

(7)

$$MFFI={\varphi }_m={\displaystyle {\int }_{T_{2B}}^{T_{2\max }}por\left(T_2\right)d{T}_2}$$

(8)

$$MFWI={\varphi }_w={\displaystyle {\int }_{T_{2B}}^{T_{2\max }}S\left(T_2\right)por\left(T_2\right)d{T}_2}$$

(9)

$$MFOI={\varphi }_o={\displaystyle {\int }_{T_{2B}}^{T_{2\max }}\left(1-S\left(T_2\right)\right)por\left(T_2\right)d{T}_2}$$

(10)

MBVI stands for bound fluid porosity; MFFI stands for movable fluid porosity; MFWI stands for movable water porosity; MFOI stands for movable oil porosity.

The construction of nuclear magnetic resonance (NMR) water spectrum involves extracting the T₂ distribution corresponding to the depth of saturation from the NMR logging, aligning it with the true saturation of the core, and then determining the water spectrum function. The basic idea is to weigh the extracted NMR logging T₂ distribution with the water spectrum function determined by different water spectrum parameters (i.e., multiplying the T₂ distribution by a weighting function) to obtain a series of water saturation values. Then, based on the true water saturation of the labeled core samples (obtained from closed core saturation) and the series of water saturation values calculated from the water spectrum function, the minimum average relative error is determined. This error is used as an index to derive the parameters of the water spectrum function. Finally, the obtained parameters of the water spectrum function are used to determine the final water spectrum function, which is then applied to other wells in the same block.

3. Simulation of Water Spectrum Function Parameters

Studying mercury injection and nuclear magnetic resonance core experiments (Figure 4), adopting the concept of the J-function, and considering the properties of fluid filling in rock pores, observe the distribution patterns of saturation of a certain wetting or non-wetting phase with pore radius. This can be understood as the pores corresponding to the saturation of a certain non-wetting or wetting phase in the rock [26]; the remaining pores are filled with either wetting or non-wetting phase fluids. It was found that transverse relaxation or capillary pressure exhibit certain distribution patterns with saturation. Through fitting with a weighting function (Figure 5), it was found that the fitting goodness was 0.997. According to the improved nuclear magnetic pore model, before 2.5 ms, the pores are completely filled with clay-bound water, while after 2.5 ms, larger transverse relaxation corresponds to larger-capillary radii, indicating more oil filling, and smaller transverse relaxation corresponds to smaller-capillary radii, indicating more water filling.

To simplify the water spectrum function and reduce the number of parameters that need to be determined, it can be inferred from the function fitting formula constructed in Figure 5 that the parameters of the water spectrum function can be simplified. Letting a1 = 1 and a2 = 0, the water spectrum function becomes:

$$\{\begin{array}{l}S(T_2)=1,T_2\le 2.5\\ S(T_2)={\displaystyle \frac{1}{1+{\left(\frac{T_2}{T_{2CW}}\right)}^m}}\end{array},T_2>2.5$$

(11)

This means that only the parameters m and T_2CW need to be determined.

Studying the effect of varying parameter m on the properties of the water spectrum function, as shown in Figure 6, when T_2CW is fixed at 33 ms, reveals the corresponding water spectrum functions for different values of m. It is observed that when m is positive, a larger m corresponds to a greater span of weights associated with the change from small to large transverse relaxation; smaller values of m correspond to weights closer to 50% for the transverse relaxation. Conversely, when m is negative, smaller values of m correspond to a greater span of weights associated with the change from small to large transverse relaxation, while larger values of m correspond to weights closer to 50% for the transverse relaxation.

Studying the effect of T_2CW variation on the properties of the water spectrum function, as shown in Figure 7, when m is fixed at 1, it depicts the corresponding water spectrum function as T_2CW changes. It is observed that as T_2CW increases, the oil–water line moves in the direction of increased transverse relaxation, indicating an increase in the T₂ value when oil and water each occupy 50% of the pore fluid.

The movable fluid cutoff value obtained from the core NMR experiment in the Satuo block is 23.53 ms (the cutoff value is obtained through NMR experiments, which include but are not limited to this parameter). In this study, the effect of changes in m and T_2CW on the weighting of the water spectrum function is investigated when the transverse relaxation is 23.53 ms. As shown in Figure 8, the relationship between the weighting distribution and m and T_2CW when the transverse relaxation is 23.53 ms is depicted. According to the properties of the water spectrum function, it is evident that when T_2CW is 23.53 ms, the weight remains at 0.5 regardless of the variation in m; as m increases, the span of the weight change corresponding to the range of T_2CW from small to large increases.

The empirical parameter for the clay-bound water cutoff value in the Satuo block is 2.5 ms, and the conventional model defaults to T₂ < 200 ms for water and T₂ > 200 ms for oil. Combined with the analysis in Figure 7, the optimization range for T_2CW is set to [2.5, 200]. Further, from the comprehensive analysis of simulation results, it is evident that excessively large values of |m| can lead to extreme weights at both ends of the oil–water line, as shown in Figure 6. Therefore, m is constrained to |m| < 10. Thus, considering the range of m within [−10, 10] and T_2CW within [2.5, 200], particle swarm optimization (PSO) is applied to optimize m and T_2CW using core saturation as the label.

Typically, the selection of the parameter T_2CW considers values between the clay-bound water cutoff and the oil–water line. This selection helps to avoid falling into the misconception of pure water or pure oil pore models [27]. As for the parameter m, the selection usually considers values within |m| < 10. From the simulation results in Figure 6, it is evident that regardless of whether the system is oil-wet or water-wet, as |m| approaches 10, the distribution span at both ends of the oil–water line increases. Excessively large values of m can lead to oil–water separation phenomena, which contradict common sense. Therefore, the selection range for parameter m is considered within |m| < 10.

4. Parameter Optimization Using Particle Swarm Algorithm (PSO)

The specific process for determining the parameters of the water spectrum function involves removing the clay-bound water portion (i.e., T₂ values after the effective cutoff value, the clay-bound water cutoff value) to obtain the remaining portion and constructing the water spectrum to calculate water saturation in the absence of known parameters m and T_2CW [28]. Using the core water saturation as the label, the parameters m and T_2CW are then determined by indexing the minimum average relative error. This is achieved by employing the particle swarm optimization (PSO) algorithm to compute the optimal values of m and T_2CW corresponding to the minimum average relative error.

Before optimizing the parameters of the water spectrum function, it is necessary to extract the nuclear magnetic resonance (NMR) logging T₂ distribution corresponding to the depth of the sealed core samples [29]. Here, the extraction is performed using nearest neighbors. The extracted T₂ distribution is then used to calculate the water spectrum saturation. With both parameters, m and T_2CW, unknown, they are optimized by calculating the average relative error between the water spectrum saturation and the sealed core saturation. A multi-objective optimization is conducted using the particle swarm optimization algorithm. When the iteration reaches the minimum average relative error, the corresponding m and T_2CW are considered the optimal solution.

Let the sealed core sample’s water saturation be S_w, and the water spectrum saturation calculated from the corresponding T2 distribution at depth be S_w’. The average relative error can be expressed as follows:

$$MRE={\displaystyle \frac{{\displaystyle \sum _i\left|S_{w_i}-s{w}_i\right|}}{n}}={\displaystyle \frac{{\displaystyle \sum _i\left|\frac{{\displaystyle \sum _{j}S(T_{2_{ij}})por(T_{2_{ij}})}}{{\displaystyle \sum _{j}por(T_{2_{ij}})}}+MCB{W}_i-s{w}_i\right|}}{n}}$$

(12)

where MCBW represents the clay-bound water saturation, n is the number of sealed core samples, S_wi refers to the water saturation calculated using the NMR water spectrum function, while s_wi represents the core-sealed saturation.

The particle swarm optimization algorithm initializes a group of random particles within the range of m [−10, 10] and T_2CW [2.5, 100], and then iteratively searches for the optimal solution for m and T_2CW. Each particle has only two attributes—velocity and position, where velocity represents the speed of movement and position indicates the direction of movement. Each particle independently searches for the optimal solution for m and T_2CW in the search space and records it as its current individual best solution. These individual best solutions are shared among all particles in the entire swarm, and the best one among them is selected as the current global best solution for the entire swarm. Based on their own current individual best solution and the current global best solution shared by the entire swarm, all particles adjust their velocity and position.

The pseudo code of the particle swarm optimization algorithm is shown in Algorithm 1, and its process can be summarized as follows [30]:

(1) Randomly initialize each particle.

(2) Evaluate the global optimal solution for each particle and terminate if the condition is met.

(3) Update the velocity and position of each particle and evaluate its fitness function.

(4) Update the historical best position of each particle.

(5) Update the global best position of the new group and terminate if the condition is met.

Algorithm 1 Pseudocode for particle swarm algorithm (PSO)

1: for each particle i

2:  Initialize velocity Vi and position Xi for particle i

3:  Evaluate particle i and set pBesti = Xi

4: end for

5: gBest = min (pBesti)

6: while not stop

7:  for i = I to N

8:   Update the velocity and position of particle i

9:   Evaluate particle i

10:    if fit (Xi) < fit (gBest)

11:    pBest = Xi

12:    end if

13:    if fit (pBesti) < fit (gBest)

14:     gBest = pBesti;

15:    end if

14:   end for

15: end while

16: print gBest

Particle swarm optimization was conducted to optimize the parameters of the water spectrum function using a total of 337 core samples from 8 wells (

Table 1. Samples for parameter optimization using particle swarm algorithm (PSO).

Well Name	Stratigraphic Horizon	Number of Labels
#1	P2	17
#2	P2	109
#3	P2	21
#4	P2	32
#5	P2	31
#6	P2	36
#7	P2	36
#8	P2	57

). Figure 9 shows the iterative process, with a final average relative error of 0.11, indicating good performance. Figure 10 illustrates the optimal positions of m and T_2CW corresponding to the average relative error. With T_2B already determined as 2.5 ms, the calculated values are m = 1.1602 and T_2CW = 24.2912. Figure 11 displays the obtained water spectrum function.

$$\{\begin{array}{l}S(T_2)=1,T_2\le 2.5\\ S(T_2)={\displaystyle \frac{1}{1+{\left(\frac{T_2}{24.2912}\right)}^{1.1602}}}\end{array},T_2>2.5$$

(13)

The purpose of optimizing the parameters of the particle swarm optimization algorithm for the water spectrum function is to improve the iteration speed and accuracy in order to obtain the optimal values of m and T_2CW. This enhances the computational time and efficiency. By optimizing the parameters using the particle swarm optimization algorithm, one only needs to provide the optimization range, thereby obtaining the optimal water spectrum function parameters suitable for the specific block [31].

The water spectrum function parameters obtained through this method, m and T_2CW, are fixed values that can be understood as empirical parameters suitable for the S_a exploration block. If saturation calculations for the water spectrum function are to be conducted for different blocks, it is necessary to first determine the water spectrum function parameters for that block through core experiments and measured nuclear magnetic resonance logging T₂ distributions, and then to perform saturation calculations for the water spectrum function for that block.

5. Case Study Verification and Application of Water Saturation Calculation

The parameters obtained from the particle swarm optimization algorithm, namely, m = 1.1602 and T_2CW = 24.2912, were combined with the improved nuclear magnetic resonance (NMR) porosity model. Together with the clay-bound water cutoff value of 2.5 ms, the water saturation was calculated using Equations (5) and (9) from the water spectrum function shown in Figure 11.

Figure 12 illustrates the application example in Well S1, where downhole measurements were conducted using the MRILP-type NMR tool. The water saturation was calculated using the water spectrum function from the 9th channel joint spectrum, and the comparison between the calculated saturation and the core saturation is presented in the 13th channel. The calculation results demonstrate good performance.

Figure 13 depicts the application example of Well S2, where downhole measurements were conducted using two-dimensional nuclear magnetic resonance (NMR). The water saturation was calculated using the water spectrum function from the T₂ distribution of the 10th channel, and the saturation was compared with the core saturation in the 12th channel, showing significant computational efficacy. Additionally, perforation was performed at the core sampling location of this well, and the well testing result indicated a daily oil production of 95.99 tons and a daily water production of 13.09 cubic meters. The water saturation calculated using the water spectrum function that ranged from 70% to 90%, which closely matched the actual situation.

As shown in Figure 14 and Figure 15, the error analysis for wells S1 and S2 indicates that the regression analysis R² is greater than 0.4, and the average relative error is less than 5%. This demonstrates that the saturation calculated from nuclear magnetic resonance logging has a certain degree of accuracy.

The validation through examples shows that the constructed water spectrum function, along with the improved porosity model and saturation calculation method based on the T₂ distribution from nuclear magnetic resonance logging, is suitable for the oil-bearing analysis of tight sandstone reservoirs. This porosity model and saturation calculation method not only enhance the accuracy of oil-bearing evaluations but also provide more effective and precise tools for oilfield development and assessment under complex geological conditions.

6. Conclusions

This paper proposes a new NMR porosity model and presents a method for calculating saturation from the T₂ distribution in NMR logging using a water spectrum function. By focusing on the fluid distribution characteristics in different pore structures of tight sandstone reservoirs, this method overcomes the challenge of the strong dependence on resistivity for oil saturation calculation, enabling accurate saturation calculations in tight reservoirs. It takes into account the coexistence of oil and water with different transverse relaxation times and identifies an NMR-based saturation calculation method suitable for the S_a Block. The water spectrum function is adaptable to different blocks, improving the accuracy of saturation calculations in tight reservoir exploration and development, thereby providing strong support for more precise reservoir production forecasts:

(1) The traditional NMR porosity model was improved to avoid the limitations of fixed transverse relaxation intervals for bound water, movable water, and movable hydrocarbons.

(2) A water spectrum function was constructed, assigning weight ratios for oil–water coexistence in pores outside the clay-bound water zone. Combined with the improved NMR porosity model, this approach enables precise saturation calculation.

(3) The particle swarm optimization (PSO) algorithm was employed to optimize the parameters of the water spectrum function. By establishing an objective function based on the average relative error of core saturation, the water spectrum function parameters were iterated quickly.