Modified Low-Field NMR Method for Improved Pore Space Analysis in Tight Fe-Bearing Siliciclastic and Extrusive Rocks Open Access

The constant demand for energy resources requires detailed knowledge of the reservoir rocks. Understanding the architecture and distribution of pores within a reservoir is key to new hydrocarbon and geothermal water discoveries [1, 2]. Recently, numerous geothermal resources have been identified globally in tight rocks with a considerable share of microporosity. For instance, the tight lower Rotliegend formations of the Gorzów Block area located in West Poland were recently recognized as an especially attractive, unconventional source of geothermal heat [3].

To locate and manage such resources, the detailed fluid storage capacity and filtration properties of the reservoir rock need to be evaluated, and the precise determination of the pore size distribution (PSD) is required. Since tight rock formations are usually characterized by low porosity and permeability, a considerable share of microporosity, and the co-existence of various pore types, the quantitative characterization of the related pore networks is challenging and requires an appropriate choice of analytical methods. Such a characterization is of considerable importance [4], especially if it involves assessing a reservoir’s utility for geothermal energy extraction and Carbon Capture and Storage, as demonstrated in Enhanced Geothermal Systems [3, 5] or optimizing oil recovery in unconventional reservoirs through hydraulic fracturing [6, 7].

The most popular techniques used in quantitative pore space assessment include mercury intrusion capillary pressure (MICP) and low-field nuclear magnetic resonance (LF-NMR) [3, 8-11]. Both methods commonly complement each other but show different aspects of the pore space [12]. MICP employs nonwettable mercury as an injection fluid, reflecting only the connected pores and their pore throat sizes down to 3 nm, meaning that micropores (<2 nm) and capillary porosity are undetectable by means of this method [10]. In contrast, the LF-NMR data reflect the total pore space saturated with a wettable fluid, usually water, and are capable of penetrating the microporous space to supply data about the capillary porosity. The LF-NMR method has found broader applications regarding its nondestructive and noninvasive character and the ability to provide a wealth of precise information about tight rock samples such as, among others, the total and effective porosity, PSD, free and bound fluid content, and wettability [13]. However, NMR is not a direct method of measurement, and thus, the key to quantitative pore space characterization is performing a conversion between NMR-derived T2 relaxation times and pore sizes [14].

Despite the availability of numerous PSD conversion procedures [11, 15-17], with a particular emphasis on tight and shale rocks [18, 19], none of them offers a simultaneous solution to the following issues during conversion: (a) neglecting physical aspects controlling the relations between pore throats and pore bodies; (b) the use of simple linear models for complicated pore systems; (c) ignoring the clay and capillary-bound water contributions of a T2 distribution used as an input for porosity estimation and further PSD conversion; and (d) assuming that bound water always resides in small pores, while it could also be wetting the edges of larger voids according to the membrane-bound water model [15].

Another important issue in LF-NMR rock studies is the mineral composition. The transverse surface relaxivity (ρ₂) is considered a key factor for transforming the T2 relaxation time distribution into a PSD as it defines the surface relaxation phenomenon [20]. Since surface relaxivity is a lithology-dependent parameter [21], the content of ferro- and para-magnetic elements that are typically found in clay and iron minerals is a major factor controlling an increase in the surface relaxation rate in porous rock media [22-24]. While the surface relaxivities of major rock-forming minerals such as quartz and calcite are well documented [22, 25-28], the variability of this parameter for clay minerals such as illite, chlorite, smectite, and mixed-layers illite-smectite is still not well understood [21, 29-33]. For this reason, the complex mineralogy of extrusive and siliciclastic rocks and the varied structure of clay minerals present in their matrix can potentially cause drastic changes to the LF-NMR PSD estimation process and their presence cannot be neglected.

Considering the above-listed issues, this article aims to improve the accuracy of PSD estimation in tight extrusive and siliciclastic rocks by setting up a new MICP-NMR conversion procedure based on the percolation theory combined with nonlinear regression analysis and a modified differential LF-NMR protocol. In this protocol, the time domain signal recorded for the dry sample is subtracted from that of a 100% water-saturated one before applying the Inverse Laplace Transform (ILT). As a result, the differential T2 relaxation distribution obtained corresponds directly to the effective porosity since the information from protons bound in a crystal lattice of matrix particles, associated with organic matter, clay matter, and closed porosity are uniformly removed throughout a 100% water-saturated spectrum. The protocol was first proposed by Habina et al. [34] but has never been applied for PSD conversion in tight rocks before and is thus tested for this purpose. It has a major advantage over the standard NMR T2 cut-off approach because the exact quantification of the clay-bound water (Clay_bw) in the latter is not possible in tight rocks. The differential distributions should correlate much better with the MICP data while including information about micropores and capillary pores beyond the reach of the porosimetry methods. To underline the strong points of using the differential NMR T2 distributions, we also compare the corresponding results obtained using both NMR protocols. The two scenarios are verified on selected tight siliciclastic and extrusive rock samples.

The tested materials comprise siliciclastic rocks (thirteen tight sandstones and three conglomerates) and extrusive rocks (seventeen samples) of the Lower Permian and/or Carboniferous age. Samples were taken from the Rotliegend formation of the Ośno IG-2 well located in West Poland, targeting the geological strata of the so-called Gorzów Block unit. The corresponding depth range of their extraction is 3212–3658 m b.s.l [3].

For the needs of the LF-NMR and MICP porosity measurements, thirty-three cylindrical cores (sixteen siliciclastics and seventeen extrusives) with a diameter of 2.2–2.5 cm and a length of 2–5 cm were cut from the drill core. Twenty-eight samples (eleven siliciclastics and seventeen extrusives) were selected for X-ray diffraction (XRD) analysis. To illustrate the T2-PSD conversion process, a subset of three samples from both the siliciclastic and extrusive groups were chosen. The selected samples represented different lithologies and featured a diversified character of the corresponding cumulative intrusion (MICP) and saturation (LF-NMR) curves.

The XRD measurements were carried out using the Panalytical X’Pert Pro apparatus equipped with an advanced ultra-fast detector (real-time ultra strip X’Celerator). The operating conditions were set at a voltage of 40 kV, a current of 40 mA, and a step width of 0.02° 2Θ. The scanning range was from 5° to 65° 2Θ.

Before the measurements, all samples were ground using a McCrone micronizing mill to achieve a grain size below 20 µm. Following [35], ZnO was incorporated as an internal standard of the analysis before grinding to ensure the thorough homogenization and quality of the analysis. Measurements were taken on randomly prepared samples by side-loading.

A quantitative mineral composition analysis was performed using the Rietveld method [36], facilitated by the SIROQUANT program, which is particularly effective for samples containing clay minerals.

MICP measurements were performed using an Auto Pore IV 9520 (Micromeritics, USA) mercury porosimeter. Before the analysis, all samples were crushed and dried at 105 °C for 24 hours to remove moisture from the pore space. The pressure was measured at 120 points in the 0.003–414 MPa range, which corresponds to pore sizes from about 500 µm to 3–4 nm. Using the MICP method, the effective porosity (φ_MICP) and PSD of the studied samples were determined. Equivalent pore throat diameters were computed based on the corresponding capillary pressure using equation (1) [37]:

where d_throat is the pore throat diameter (µm); γ_Hg is air/mercury surface tension (γ_Hg = 485 mN/m); θ_Hg is the contact angle between mercury and a porous material (θ_Hg = 130°); and P_c is capillary pressure (MPa).

To conduct the LF-NMR analysis, core samples were prepared in both a dried and saturated state. Dry samples were obtained by vacuum oven drying at 110 °C and under the pressure of 0.007 bar for 12 hours to remove free and capillary-bound water [38]. Samples in a saturated state were obtained by saturating them in demineralized water in analogous vacuum conditions for 12 hours.

The LF-NMR measurements were performed using a 2 MHz Magritek Rock Core Analyzer (Magritek, Germany) with a 29 mm probe. Experiments of the T2CPMG (The Carr–Purcell–Meiboom–Gill) sequence were carried out using the following acquisition parameters: repetition time, RT = 2000 milliseconds; echo time, TE = 60 μs; number of scans, NoS = 512; and number of echoes, NoE = 20,000. The LF-NMR output data were analyzed using the Prospa software package. T2 relaxation times for all samples in both dry and 100% water-saturated states were measured. Differential T2 distributions were obtained using a modified differential T2 protocol.

In the standard T2 differential protocol, the “background” signal registered for the probe is subtracted from the signals registered for samples in dry and saturated states. Then, the T2 distributions are calculated using the ILT. The differential T2 distributions used in this work are obtained by applying a modified differential protocol. In this method, the T2 spin-echo signal originating from the rock matrix and closed pores (dried sample signal) is subtracted from the signal of the 100% water-saturated sample and then converted using the ILT. Assuming the same background noise is maintained for samples in different saturation states, this approach preserves higher signal-to-noise ratio and thus provides more accurate ILT results [34]. In addition to this study, the employed modified differential protocol has been widely adapted in pore space characterization of various lithology types, including natural silica [26], shales [10], cherts [27], and carbonates [39, 40].

The differences in the distributions obtained using the standard T2 cut-off approach and modified differential methodologies are illustrated in Figure 1. The conventional technique entails the utilization of T2 distributions and cumulative saturation curves derived using dried and 100% water-saturated samples following the standard T2 protocol. Using this technique, the assumption is made that chemically bound water resides exclusively in small pores (according to the small pore-bound water model [15]). Consequently, to obtain the effective porosity (φ_std), the T2_clay cut-off value determination was performed to quantify the contribution of Clay_{bw_std} to the overall 100% water-saturated state distribution (Figure 1(a)). Accordingly, employing the standard cut-off approach and the oven-drying method, the irreducible water saturation index (SWI_std) was calculated using the following equation (2):

In the modified differential approach, we utilized a new T2 distribution obtained by subtracting the signal of the dry sample from the signal of the 100% water-saturated sample. By the nature of this protocol, we obtain the parameters by the assumption of the membrane-bound water model [15]. Importantly, the resulting differential T2 distribution is devoid of clay-bound water signals and only contains information about effective porosity. Hence, for the cumulative T2 curve of the differential distribution, effective porosity (φ_diff) could be obtained as its maximum value. To determine the exact volume of Clay_{bw_diff}, the corrected T2′_clay cut-off times were derived using the methodology described by Habina et al. [34] (Figure 1(b)). Therefore, to utilize the data acquired through the modified differential approach in the appropriate relation to φ_total, the irreducible SWI_std for the T2 differential distribution (SWI_diff) was calculated using equation (3):

The application of NMR in pore system analysis is based on the assumption that, in a uniform magnetic field, the diffusion relaxation time T2_D of fluids saturating geological media in a fast diffusion limit can be neglected, and the bulk relaxation time T2_B of bulk water is much larger than T2_S. Therefore, the relaxation time is mainly influenced by surface relaxation T2_S [41] as shown in equation (4):

where T2_B is the bulk relaxation of saturating fluids (ms); T2_S is the surface relaxation caused by the interactions between saturating fluids and pore surfaces (milliseconds); T2_D is the diffusion relaxation (milliseconds). In our case, we applied a very short echo time (see Section 2.3) and induction of the magnetic field B₀ was equal to just 0.05 T. Therefore, the T2 relaxation time could be related to pore geometry through the surface relaxivity ρ₂ and pore surface to pore volume ratio, in a relation commonly expressed by equation (5) equation (5)[20]:

where ρ₂ is the transverse surface relaxivity (μm/s); and S/V is the pore surface to pore volume ratio (1/µm), which, assuming that the pore shape is regular, can be related to shape factor Fs and a pore body size d_por as shown in equation (6):

where d_por is the pore diameter (µm); and Fs is the geometric factor of the pore shape; here, we assume Fs for spherical geometry such that Fs = 6.

The proposed methodology for the PSD conversion is outlined in Figure 2. Once the input data from the LF-NMR and MICP are prepared for PSD estimation, a framework for T2 and pore size correlation needs to be established. This framework aims to determine the conversion points based on the differences between the two experimental methods and the specific characteristics of the investigated tight pore space. In this study, we adopted the percolation theory, as suggested by Daigle and Johnson [16] and Zhao et al. [19], to determine the conversion points.

According to the percolation theory, we assume that for pores of diameter d, the mercury intrusion cumulative curve represents the occupation probability of bonds, denoted as X(d), and the NMR T2 inverse cumulative porosity curve represents the accessible fraction, denoted as X^A(d). Then, based on these assumptions and the simulation results obtained by Liu et al. [42], we use the following relation equation (7):

where $Z$ is the coordination number, computed from the relation to the percolation threshold $fc$ [43] shown in equation (8):

The percolation threshold f_c was determined at MICP curves at a pressure for which mercury managed to penetrate the connected pore space of the sample for the first time, hence forming the first sample-spanning cluster. After this step, the MICP occupation probability of bonds and LF-NMR accessible fraction could be calculated using equations (7) and (8). Thus, the correlation framework was determined with the consideration of SWI_std and SWI_diff in specifying the proportion of the MICP cumulative curve to LF-NMR 100% water-saturated and differential cumulative curves, respectively.

Once the correlation framework is established, the selected conversion points were correlated with consideration of the fact that the NMR T2 time distribution corresponds to the pore size (d_por), whereas MICP is primarily influenced by pore throat size (d_throat). Addressing this difference, we adopted the approach introduced by Marschall et al. [44], utilizing the effective surface relaxivity coefficient ρ_e to replace the transverse surface relaxivity ρ₂ in equation (5) based on the relation presented in equation (9):

Therefore, transverse surface relaxivity ρ₂ is proportional to the product of the effective surface relaxivity ρ_e and the ratio of pore body size to pore throat size.

Then, based on equations (6) and (9), equation (5) could be rewritten to the form shown in equation (10):

To establish the relationship between LF-NMR pore size (d_por) and MICP pore throat size (d_throat), following the works of [11, 45, 46], we employed the power function shown in equation (11):

where m and n are the proportionality coefficients. Then, based on equations (10) and (11), assuming C = m(ρ_eFs)ⁿ, the pore size d_por could be calculated using equation (12):

where C is the conversion coefficient for the nonlinear T₂-PSD conversion. Then the unknown coefficients C and n were determined through nonlinear regression analysis between the LF-NMR pore size (d_por) and MICP pore throat size (d_throat). To achieve this, we calibrated the model function by taking logarithms on both sides of equation (12) and minimizing the sum of squared errors between the converted pore sizes and pore throat sizes for m conversion points using equation (13):

After calibration, with the obtained conversion coefficients C and n, the NMR T2 distribution could be converted to absolute PSD using equation (12) with corresponding pseudo capillary pressures calculated using transformed equation (1). The obtained PSDs were divided into microporosity, mesoporosity, and macroporosity according to the IUPAC classification [47].

As the reference for LF-NMR porosity determination methods, mass–volume measurements obtained during the sample preparation process were used. Mass–volume porosity (φ_mass-vol) was calculated by comparing the bulk density of a sample (mass per unit volume of 100% water-saturated sample) to its true density (the mass per unit volume of a dry sample), according to equation (14):

where m_sat and m_dried are the mass of saturated and dried samples, respectively (g); p_w is the density of water (g/cm³); and V_sample is the volume of the sample (cm³).

The relative porosity error (φ_RE) was calculated between the effective porosity from the modified differential method and the mass–volume measurements:

A detailed geological analysis of the study area has already been presented by Sowiżdżał et al. [3], and therefore, we only briefly discuss the XRD characteristics of the samples here. In the description below, labels P1–P14 refer to tight sandstones, while P16–P17 signify conglomerates and numbers above 23 correspond to extrusive rocks.

The lithological analysis performed by Gajewska [48] showed that the tight, Rotliegend sandstones from the Ośno IG-2 well are well cemented, exhibit a fine-grained texture and contain numerous siltstone intercalations often observable as laminations, as well as lenses of coarse-grained, heterogeneous sandstones. Two of the samples chosen for PSD analysis, P02 and P06, were classified as such. They could be considered subarcosic, tight sandstones, as showing less than 15% of feldspars (see Figure 3). They are cemented by calcite, the content of which is between 15% and 20% and contain negligible admixtures of clay minerals, mainly micas and chlorite. Although hematite was detected in those samples, its share does not exceed 6%, while micas are quite an important constituent of the tight sandstones group, with share up to 10%. As shown in Figure 3, however, the admixture of hematite can be significant in the sandstone rocks considered and sometimes even exceeds 15%. The rocks of this group also have a very variable content of calcite. Among the clay minerals, chlorite occurs in all the samples, but in rather low amounts compared to illite, another common mineral in the whole siliciclastic group. Pure smectite and illite–smectite occur occasionally in exceptionally low amounts.

The conglomerates mainly consist of angular to subangular clasts of volcanic rocks, quartz, and feldspars and are characterized by a poorly cemented matrix. This lithology has a chaotic character and displays poor sorting. This rock group is characterized by a high illite share (often above 10%) and more exposed admixture of smectite while maintaining similar amounts of chlorite in comparison with the tight sandstones under study. Another characteristic feature of this group is a growing share of hematite. The P17 sample chosen for PSD analysis represents this lithological type (Figure 3).

The extrusive rocks encountered in the well mainly comprise volcano-clastic breccias and dacites with a vesicular texture and varying degrees of hydrothermal alteration and iron oxide impregnations. The intensity of fracturing is variable within the dacites, and some areas display cataclastic textures. Moreover, carbonate veins are present, filling fractures up to 2 mm in thickness. The share of quartz in some cases is much lower than for siliciclastic rocks, and the dominant mineral becomes plagioclase. Among the clay minerals, illite is the most common for the extrusive rocks under study, with a maximum content even approaching 23%. In this group, the content of smectite minerals is also significant relative to the tight sandstones, occurring in most samples and reaching even up to 19%.

We chose P29, P35, and P40 samples to represent this group in the PSD analysis. They are composed of various mineral phases ranging from quartz, trough feldspars, hematite, micas, and sulfates to clay minerals. Overall, the share of hematite seems higher in comparison with siliciclastic rocks, and it often comes close to or exceeds 10%. In the P29 and P35 samples, quartz is the dominant phase which is not the case for the P40 sample of the most complex mineralogical composition. In addition to the already mentioned minerals, around 2% of dolomite occurs here as well, accompanied by subordinate amounts of ankerite.

In summary, the amount of quartz is subject to dynamic variations in samples from all lithological groups and most in extrusive rocks. All samples contain some admixture of clay minerals and hematite. Those phases are the most common for extrusive rocks and conglomerates.

The T2 relaxation time distributions of the 100% water-saturated and dry samples (Figures 4(a) and 4(b)) as well as corresponding differential distributions (Figure 4(c)) and the MICP-derived PSD (Figure 4(d)) are presented for the chosen representative sample set of tight sandstones (P02 and P06), conglomerate (P17), and extrusive rocks (P29, P35, and P40). The quantitative details of the T2 NMR distributions and MICP PSD are included in Table 1. Complete NMR porosity estimation results and their statistical analysis for all samples are shown in Table 2, while Table 3 shows the statistical analysis for the experimental results obtained using MICP.

The LF-NMR measurements provided valuable data used for the quantitative analysis of the pore space in the selected tight rocks, covering a wide range of T2 relaxation times. The T2 relaxation time distributions of a 100% water-saturated state exhibited multimodal characteristics in the case of siliciclastics and bi- to uni-modal characteristics in the case of extrusives (Figure 4(a)). The corresponding logarithmic mean T2 (T2lm_std) values ranged from 0.143 to 1.607 milliseconds, with extrusive rocks showing shorter T2lm values compared to siliciclastics. This trend is also observable in differential distributions, where the T2lm_diff shifts toward longer times ranging from 0.287 milliseconds for extrusive rocks to 2.299 milliseconds for siliciclastic rocks (Table 1). Differential distributions also further highlighted the uni- to bi-modal characteristics of extrusive rocks and multimodal characteristics of siliciclastics (Figure 4(c)). The dry-state samples displayed dominant peaks shifting toward shorter T2 times (<0.15 milliseconds), and in all cases exhibited unimodal distributions with many small (<0.05% incremental porosity) and separated peaks within longer T2 time ranges (Figure 4(b)). The corresponding T2lm_dry values ranged from 0.057 to 0.123 milliseconds, with extrusive rocks characterized by both the shortest and the longest logarithmic mean times (Table 1).

The MICP data analysis also revealed the distinct characteristics of different rock lithotypes. In tight sandstones, the pore space exhibited uni- to bi-modal characteristics, with pore systems observed across a wide range of pore sizes. The pore space in conglomerate sample P17 and extrusive samples displayed a unimodal character, dominated by mesopores (<0.05 μm; Figure 4(d)). The presented results were consistent with the median pore diameter values ranging from 0.008 to 0.026 µm and 0.026 to 0.164 µm for extrusive and siliciclastic rocks, respectively (Table 1).

In the case of the standard LF-NMR approach, the pore space components were derived using the T2_clay cut-off values with a mean of 0.27 milliseconds for all the investigated lithological groups. Among the studied samples, siliciclastics exhibited the lowest mean total porosity (φ_total) value of 8.9%. In contrast, extrusives displayed much higher values of φ_total equal to 12.2%. When considering the Clay_{bw_std} parameter, the extrusives had the highest mean value of 6.1%, followed by siliciclastic rocks (3.0%). The effective porosity (φ_std) demonstrated the highest mean value for extrusive rocks (6.1%), while tight sandstones and conglomerates of the siliciclastic group showed similar, significantly lower mean values of 5.9% (Table 2).

For the modified differential approach, T2′_clay was incorporated as the cut-off value for deriving porosity parameters. In this case, the average values reached significantly shorter times, equal to 0.17 and 0.19 milliseconds for siliciclastic and extrusive rocks, respectively, while having a significantly lower standard deviation reaching a maximum of only 0.08 milliseconds for siliciclastics. The Clay_{bw_diff} parameter derived using the differential method also exhibited the highest mean value for extrusive rocks, despite its relatively low value of only 4.6%, followed by siliciclastics (2.3%). The differential effective porosity (φ_diff) parameter, analogous to φ_std obtained from the standard approach, also demonstrated the highest mean value for extrusives (7.6%) and the lowest mean value for siliciclastics (6.6%). However, both values are significantly higher than those that would be obtained using the standard method (Table 2).

In terms of the MICP experimental results, siliciclastics exhibited the lowest mean total pore area and effective porosity. The average total pore area for sandstones and the conglomerates group was 6.619 m²/g, while the mean φ_MICP was 5.8%. In contrast, extrusive rocks had a higher mean total pore area of 10.417 m²/g, and similarly to the NMR results, they displayed the highest mean φ_MICP of 7.7% (Table 3).

The mean threshold pressure measured during the experiments was relatively low for all lithologies and ranged from 0.35 MPa for extrusives to 0.371 MPa for siliciclastics. The corresponding mean permeabilities were found to be the highest for siliciclastic rocks (0.114 mD) and the lowest for extrusive rocks (0.061 mD). However, it should be noted that the permeability values for siliciclastics had the highest standard deviation of 0.067 mD. Conversely, extrusives displayed the lowest SD of only 0.020 mD (Table 3).

A comparison of the obtained porosity values with the reference mass–volume measurements showed that φ_std achieves the highest coefficient of determination, as well as the lowest RMSE for the x = y curve (Figure 5(b), Table 4). The situation is different for the differential method, where the results of the linear regression analysis indicated that the porosity data become less scattered and more consistent around the fitting line. Another effect associated with applying the differential method is related to narrower prediction and confidence bands. Consequently, the corresponding regression allowed for increasing the coefficient of determination by as much as 0.31 reaching an R² of 0.88 compared to φ_std with a change of only 0.02 and an R² of 0.72. Simultaneously, φ_diff achieves the lowest RMSE during the regression, reaching only a porosity prediction error of just 0.8%. Nevertheless, one can still see a clear, almost parallel shift (slope equal to 1.02) of the fitted function with respect to the x = y curve (Figure 5(c); Table 4). In such a case, the slope coefficient of the linear function can be ignored, which means that the intercept determines an almost systematic porosity error of about 0.93%. At the same time, the comparison of MICP and mass–volume data showed the weakest correlation, confirmed by the lowest R² and the highest RMSE for both the y = x curve and the linear fit (Figure 5(a); Table 4).

Since we noticed exceptionally low T2 relaxation times in our data (Table 1) and almost systematic errors in the differential effective porosities obtained (Table 4), we decided to analyze the influence of mineral composition on the results of LF-NMR measurements. Of all the minerals analyzed, only hematite and clay minerals: illite and smectite significantly affected the variation in parameters obtained by LF-NMR (Figures 6(a)–6(f)). No correlation was observed between chlorite content and NMR parameters. In the case of micas, only a weak, negative correlation with relative porosity error was identified (online Supplementary Material Figure S1).

It was found that the T2 relaxation time of the analyzed samples was influenced the most by the content of hematite—a ferromagnetic mineral that contributes to the disappearance of total magnetization (Figure 6(a)) and clay minerals (illite and smectite)—paramagnetic phyllosilicates known for their ability to swell when exposed to water (Figure 6(b)). In the case of siliciclastic rocks, a power relationship between the sum content of these minerals and the results of T2 measurements was found, explaining as much as 85% of the variance in T2 log-mean relaxation time. In the case of extrusive rocks, there was only a linear influence with a two times lower coefficient of determination equal to 0.44 (Figure 6(c)).

A different situation was observed in the case of relative porosity error. We observed a linear increase in porosity errors which was dependent on hematite content for both extrusive and siliciclastic rocks (Figure 6(d)). However, for clay minerals, only smectite affected this parameter. We noted a power-law increase in the relative effective porosity deviation for extrusive rocks. Siliciclastic rocks, on the other hand, are characterized by a low content of these minerals (see Figure 3), so in their case, the effect is not apparent (Figure 6(e)). Nevertheless, the aggregated contents of hematite and smectite explain as much as 59% and 54% of the variance in relative effective porosity deviation for siliciclastic and extrusive rocks, respectively (Figure 6(f)).

The analysis of the LF-NMR data involved determining the optimal conversion coefficients C and n for both the differential and 100% water-saturated T2 relaxation time distributions (Table 5) based on conversion points determined according to the framework established using percolation theory and LF-NMR SWI data. These conversion coefficients were crucial for further conversion of the T2 relaxation time into the PSD comparable with the results obtained from MICP measurements.

The results of the regression analysis yielded the values of conversion coefficients (C) ranging from 0.007 to 0.1 and from 0.023 to 0.104 for the 100% water-saturated and differential distribution, respectively. The values of the conversion exponents (n) varied between 0.96 and 1.83 for the 100% water-saturated distribution and between 1.02 and 2.70 for the differential distribution. Overall, the obtained coefficients of determination between estimated and measured pore sizes are high, ranging from 0.90 to 0.98 in all cases (Table 5). This proves that the estimation of the parameters of the nonlinear model was successful, and the created models in every case explain more than 90% of the variance in pore throat diameters measured by the MICP in the regime of the designated for each sample percolation theory and SWI data frameworks.

By creating LF-NMR PSDs based on the obtained conversion models, a direct comparison with the distributions obtained by the MICP method was possible. In siliciclastic rocks, the estimated PSD showed a bi- or multi-modal pore distribution, with the most numerous populations coinciding with the mesopore regime and often multiple populations falling into the category of macropores. The pore distribution for clastic rocks ranged from 0.003 µm (2 nm) to 300 µm for the differential curves and from <0.0001 μm (<0.1 nm) to 600 µm for the standard curves (Figures 7(a)–7(c)).

Nevertheless, the differential curves showed a good agreement with the MICP measurements. However, a weak agreement was observed for sample P06, where the PSD significantly exceeded the distribution registered by MICP in the regime of large pores (Figure 7(b)). Still, a very good fit was observed in the regime of small pores. On the other hand, the 100% water-saturated curves did not provide a reliable PSD estimation for most siliciclastic rocks, especially in the case of the abovementioned sample, where the estimated PSD covered unrealistic ranges for both large and small pore spaces (Figure 7(b)).

The extrusive rocks displayed distinct pore size characteristics, distinguished by unimodal distributions in the mesopore regime. The pore distributions for samples in this lithology group ranged from 0.0009 µm (0.9 nm) to 70 µm for the differential curves and from <0.0001 μm (<0.1 nm) to 70 µm for the 100% water-saturated curves (Figures 8(a)–8(c)). In such cases, the differential curves also yielded PSD results that aligned better with the MICP measurements, except for sample P40 where the MICP results revealed residual amounts of macropore spaces not observed in the LF-NMR data (Figure 8(c)). In contrast, the PSD obtained from the 100% water-saturated curves for extrusive rocks exhibited unrealistic diameter ranges, particularly close to and below 0.0001 µm (0.1 nm) for all samples (Figures 8(a)–8(c)).

The pseudo capillary pressure curve obtained from the 100% water-saturated LF-NMR PSD always corresponded with a maximum saturation level of 100%, reaching pressure values >1000 MPa. However, the pseudo-Pc curve obtained from the differential approach showed much better agreement with the actual MICP results, reaching a maximum saturation level of 83% for sample P02 (in the case of siliciclastic rocks) and 70% for sample P35 (in the case of extrusive rocks) at the maximum pressures of 181 and 456 MPa, respectively. Similarly to PSD observations, the comparison of capillary pressure curves also revealed distinct misfitted spots, particularly at low-pressure ranges corresponding to larger pore throats and long T2 times, especially for samples P06 (Figure 7(b)) and P40 (Figure 8(c)). Moreover, one can see large differences in the maximum saturation level of capillary pressure curves measured by MICP, depending on the SWI_std or SWI_diff value used in the calibration process, reaching up to about 15% in the case of samples P02 (Figure 7(a)) and P35 (Figure 8(b)).

As highlighted above, NMR relaxometry yields relative PSDs expressed in time units by default [20]. Such a representation is useful for a tentative evaluation and clustering of pore groups. However, knowledge of the absolute dimensions of voids denoted in specific units of length is often required for petrophysical reservoir characterizations. This requirement stems from the fact that the lower limit of pore sizes is crucial in stating whether a given fluid such as crude oil, natural gas, or formation water can still pass through a particular pore network [49]. Another aspect is that the larger the pores, the easier it is to assess the porous media they host without the need for excessive artificial reservoir stimulation, such as hydraulic fracturing or acidization.

The topic of pore size conversion from time scale to metric units of length has long been critical in the scientific literature related to petrophysics. Some researchers have proposed employing diffusion NMR measurements for that purpose and extracting S/V ratios by observing diffusion constant changes under different echo times [50, 51]. However, the accuracy of diffusion-based conversion is debatable since it is often assumed that surface relaxivity remains the same for all grains and that the influence of internal magnetic field gradients is negligible [24]. Because of this fact, many researchers selected another way to convert the time-domain PSD into absolute distribution by combining NMR and MICP data, for example [52]. A potential issue in the MICP-based approach is that both methods detect different structures. As mentioned earlier, the MICP method is sensitized to pore throats, whereas NMR mostly shows the distribution of pore bodies [53]. However, Volokitin et al. [14] observed that the shape of the NMR T2 distribution curve for 100% water-saturated, water-wetted rock is similar to the pore throat size distribution converted from capillary pressure curves. Without additional mathematical data manipulation, such a conversion would only be reliable provided that the size differences between the apparent pore bodies and throats are relatively small and proportional [54]. To partially avoid this restriction, one solution is to establish a specific nonlinear conversion function between pore throat and pore body sizes before transforming the NMR T2 distributions into a PSD or pseudo-Pc curve using the Washburn equation [11, 16, 18, 19, 44, 55, 56]. During this step, it is important to identify which points on the MICP and NMR intrusion and saturation curves should be used for the conversion process. For instance, Xiao et al. [11] used piecewise conversion functions, differentiating between small and large pore sections on saturation curves by means of so-called “turning points.” Zhou et al. [18] developed this method further by employing a Bayesian Regularization Neural Network. However, as finding precise locations of the turning points for heterogeneous samples is a difficult task, other authors have presented a completely different approach. Daigle and Johnson [16], followed by Zhao et al. [19], established the correlation framework between the T2 and d_throat in the high-pressure partition of the MICP curve and short NMR T2 times using the percolation theory, the use of which minimizes the potential conversion errors and improves the accuracy of the whole process, especially in small pore sections of PSD.

Although the fundamentals of both percolation theory and MICP-based NMR PSD conversion are already well established in the literature, to the best of our knowledge, this concept has not yet been commonly applied directly to derive the absolute PSDs from NMR in tight rocks. This was only realized by Zhao et al. [19], but the cited authors decided to use a linear regression approach for small pores and surface relaxivity computed from Kozeny’s equation for large pores. On the other hand, some authors have used a more precise, nonlinear approach which was not supported by the percolation theory itself [11, 45, 46]. Therefore, in our opinion, filling this gap is not only crucial from a scientific point of view but also from the industrial perspective since spotting potential resources with wrongly computed PSD conversion models would generate overwhelming operational costs. At this point, we should recall that another enhancement proposed in this article is the application of the modified LF-NMR protocol to obtain the T2 differential curves used for the conversion itself as well as SWI parameters used for the MICP calibration process.

Conversion itself plays an unprecedented role in obtaining the correct absolute PSD but there is also another major problem which seems to be often ignored during this process, namely the choice and preparation of a suitable T2 distribution. In many cases, the distribution of the sample in a 100% water-saturated state is chosen, which is acceptable for conventional reservoir rocks generally devoid of organic matter and significant amounts of capillary and clay-bound water as well as para- and ferro-magnetic mineral content. However, since NMR is sensitive to ¹H hydrogen atoms, the T2 time distribution curves of a 100% water-saturated sample used in the standard methodology also contain information about the chemically bound water, organic matter, and closed pores and is affected by the influence of field gradients induced by the presence of minerals with magnetic properties [10, 21, 22, 24, 34]. Therefore, using these distributions and the related cumulative T2 curves obtained during the T2 experiments as the inputs for the PSD conversion can lead to large estimation errors. In other words, PSD obtained this way may be overprinted by the signal from hydrogen nuclei which are not necessarily associated with the effectively connected pores such as those depicted by MICP and thus significantly shifted from their true ranges.

In rocks with relatively simple anatomy, oven drying can remove nearly the entire volume of capillary-bound water [38]. Therefore, in this article, we decided to use differential T2 distribution obtained according to Habina et al. [34] by means of a modified differential protocol, which in subsequent years was successfully applied to a set of shale rock cores with ultra-small pores, with dominant sizes of the order of just 2 nm [10, 57].

The results of the pore characterization obtained for the studied tight sandstones, conglomerates, and effusive rocks revealed notable differences in the measurement and interpretation capabilities of the effective porosity obtained by the two LF-NMR data preparation protocols employed. While the standard protocol for porosity evaluation is commonly chosen, it faces obstacles in distinguishing signals from chemically bound water in rocks with tight pore spaces [10, 34].

In general, the standard protocol proves to be unsuitable for porosity measurements in tight rocks with high clay content due to the following reasons: (a) residual water that was not removed during the drying process (closed pores) remains visible on the T2 distribution graph, at T2 > T2_clay cut-off and (b) remains of water that was assumed to be removed during the drying process (free water and capillary-bound water) is still detectable on the T2 distribution graph, at T2 ≤ T2_clay cut-off (see Figure 1).

The modified differential protocol incorporated a more accurate clay-bound water porosity and provided additional insights into the presence of clay minerals in the pore spaces. For instance, remarkably high mean values for Clay_bw in extrusive rocks obtained by the standard protocol indicate a higher content of clay minerals in their pore spaces, while the differential method showed lower amounts of approximately 1.5% (Table 2). This affected the estimation of effective porosity and the quality of its correlation with the reference mass–volume porosity measurements. For the standard method, the obtained values showed a maximum coefficient of determination equal to 0.72, while this coefficient reached as much as 0.88 for the modified differential method (Figure 5).

At the same time, a near-parallel shift of the linear fit with respect to the x = y curve was observed for the modified differential method-derived effective porosities. This phenomenon was not observed for the values obtained by the standard method. As noted by Elsayed et al. [30, 31], such an effect may be due to the appearance of nonphysical peaks related to a change in the diffusive regime (due to strong internal gradients), especially in nano/micropores, with the presence of minerals containing para- and ferro-magnetic elements. For this reason, signals from pores in those size regions for dry and 100% water-saturated distributions no longer correspond to each other. This, in turn, can lead to the generation of artificial noise during signal subtraction in the modified differential protocol and thus signal amplification during the ILT transform.

In addition, we observed a relationship between the content of clay minerals and hematite, and the relative error in estimated differential effective porosity values (Figure 6). We also found that this error is almost systematic, which makes it straightforward to correct by using a linear fit function to mass–volume data (Figure 5; Table 4). Both corrections, based on mass–volume, as well as hematite and smectite content equations, could be executed, yielding comparable results (online Supplementary Material Figure S2).

The regression analysis showed that the conversion points from both differential and standard LF-NMR curves successfully matched the PSD obtained from the MICP measurements. However, it should be noted that the ultimate results of the PSD conversion generated using the data from both standard and differential modified methods did not yield comparable outcomes, despite the relatively high values of the determination coefficients for conversion functions. This is particularly evident in the case of samples P02 and P06 (tight sandstones) and P29 and P40 (extrusive rocks) (Figures 7(a) and 7(b); Figures 8(a) and 8(c)).

This discrepancy was also evident when comparing the capillary pressure curves transferred from the LF-NMR standard distribution with the MICP capillary pressure curves. The length and course of the curves showed clear variations, especially in the small pore section. On the other hand, the capillary pressure curves transferred from the differential distribution could be calibrated by SWI_diff data and therefore exhibited better fitting with the MICP curves which is particularly discernible for small pores.

Moreover, the estimation of PSD from standard curves failed in most cases, particularly in samples P06 (tight rock), P29, and P40 (extrusive rocks). The resulting pore diameters obtained from the standard curve often yielded unrealistic values and included unwanted contributions from chemically bound water. This resulted in shifting the distribution to the left relative to the PSD obtained from the MICP measurements.

It is also noteworthy that the LF-NMR distributions of the P06 sample are characterized by the presence of large pores invisible on the MICP curves (Figure 7(b)). Unfortunately, the presence of such phenomena seems unavoidable, as MICP and NMR differ in the context of sample preparation and the physics of their measurement. As a result, regardless of the PSD estimation method used, there may always be cases of nonmatching distributions that will negatively affect the data conversion process. To minimize this effect, it seems beneficial to employ PSD models using polynomial or combined functions containing more than one value of surface relaxivity [11, 19]. Nevertheless, this effect was only observed for one tight sandstone sample and is not expected to occur for extrusive rocks characterized by rather unimodal distributions of relaxation times and pore sizes.

In this study, a strong diversification of fluid storage and filtration properties was noted between the analyzed rock lithotypes. Different PSDs and pore distribution patterns emphasize the importance of considering the heterogeneity of the pore space when assessing the reservoir properties. Employing a modified differential approach should overcome most of the issues related to porosity and PSD determination using the LF-NMR technique, especially in tight rocks. Even though para- and ferro-magnetic minerals interfere with the execution of the signal subtraction procedure, the generated porosity error can be considered systematic and may be easily corrected. A much better correlation between the differential LF-NMR effective porosity and the mass–volume porosity compared to the standard NMR approach and MICP data indicates that the proposed approach is very efficient in delivering the absolute PSD in tight rocks. The regression analysis between NMR and MICP data suggests that the T2 times from both NMR protocols are well correlated with the pore throat distributions. Despite this, the PSD curve obtained based on the standard NMR data may be deceptive as it does not differentiate between the actual open pores and the chemically bound water. Consequently, the PSD conversion results derived by the LF-NMR differential protocol matched the MICP data much better in comparison with the standard approach, especially in the case of conglomerates and extrusive rocks characterized by unimodal PSDs. Consequently, the output of the LF-NMR-based PSD estimation is only compatible with MICP when differential distributions are used. They are devoid of unwanted hydrogen nuclei contributions and represent a full spectrum of pores accessible to wetting fluids such as hydrocarbons and formation waters, including capillary pores inaccessible to mercury.

Our research further indicates that the relation between the T2 NMR relaxation times and MICP pore sizes is not linear and difficult to describe with one universal equation, especially for tight rocks with high para- and ferro-magnetic mineral content. Unfortunately, as observed by other authors, the use of the percolation theory in PSD conversion was only successful for rocks with unimodal T2 and PSD distribution containing predominantly micropores and mesopores. For rocks with larger voids, the conversion using a single surface relaxivity approach usually failed, especially in the case of tight sandstones. Moreover, even the use of the LF-NMR method with very low echo times and the implementation of the modified differential procedure does not completely eliminate the influence of para- and ferro-magnetic minerals on the obtained porosity and PSD results. However, as shown in our study, it seems to mitigate them significantly.

Therfore, our research contributes to improving the understanding of the potential fluid storage and flow capabilities of different tight rock lithotypes, supporting the exploration of unconventional geothermal and hydrocarbon reservoirs. However, the direct application of the demonstrated differential curves for PSD estimation does not seem to be possible under in-situ studies, as controlling the saturation level of the investigated rock formations in such conditions is insurmountable. Thus, it is impossible to carry out the signal subtraction procedure crucial to the modified differential method. Moreover, using a nonlinear model for PSD estimation, it is impossible to estimate the corresponding surface relaxivity values used in linear models without knowing the proportionality coefficients between pore throats and pore bodies. In such cases, only the T2_clay cut-off values obtained in laboratory measurements using the modified differential protocol can still be employed to determine effective porosity and SWI using the NMR distributions from well-logs more accurately.

The data supporting the conclusions of this study are provided at Mendeley Data repository via https://doi.org/10.17632/8zchvckp4h.4 with a CC BY 4.0 license.

The authors declare that there is no conflict of interest regarding the publication of this paper.

The research leading to these results has received funding from the Norway Grants 2014–2021 via the National Center for Research and Development, Poland (no. NOR/POLNOR/EnerGizerS/0036/2019) and Polish Ministry of Science and Higher Education via AGH University of Krakow statutory funds at the Department of Energy Resources at Faculty of Geology, Geophysics and Environmental Protection (no. 16.16.140.315/05).

We would like to thank W. Mazur-Rosmus for the extensive discussions about the research results and the support provided during the LF-NMR measurements.

Supplementary materials include quantitative X-ray diffractometry results and scatter plots illustrating correlations between para- and ferro-magnetic mineral contents and LF-NMR parameters not included in the main text. In addition, a validation of effective porosity measurements after error correction is included, along with a list of abbreviations used throughout the article.