Permeability Calculation of Sand Conglomerate Reservoirs Based on Nuclear Magnetic Resonance (NMR)

PetroChina Research Institute of Petroleum Exploration and Development, Beijing 100083, China School of Civil Engineering, Hebei University of Engineering, Handan 056000, China School of Energy Resource, China University of Geosciences (Beijing), Beijing 100083, China Beijing Key Laboratory of Unconventional Natural Gas Geological Evaluation and Development Engineering, Beijing 100083, China Earth Science Faculty, China University of Geosciences (Wuhan), Wuhan 430000, China


Introduction
The pore structure of reservoir rocks has fractal characteristics [1]. Because of this, a number of scholars have recently adopted fractal theory to establish the porosity or the permeability model of reservoir rocks and have achieved some productive results. For instance, Li et al. [2], Zheng et al. [3], and Bai et al. [4] have successfully applied fractal theory to build models for forecasting the permeability of tight sandstone reservoirs. With the development of fractal theory, however, it has become apparent that the theory has limitations in characterizing a real reservoir space. When computing fractal dimensions using capillary data, for example, it has been noted that in a log-log coordinate system, the lg S Hg to lg P c plot sometimes exhibits multifractal characteristics for some kinds of rocks, instead of being linearly fitted, which significantly influences the accuracy of computation [5][6][7][8][9]. Meanwhile, Sima et al. [10] have pointed out that sand conglomerate rocks with high heterogeneity tend to exhibit more evident multifractal characteristics or even nonfractal characteristics, on log-log coordinates, illustrating that the probability of multifractal characteristics is greater in the conglomerate samples. This section analyzes the fractal feature gained from the capillary buddle subsection model and draws similar conclusions for sand conglomerate rocks ( Figure 1). Figures 1(a)-1(c) show that the fractal characteristics of the capillary pressure curve of the sand conglomerate reservoir generally need to be fitted with two or more straight lines, which indicates that different intervals of pore size have different fractal dimensions. Figure 1(d) shows that the fractal dimensions of some large intervals of pore size, calculated by the graphical method, are greater than 3, and this phenomenon is interpreted as being attributable to nonfractal characteristics.
To address this problem , Li et al. [11] and Lai et al. [12] proposed to run a weighted average of the mercury saturation of each fractal segment to obtain a weighted fractal dimension. However, this method neglects the detailed features of the pore space and hence only reflects the overall fractal characteristics of rock pores. Moreover, although the fractal dimension is a highly significant parameter in fractal theory, it has limitations in characterizing a real rock pore space. This is because porous media with different pore structures may have the same or similar fractal dimension [13]. Figure 2 shows that although the pore structures of the four samples are different, their fractal dimensions tend to be similar. The specific calculation process for threedimensional fractal dimensions is set out by Zhang et al. [15]. Therefore, utilizing a single fractal dimension to establish a porosity or permeability model might inevitably constrain subsequent research on permeation performance and reservoir engineering.
Given the limitations of the single fractal dimension, Italian scholar Pia et al. [16][17][18][19][20][21][22][23] introduced an intermingled fractal unit (IFU) model, which has proven to be useful in simulating porous media. Pia and his research team applied the model to predict the thermal conductivity, electrical conductivity, mechanical properties, and permeability of porous media materials. Intuitively, the IFU model can be considered to be several superimposed fractal units. Each of the small fractal units in the IFU model is a modified Sierpiński carpet, with each fractal unit having different iteration parameters. Thus, the IFU model can be regarded as the Sierpiński Carpet modified by altering its iteration rule. The advantage of this model is that it can present detailed characterizations of internal spaces of complicated objects.
At present, applications of the IFU model have relatively little engagement in the field of petroleum engineering. Since Pia first used the IFU model to predict permeability based on mercury intrusion capillary pressure (MICP) data [20], some scholars have also done a number of works to widen the application of the model in recent years. Tan et al. [24] and Chen et al. [25] applied the IFU model to forecast the permeability with consideration of tortuosity according to the MICP data. Tan et al. [26] used X-ray computed tomography (CT) and scanning electron microscopy (SEM) data to make permeability prediction based on the model. Li et al. [27,28] applied a 2D and 3D IFU model to predict permeability using SEM images for shales. Objectively, at present, the applications of the IFU model in reservoir physics basically focus on the prediction of permeability, and the IFU model based on NMR data has not been reported.
This paper, based on the Pia IFU model, constructs a model to fit the pore size distribution spectrum according to NMR measurements and makes permeability predictions using the free fluid T 2 spectrum. Meanwhile, some key issues in the process of modeling using NMR data are also pointed out. Also, new intermingled fractal expressions of the Kozeny-Carman equation and relative permeability are derived to broaden the application of the IFU model in reservoir physics.   [20]. Without derogation from general principles, a specific IFU is used as an example to examine the relationships between its iteration parameters ( Figure 3).

2.5
Let d max be the maximum pore side length of this fractal unit. The side length of the fractal unit will then be 3d max . N p is the number of the child squares removed from each parent square in every iteration, and b is the ratio factor (a multiple of equal divisions of the side length in each iteration). Thus, N p = 1 and b = 3. N soild is the number of the child squares that do not take part in the iteration. Thus, N soild = 0. The value i is the iteration time.
According to fractal geometric theory [29], the fractal dimension ðD f Þ is the value of the number of squares that participate in the next iteration ðN iteration Þ divided by the ratio factor ðbÞ: For the Sierpiński carpet fractal model, the fractal dimension of the 3D spatial model is the value of the 2D plane model plus 1.

Geofluids
To produce a Sierpiński carpet similar to the unit shown in Figure 3, the side length of the new pores formed in the i-th iteration is defined by Pia et al. [20] as follows: The number of new pores formed in the i-th iteration ðN i Þ is given by Li et al. [27]: For ease of theoretical calculation, it is assumed that the real percolation section of every square pore is its maximum incircle and that the diameter, position, and quantity of pores will be the same on any rock section when the 2D IFU model is extended to 3D. If L 0 is the length of the model and the total iteration time for the fractal unit is k, then the generated total pore volume V p is given as: The IFU model can then be utilized to reproduce the pore size distribution spectrum and porosity. The procedure used to reproduce porous microstructures is shown in Figure 4. For a detailed explanation of the process, see Pia's series of papers.

Permeability Calculation Using Traditional Classic
Models Based on NMR Experiments. NMR experiments can quantitatively characterize the pore structures of rocks and can be applied directly to calculate permeability. The current classic models for calculating permeability using NMR data are the Timur-Coates and SDR models.
The Timur-Coates model [30] is expressed as follows: The SDR model [31] can be expressed as follows: where K is the permeability (10 -15 m 2 ), ϕ is the effective porosity (%), T 2g is the geometric mean of T 2 (ms), FFI is the saturation of movable fluid (%), BVI is the saturation of immovable fluid (%), and m 1 , m 2 , m 3 , n 1 , n 2 , and n 3 are the corresponding parameters of each model, which are determined by the results of actual core data used to fit experimental permeability. In the Coates model, in order to facilitate the calculation, the values of n 2 and n 3 are sometimes directly regarded as 4 and 2, respectively, as well as the values of m 2 and m 3 in the SRD model [32,33].
In recent years, some scholars have pointed out that the classic Timur-Coates and SDR models have certain limitations for calculating the permeability of complex reservoirs. For example, the classic model has good applicability for the calculation of permeability in medium and high permeability reservoirs. However, for reservoirs with low porosity and low permeability, it is difficult to correlate the calculation results with experimental data [34]. This paper applied the classical model to predict the permeability of sand conglomerate reservoirs in the Mahu rim region of the Junggar Basin, and the calculation results exhibit significant discrepancies from the experimental results (see Section 3.2). This is because the sand conglomerate rocks have extremely high heterogeneity due to their complex gravel-cement configuration. By using micro-CT scanning technology on some typical sand conglomerate rock samples, it is indicated that the sand conglomerate rocks in the research area have complex and diverse microscopic pore structures ( Figure 5).
Meanwhile, the actual calculations in the classic models are complicated, resulting in large volumes of calculations for the process of converting the NMR T 2 spectrum into a pseudocapillary pressure curve using the nonlinear method. Moreover, it is found that the conversion coefficients and conversion indices calculated by the classic model vary greatly, even for rock samples in the same study area. The average value of the region does not satisfy the requirement for calculation accuracy, and the nonlinear fitting method does not explain the theoretical relationship between T 2 spectral relaxation time and real pore sizes. Because of this, this paper utilizes the "centrifugation T 2c method" [35] to convert the T 2 spectrum into a corresponding pore size distribution for subsequent permeability prediction based on the IFU model.
There is an issue worth pointing out here. Since the sand conglomerates are easy to split in the process of centrifugal experiment, the cores are wrapped with a special heatshrinkable film in the experiments. Meanwhile, in order to determine the centrifugal force suitable for the samples in the study area, five representative samples were selected to test six different centrifugal forces of 100 psi, 150 psi, 200 psi, 250 psi, 300 psi, and 350 psi. By analyzing the descend range of irreducible water saturation before and after centrifugation, it was found that until the centrifugal force 4 Geofluids increased from 300 psi to 350 psi, the reduction of irreducible water saturation of five samples with different physical properties tended to be stable. Therefore, we select a centrifugal force of 300 psi to calibrate a T 2 cutoff value, and then, the corresponding lower limit of the movable fluid throat radius is 0.0704 μm. The experiment adopts a GeoSpec2 NMR core analyzer produced by Oxford Instruments Magnetic Resonance in England. The working frequency is 2 MHz, the temperature of the instrument is 35°C, the room temperature is 25°C, and the humidity is 50%-60%.

Permeability Calculation
Using the IFU Model Based on NMR Experiments. The expression of permeability based on the IFU model is derived in this section. The flow of a single pore formed after the i-th iteration is given by the Hagen-Poiseuille equation: where q is the flow (cm 3 /s), Δp is the pressure difference between the ends of the capillary tube (atm), μ is the viscosity of fluid (mPa·s), L 0 is the characteristic length of the capillary tube (cm), and τ is the tortuosity (dimensionless parameter). The ratio of the capillary tube's real length ðL t Þ to the characteristic length ðL 0 Þ is expressed: Based on the average tortuosity-porosity analytical model for porous media consisting of square particles introduced by Yu et al. [36], Yun et al. [37] proposed an enhanced expres-sion of the average tortuosity of porous media consisting of circular particles, expressed as follows: where ϕ is the porosity. Thus, the total flow (Q) of the pores formed from the first to the k-th iteration is given as: According to the equivalent percolation principle, Finally, the expression for computing permeability can be derived as follows: Thus, the permeability of a single fractal unit is given in formula (12), and the formula is similar in form to the   [20]. If the IFU model contains several fractal units, formulas (9) and (10) can be used to calculate its tortuosity and flow rates, respectively. The total permeability of the IFU model is easily obtained according to the equivalent seepage principle. 6 Geofluids determined by the centrifugal experiment, is a constant, and thence, the pore size r i , corresponding to the i-th relaxation time T 2i , is given as [38]

Results and Discussion
During the NMR experiment, the centrifugal pressure for testing is 300 psi, so the calculated pore size r 0 is 0.0704 μm, corresponding to T 2c . By applying formula (13), the pore size corresponding to different relaxation times can be calculated, and the pore size distribution based on T 2 spectra can then be constructed as shown in Figure 6.
The T 2 spectrum of 100% saturated water is the reflection of all pore spaces in the rocks [39]; i.e., it contains connected pores, dead pores, and micropores storing bound water. Therefore, for permeability calculations based on NMR data, the free fluid T 2 spectrum should be used rather than the 100% saturated water T 2 spectrum. This is because only connected pores with percolation ability contribute to permeability in the pore space, and thus, only the free fluid T 2 spectrum has good correspondence with the real seepage space. The acquisition method for this value is to subtract the centrifugal T 2 spectrum signal from the saturated water T 2 spectrum sig-nal, and the remainder is the free fluid T 2 spectrum ( Figure 7). Also, an abnormality persists for some samples, in that the centrifugal T 2 spectrum in some samples is higher than the saturation T 2 spectrum (Figures 8(a) and 8(b)), resulting in the appearance of a "negative accumulation area" when the free fluid accumulation curve is obtained (Figures 8(c) and  8(d)). These anomalies in the micropore intervals are difficult to avoid, because they are caused by system errors in the instrumentation. Therefore, when using the IFU model to fit the pore size distribution spectrum of the free fluid, it is suggested that the model is fitted to the highest point of the cumulative curve to circumvent system errors.
To summarize, the procedures of conversion of the T 2 spectrum to pore size distribution are as follows: (1) Subtract the centrifugal T 2 spectrum from the 100% saturated water T 2 spectrum to obtain the free fluid T 2 spectrum (2) Eliminate the "negative accumulation area" to obtain an effective free fluid accumulation curve (3) Use formula (13) to convert the effective free fluid T 2 spectrum to the pore diameter cumulative distribution curve (1) Select the fractal units with different parameters for building the IFU model (Table 1), and change the iteration parameters to assure a consistency between the pore diameter distribution spectra of the structural model and the experimental data. Figure 9 reveals fitting results for some representative samples (2) Based on the completion of the above fitting process, formula (12) can be applied in the IFU model to calculate permeability (Table 2)

Geofluids
In order to make a better comparison between the computed results of the IFU model and the classical models, we subdivide the two classical models into the following four types according to whether the undetermined parameters of the models are constant or not. The relevant parameters of the experiment and models are shown in Tables 3 and 4.
The parameters in Table 3 are substituted into each model in Table 4, and the best model coefficient is determined by programming using the method of three-coefficient regression analysis. It can be seen that in Coates model 1 and SDR model 1, except that the coefficient is a variable, the index is also not a constant. Actually, when the regression statistical method is used to determine the coefficients and exponents in the equation, due to the lack of data points, the variation range of model parameters is closely related to the selection of initial values, which requires some prior knowledge. Therefore, a more feasible method is to take logarithms on both sides of the equation and then calculate each undetermined coefficient by regression statistical analysis. Table 5 and Figure 10 show that the classical models cannot accurately calculate the permeability of the sand conglomerate reservoirs in the study area. The average relative error of the results calculated using the IFU model is 23.1%, within a certain tolerance for error, indicating that the IFU model is reliable in simulating the movable fluid pore space based on NMR data.
Also, it is noted that the permeability calculation results, based on the IFU model, are mostly smaller than the experimental ones. This is mainly due to the fact that, for reservoirs with complex pore structures such as sand conglomerate reservoirs, it is difficult to replicate the real conditions of a bound water state in centrifugal experiments. Centrifugal parameter settings often have strong regional characteristics. However, in fact, in order to truly replicate a bound water state and obtain the free fluid T 2 spectrum, the centrifugal force and centrifugation time must be markedly different for each sample.
Therefore, it is not realistic to increase the centrifugal force and the centrifugation time without reference to other parameters, because "grain shattering" occurs easily during the centrifugation procedure for sand conglomerate rocks. Excessive centrifugal force will damage the samples and then affect the experimental results. Therefore, since it is difficult to accurately measure the content of the movable fluid in centrifugal experiments, the computational accuracy of the IFU model based on NMR data is inevitably somewhat compromised.

Intermingled Fractal Expression of the Kozeny-Carman Equation.
The Kozeny-Carman equation is a semiempirical formula widely used in calculating permeability, whereas the determination of the KC constant is highly empirical. In fact, the KC constant is closely related to the microscopic pore structure and is not a fixed value. Since the IFU model can accurately quantify specific pore throat information in the porous medium, the KC constant, with a practical physical manifestation, can be obtained using the IFU.

Geofluids
In order to facilitate the presentation of the formula derivation process, this section takes the example of an IFU model with only one fractal unit. Taking tortuosity into consideration, the total internal surface area (A in ) of the pores produced after the k-th iteration is straightforwardly obtained according to equation (3): Thus, the specific surface (S) of the IFU model is From equation (4), the porosity (ϕ) of the IFU model can be rendered as where A is the cross-sectional area of the model. Note that for some tight reservoirs with extralow porosity and permeabil-ity, the pore sizes converted using the NMR free fluid T 2 spectrum are the effective flow radius, without considering boundary layer effects. However, for pore sizes converted by mercury injection capillary pressure (MICP) measurements, it is necessary to deduct the thickness of the boundary layer to obtain an effective flow radius, since mercury does not produce a boundary layer in MICP experiments. For quantitative calculation of the thickness of the boundary layer, please refer to the references made by Tian et al. [40] and Meng et al. [41]. Thus, permeability of the IFU model, taking the thickness of the boundary layer (h i ) into consideration, is given by According to the classical Kozeny-Carman equation, where C is the KC constant. By substituting formulas (15)- (17) in (18), the expression of the KC constant is obtained: It can be seen that the KC constant is not an empirical constant without precise physical meaning but is controlled by factors such as iteration parameters, tortuosity, and the thickness of the boundary layer in the IFU model. Once these parameters are given, the value of KC for each sample can be calculated simply by programming.

Intermingled Fractal Expression of Relative Permeability.
Based on the IFU model, the expression of relative permeability can also be derived. Similarly to Section 3.3.1, the IFU model with only one fractal unit is taken as an example in this case. The assumptions for the model are as follows: (a) The reservoir is composed of unequal diameter capillaries with fractal features (b) The fluid flow conforms to Darcy's law, and the wet phase is distributed in the capillaries whose pore diameters are smaller than the critical capillary radius (r c ) and the nonwet phase is distributed in the capillaries whose pore diameters are larger than the critical capillary radius. Therefore, only one fluid exists in a single capillary (c) The boundary layer fluid adheres to the inside of the tube wall in the form of a bound water film (d) The fluid is incompressible, the viscosity is constant, and the percolation process is isothermal Thus, considering the thickness of the boundary layer, formula (12) can be rewritten as According to the definition of saturation, 10 Geofluids

Geofluids
The effective permeability can then be rewritten as It is indicated from the previous study that, as the iteration time (i) increases, the diameters of the generated pores reduce. Therefore, if c is the critical iteration time, then the nonwet phase exists in the pores of i < c and the wet phase exists in the pores of i ≥ c.
The effective permeability of the wet phase is then given by The effective permeability of the nonwet phase is given by The relative permeability of the wet phase can then be derived as follows: And the relative permeability of the nonwet phase is where τ, τ w , and τ nw are the tortuosities of 100% saturated single-phase, wet-phase, and non-wet-phase fluids, respectively, and they are determined by formula (9) according to their respective saturations. Also, it is found that the forms