A Novel Model for Predicting Permeability Using Porosity Frequency Spectrum in Fractured Deep Metamorphic Rock Reservoirs: A Case Study of the Archaeozoic Formation in BZ 19-6 Region of Southwestern Bohai Bay Basin

1. Introduction

Permeability played an important role in indicating formation productivity [1]. Generally, reservoirs with high permeability always had better filtration capacity and higher hydrocarbon production. On the contrary, if a formation contained ultra-low permeability, the effectiveness would be worst, and it would be low-yield or dry layer. Precisely predicting permeability was indispensable in reservoir evaluation and effectiveness characterization, especially in ultra-low to tight reservoirs [2,3].

In generation, total porosity was always used as a single input parameter to predict permeability after power function or exponential function was applied by using plenty of routine core analysis data [4,5]. This method was available only in formations with relatively simple pore structure, and the relationship between total porosity and permeability could be easily expressed. Permeability prediction faced great challenge in formations with complicated pore structure, such as low-permeability to tight sandstone or carbonate reservoirs [6]. Factors affecting formation permeability in such types of reservoirs were multiple. Besides total porosity, the ratio of movable pores, the pore morphology, the irreducible water saturation and so on, also heavily affected formation permeability. Hence, many models, e.g., the Timur -Coates model, the Schlumberger Doll Research Center (SDR) based model, the Swanson based model, and the R10 (defined by the pore-throat radius corresponding to 10% mercury injection saturation) based model, were proposed to involve these parameters and used to predict permeability [7,8,9,10,11,12,13]. In addition, the hydraulic flow unit (HFU) approach was introduced to classify formations into several clusters, and relationships between total porosity and permeability were established to extract permeability from total porosity for every cluster of formation [14,15]. These methods played great important roles in estimating permeability in formations with relatively simple pore structure. In such type of formation, the pore spaces were dominated by intergranular pores, and the main flow paths were pores and throats, and the throat-to-pore coordination number was positive to permeability [16]. However, in reservoirs with secondary pores, e.g., carbonate, volcanic rock or metamorphic rock, all the current methods could not operate. This is because there were two to three types of pore spaces in such types of reservoirs [17,18,19,20]. Apart from the intergranular pores, secondary pores, such as fractures or caves, were widely developed in these types of formations. Although the fracture porosities were not high, they could effectively connect the intergranular pores, thereby enhancing rock permeability. To effectively predict permeability in fractured reservoirs, the contribution of secondary pores should be considered.

For permeability evaluation in fractured formations, various types of models have been proposed based on electrical imaging logging [21,22,23,24,25,26]. The most classic one is the slab model [27]. In the slab model, formations were cracked and classified into several parallel clusters by fractures, and relationship between fracture aperture and permeability were established based on exponential function. Formation permeabilities were overestimated in many fractured reservoirs based on the slab model, the reason was that the fracture surface was considered smooth and parallel. The real form of fracture surface was complex and irregular, this would reduce rock permeability. In addition, the calculation of fracture aperture from electrical imaging logging was a worldwide challenge, all these limited the widely applicability of slab model. Li et al. (2024) raised a model to predict permeability from porosity frequence spectrum based on the morphological method [23]. In this model, the porosity frequency spectrum was classified into two parts based on the shape, and a normal distribution function was used to simulate the distribution of matrix pore. The ratio of the area of the distribution of matrix pore to total area of porosity frequency spectrum was defined as an input parameter in this model. This model had been verified to be valuable in fractured tight sandstone reservoirs in Ordos Basin. However, the operability of this model was questionable. Because fractures were random and widely developed in metamorphic rock reservoirs, especially the reticular fractures. This made the sharp of porosity frequency spectra was multimodal, it was difficult to extract the distribution of matrix pore from total porosity frequency spectrum. Xiao et al. (2000) pointed out that the porosity frequency spectrum could be first translated as pseudo capillary pressure (Pc) curve to characterize formation pore structure, and then the Swanson based model or R10 based model were applied to calculate permeability [3]. The applicability of this method was determined by the accuracy of extracting pseudo-Pc curve from porosity frequency spectrum. This resulted in this model being effective only in formations where intergranular pores were developed. The array acoustic logging was considered to be effective in indicating the development of fracture, and the corresponding permeability prediction model based on the Stoneley wave was raised and widely used all over the world in the last several decades [28]. This model was usable in oil-bearing or water saturated layers. However, the effect of pore fluids could not be removed in gas-bearing reservoirs, this made that the predicted permeability was underestimated. Comparing with nuclear magnetic resonance (NMR0 logging, the electrical imaging logging was considered to have advantages in reservoir characterization in the following three aspects: (1) the extracted parameters did not be affected by the pore fluids; (2) the high longitudinal resolution enabled the effective identification of thin layers; (3) the effectiveness was high in fractured reservoir evaluation [29].

The purpose of this study was to raise a novel model to directly predict permeability from porosity frequency spectrum in deep metamorphic rock reservoirs after the shapes of porosity frequency spectra in various types of formations were compared. Two parameters of the logarithmic mean and standard deviation between two golden section points were extracted, and relationships among these two parameters and permeability were established. After this raised model was applied in field applications in the Archaeozoic Formation in BZ 19-6 Region of southwestern Bohai Bay Basin, much reasonable permeability was predicted.

2. Geological Setting

The Bohai Bay Basin was a superimposed basin developed on the ancient Paleozoic cratonic basement of North China, formed through strike-slip and rifting processes [30]. The BZ Sag was a primary structural unit of the Bohai Bay Basin and represented the largest hydrocarbon-bearing sag in the Bohai Bay Basin. BZ 19-6 Region was located in southwestern Bohai Bay Basin, characterized by an anticline structure formed between the southwestern sub-sag and the southern sub-sag. The geological reserves of natural gas in BZ 19-6 Region had been verified to exceed 1.0×1011 m3, making it the largest gas field discovered in eastern China [31]. The BZ Sag was located on the basement rocks, with complete absence of Paleozoic and Proterozoic strata, and partial absence of Mesozoic strata. From top to bottom, formations encountered during drilling contain the Quaternary Pingyuan Zone, the Neogene Minghuazhen Zone and Guantao Zone, the Paleogene Dongying Zone, Shahejie Zone, and Kongdian Zone, as well as the Mesozoic Zone and the Archean metamorphic basement (Figure 1). Our target formation was the Archaeozoic metamorphic rock formation. The average burial depth was higher than 4400 m，and the formation temperature was higher than 140 ℃, belonging to typical deep metamorphic rock reservoir. The main lithology contained gneiss, ganite, cataclasite and granodiorite. The pore types could be classified into three major categories, namely fracture type, fracture-pore type and pore-fracture type. Complicated lithology and pore types made formation evaluation facing great challenge in the Archaeozoic metamorphic rock reservoirs. In many regions, formations with low porosity but high hydrocarbon production were frequently found. The relationship between core-derived porosities and permeabilities was unusual. As the porosities increased, permeabilities did not significantly rise (Figure 2). It was urgent to predict permeability in fractured formations.

3. Analysis of Morphological Characteristics of Porosity Frequency Spectrum

3.1. Extraction of Porosity Frequency Spectrum from Electrical Imaging Logging

The electrical imaging logging had plenty of button electrodes evenly distributed within a 360-degree range around the wellbore (e.g., formation micro-resistivity scanning image logging (FMI) tool had 192 electrodes, star II tool had 144, whereas the number of electrodes for electrical resistivity micro imaging (ERMI) tool was 150). These button electrodes divided formation (pores) into several small units. For each unit, its conductivity followed the inter-granular flow law. Meanwhile, electrical imaging logging had a relatively shallow invasion radius, the measured results only reflected the conductivity of flushed zone [32]. Generally, the pore fluids were mud filtrate and small amount of negligible residual hydrocarbons in flush zone. Hence, the Archie’s equation, which was raised based on the inter-granular flow law, could be used to reflect the relationships of conductivity, porosity and water saturation in flushed zone [33]. The Archie’s equation was transformed and expressed as follow:

$$\phi =\sqrt[m]{{\displaystyle \frac{a\times b\times R_{mf}}{S_{xo}^n\times R_{xo}}}}$$

(1)

Where, φ was the formation total porosity in v/v. R_mf and R_xo were the resistivity of mud filtrate and flushed zone, separately. The unit of them was Ω.m. S_xo was the water saturation in flushed zone in v/v. a, b, m and n were the constants related to conductivity, and they were named as the rock resistivity parameters.

If R_xo was replaced by conductivity, based on the theoretical relationship between resistivity and conductivity (eq. 2), we could derive a transformed equation from eq. 1:

$$C_{xo}={\displaystyle \frac{1000}{R_{xo}}}$$

(2)

$$\phi =\sqrt[m]{{\displaystyle \frac{a\times b\times R_{mf}}{{1000\times S}_{xo}^n}}\times C_{xo}}$$

(3)

Where, C_xo was the measured conductivity in flushed zone in S/m.

Combining with eq.s 2 and 3, and we inserted the measured conductivities from button electrodes, a derived formulate was:

$$\begin{gathered} \phi \left(i\right)=\sqrt[m]{{\displaystyle \frac{a\times b\times R_{mf}}{{1000\times S}_{xo}^n\times R_{xo}}}\times R_{xo}\times C_{xo}\left(i\right)} \\ =\phi \times \sqrt[m]{{\displaystyle \frac{R_{xo}\times C_{xo}\left(i\right)}{1000}}} \\ =\phi \times g\times \sqrt[m]{R_{xo}\times C_{xo}\left(i\right)} \end{gathered}$$

(4)

Where, φ(i) was the calculated i^th imaging porosity in v/v. C_xo(i) was the extracted conductivity from i^th button electrode in S/m. g was a constant, and its value was determined by m.

Eq. 4 illustrated that if our target formation was homogeneity, the measured conductivities C_xo (i) within a 360-degree range around the wellbore would be a fixed value, and the values of φ(i) would be equal to φ. However, true formation was strong heterogeneity, especially for low permeability to tight reservoirs, the measured C_xo (i) were various. This caused the values of φ(i) were not fixed, but fluctuated around total porosity. If we conducted a statistical analysis of the calculated φ(i) with given range of porosities, we would get the porosity frequency spectrum [34,35]. The x-axis of the porosity frequency spectrum would be the given range of porosities, whereas the y-axis was the statistical distribution frequency. The principle of obtaining porosity frequency spectrum from electrical imaging logging was exhibited in Figure 3.

3.2. Analysis of the Morphological Characteristics of Porosity Frequency Spectrum

As was expressed in eq. 4, the morphological characteristics of porosity frequency spectrum was determined by the measured conductivities C_xo(i) from electrical imaging logging. For formations with simple pore shape, e.g., intergranular pore or fracture, the measured conductivities were fixed or fluctuated within a very narrow range, and the extracted porosity frequency spectrum would exhibit as unimodality. However, once the pore type and the configuration relationship of pore and fracture were complicated, e.g., fracture-pore type or pore-fracture type, the rock conductivity would be entirely different. Because the mud filtrate easily invaded into the fractured pore space to improve the conductivity, this made the measured conductivity for every button electrode varied in a wide range, and the extracted porosity frequency spectrum would exhibit as bimodality or multimodality.

The cast thin sections of core sample illustrated that the pore spaces of our target Archaeozoic Formation in BZ 19-6 Region, southwestern Bohai Bay Basin, contained three clusters, which were single fracture type, fracture-pore type and pore-fracture type (Figure 4). In Figure 5a to c, we exhibited the pore frequency spectra of three types of formations. For every type of formations exhibited in this figure, they contained similar porosity but different lithology. Comparisons of the morphologies of porosity frequency spectra for different types of formations illustrated following characteristics: (1) Formations with single pore type (fractured formations with ultra-low porosity) exhibited as unimodal distributions in the porosity frequency spectra, even they contained various lithology. For example, four porosity frequency spectra displayed in Figure 5a, the corresponding formations covered four types of lithologies, but their morphologies were similar. The reason was that they contained similar porosities. (2) Once formation pore types were diverse, the shapes of porosity frequency spectra were absolute different. They exhibited as bimodal or multimodal distributions, and the larger the matrix porosity, the wider the porosity frequency spectrum (Figure 5b and c). (3) The shapes of porosity frequency spectra were solely determined by the pore type, but were independent of lithology. Like the respectively exhibited results in Figure 5a, b and c, formations with similar porosities but four different types of lithologies had similar porosity frequency spectra.

4. Establishment of Permeability Prediction Model Based on Porosity Frequency Spectrum

4.1. Extraction of Input Parameters from Porosity Frequency Spectrum

Figure 5a to c notably illustrated that the shapes of porosity frequency spectra were heavily associated with the pore type and pore-fracture configuration relationship. Formations with diverse pore types, such as fracture-pore type or pore-fracture type, often had better pore connectivity, more superior pore structure, and correspondingly better permeability, and vice versa. Hence, the shapes of porosity frequency spectra can be used to reflect rock permeability. The wider the porosity frequency spectra, the higher the permeability.

In order to quantify the correlation between porosity frequency spectrum and permeability, we extracted two parameters from porosity frequency spectrum to characterize the form, and defined them as the logarithmic mean and standard deviation between two golden section points (GSP). The physical significance and calculation method of these two parameters for different types of porosity frequency spectra were exhibited in Figure 6 and expressed in eq.s 5 and 6, separately.

$${\phi }_{gm}={10}^{{\displaystyle \frac{\sum _{i=1}^N\left({log}_{10}\right(\phi \left(i\right))\times amp(i\left)\right)}{\sum _{i=1}^N\left(amp\right(i\left)\right)}}}$$

(5)

$${\phi }_{gsr}={\phi }_r-{\phi }_l$$

(6)

where, φ_gm was the logarithmic mean of porosity frequency spectrum in %. Its value represented the position of the main peak of the porosity frequency spectrum. For example, formation had a more superior pore structure, and the position of the porosity frequency spectrum was relatively to the right, with the corresponding logarithmic mean value being larger. φ_gsr was the standard deviation between two golden section points in %. Its value represented the width of porosity frequency spectrum. The more complex the pore structure of the formation was, the wider the distribution of the porosity frequency spectrum would be, and the larger the corresponding φ_gsr would be.

After electrical imaging logging was processed, the porosity frequency spectra could be consecutive acquired, and φ_gm and φ_gsr were calculated. To analyze the relationship among these two parameters with permeability, we extracted them corresponding to the same depth with core samples, and crossed them with permeability, as were shown in Figure 7.

Comparison of Figure 2 and Figure 7 clearly indicated that φ_gm and φ_gsr had a much stronger correlation with permeability, and both of them were all positively correlated. As the φ_gm and φ_gsr increased, the corresponding permeability also raised. They were two sensitive parameters in indicating the change of permeability.

4.2. Establishment of Permeability Prediction Model in Fractured Metamorphic Rock Reservoirs

Figure 7 had clearly exhibited the positive correlation among φ_gm, φ_gsr and permeability. Hence, we choose φ_gm and φ_gsr as the input parameters to establish a permeability prediction model, and this model was expressed as:

$$K=C\times {\phi }_{gm}^p\times {\phi }_{gsr}^q$$

(7)

Where, K was the rock permeability in mD. C, p and q were the involved model constants, and their values needed to be first calibrated by using the experimented results of core samples.

It needed to noted that, comparing with porosity, φ_gm and φ_gsr were much positively associated with permeability, but the correlation that displayed in Figure 7 was not high enough. The reason was that the displayed core samples in Figure 7 were all collected from different types of reservoirs in the Archaeozoic Formation in BZ 19-6 Region, the pore space type was not separated.

4.3. Calibration of the Involved Model Constants Based on Reservoir Classification

To observe the correlation among φ_gm, φ_gsr and permeability, we classified the data sets that exhibited in Figure 7 into three clusters based on the difference of pore space type. Figure 8 exhibited the cross-plot of φ_gsr and permeability for three different types of formations. Comparison of Figure 7 and Figure 8 indicated that the relationship φ_gsr and permeability was absolutely different for different types of formations, formations should be first classified to calibrate reasonable input constants of C, p and q. Meanwhile, in reservoirs with fracture-pore type and pore-fracture type, the variation tendency between φ_gsr and permeability was similar, there is no need to further divide these two types of formations to calibrate the involved model constants.

Based on above analysis, we classified the exhibited data sets in Figure 7 into two categories, which were single fracture type and double-pore type (contained fracture-pore type and pore-fracture type), and they were used to determine the involved input constants in eq. 7. The calibrated models were expressed as follow:

$$\mathrm{Formation} \mathrm{with} \sin \mathrm{gle} \mathrm{fracture} \mathrm{type}\text{:} K=0.096\times {\phi }_{gm}^{1.69}\times {\phi }_{gsr}^{0.22}$$

(8)

$$\mathrm{Formation} \mathrm{with} \mathrm{double}-\mathrm{pore} \mathrm{type}\text{:} K=0.104\times {\phi }_{gm}^{0.73}\times {\phi }_{gsr}^{2.03}$$

(9)

5. Case Study

To use eq.s 8 and 9 to consecutively predict formation permeability, the formation pore type identification was a key. As was well known to all that the most direct data for indicating reservoir type was the cast thin section of core samples [36]. However, the limited number of core samples had restricted the use of them to continuously indicate reservoir types. In this study, we employed the method of core calibration logging, and analyzed the response characteristics of porosity frequency spectrum for different types of formations that indicated by the cast thin section. Finally, we found that the cross-plot of φ_gm and φ_gsr was the optimal indicator for formation type in the Archaeozoic Formation in BZ 19-6 Region (Figure 9). By using Figure 9, our target formations can be precisely classified into three types. For single fracture-type formations, eq. 6 can be used for permeability prediction, and for formations with double-pore type (fracture-pore type and pore-fracture type), unified model that expressed in eq. 7 should be used.

Combining with Figure 9 and eq.s 8 and 9, we processed the electrical imaging logging in several wells in BZ19-6 Region, consecutive porosity frequency spectra were acquired, and the permeability curves were predicted. Figure 10 exhibited the comparison of predicted permeabilities with core-derived results in a well. In this figure, electrical imaging logging acquiring from the Schlumberger’s FMI tool was displayed in the sixth track, and the extracted porosity frequency spectra were displayed in the seventh track. Meanwhile, we also displayed the extracted φ_gm and φ_gsr from porosity frequency spectra in the eighth track. Observation of the shapes of electrical imaging logging and corresponding porosity frequency spectra indicated that the interval of 4594.0-4600.0m belonged to pore-fracture type, the intervals of 4620.0-4623.0m and 4649.0-4653.0m were fracture-pore types, and the other intervals could be categorized as fractured reservoirs. The corresponding permeability prediction models were respectively used to consecutively calculate a permeability curve, as was displayed in the ninth track. Comparison of predicted permeability (PERM) with core-derived results (CPERM) exhibited a high consistency. The reliability of our raised model was verified. In addition, the extreme volatility of the predicted permeability curve indicated its high resolution, and this was derived from the electrical imaging logging.

6. Discussion

The shape of the porosity frequency spectrum was heavily determined by the pore space type. In formations with secondary pores (fracture and/or cave) developed, the porosity frequency spectrum exhibited as a relatively broad double-peak or multi-peak distribution, the corresponding φ_gm and φ_gsr would vary in a large range. However, in formations where only intergranular pores were present, the porosity frequency spectra exhibited as single unimodality, and their shapes were approximately that of standard normal distributions. In this case, the value of φ_gsr would be very close to each other in the whole intervals, even if they contained various porosities. The relationship between φ_gsr and permeability would be chaotic and disorderly. Hence, our raised model in this study was only valuable in formations where secondary pores had developed. In formations where the matrix intergranular pores play a dominant role, it would lose its stage.

7. Conclusion

Our target deep metamorphic rock reservoirs in the BZ19-6 Region contained three different types of pore spaces, which were single fracture type, fracture-pore type and pore-fracture type. This leaded to complicated pore structure and extremely strong heterogeneity. The relationship between core-derived porosities and permeabilities lost its role in permeability prediction. The electrical imaging logging could play an important role in these types of reservoirs. The extracted porosity frequency spectrum could well reflect the pore structure in fractured formations.

Based on the analysis of the morphological characteristics, two parameters, which were named as the logarithmic mean and the standard deviation between two golden section points of porosity frequency spectrum, were extracted. The relationships among these two parameters and permeability were analyzed, and a novel model to predict permeability in deep metamorphic rock reservoirs was established. The involved input coefficients were calibrated after formations were classified into two clusters.

Field applications in the Archaeozoic Formation in BZ 19-6 Region clearly indicated that our raised model was valuable in predicting permeability in fractured reservoirs. The predicted permeabilities matched well with core-derived results. Besides the metamorphic rock reservoirs, this model would also play a good role in any types of fractured reservoirs, e.g., the volcanic rock reservoirs, the carbonate rock reservoirs and fractured tight sandstone reservoirs.