Analysis of ZTE MRI Application to Sandstone and Carbonate

,


Introduction:
For many years, researchers have employed various experimental techniques to characterize hydrocarbonbearing porous media.One of the aims behind characterizing porous media is to be able to extract hydrocarbons from tight gas and oil reservoirs more efficiently.In this context, the word 'tight' is used to label media with pore size smaller than 100 μm.Efficient hydrocarbon production is aided through the understanding of fundamental transport mechanisms and the bridging among elements of the pore space architecture, namely geometry, connectivity and transport.In tight media, the latter includes not only understanding the diverse pore space, but also possible pore-fracture connectivity.Computed micro-tomography (μCT) and various nuclear magnetic resonance (NMR) methods have been used for this characterization.Magnetic resonance imaging (MRI) has been used for characterization with relative success in recent years.
μCT is one of the most popular techniques used to characterize geological samples.This technique is nondestructive and can resolve features at the micro or even nano scale, making this method advantageous for high resolution.However, μCT utilizes density contrast, making it a complex process to discern pore space from solid materials, which requires rigorous and elaborate modeling, image processing and analysis [1].
Higher-resolution techniques incur in the cost of smaller sample volume.In contrast, MRI offers direct imaging, namely direct spatial encoding, not tomography or inversion protocols, of the three-dimensional (3D) features of porous media.This is possible because MRI draws contrast by encoding the NMR signals to specific nuclei, e.g.proton spins in hydrogen-bearing molecules such as water and magnetic spatial gradients [2].When compared to μCT, MRI also offers the advantage of characterizing a more encompassing image of the pore system by imaging a comparatively larger sample size.The larger size allows resolving multiscale features in a porous system.MRI and μCT can be complimentary to each other through digital registration algorithms [1].NMR methods are not new to geological media characterization, but NMR's evolving versatility offers significant advantages.For instance, Time-Domain NMR (TD-NMR) data can be collected at low-field strength (0.05 T) [1,3,4].With this low field, the induced internal gradients that compromise NMR signals at higher filed strength are minimized.However, TD-NMR spectra are known to have poor signal-to-noise ratios (S/N) and it can be challenging to detect signals with very short (ns) transverse-relaxation times (T 2 ) [1,3].Fluid-flow experiments are also absent in the literature devoted to TD-NMR, because dynamic experiments and the measurement of rheological quantities cannot be performed simultaneously with TD-NMR systems [5].Despite these drawbacks, TD-NMR can be complimentary to MRI methods in that the former can determine volume of water within the sample, T 2 relaxation distributions, and diffusion measurements on a relatively fast timescale.Single-point imaging (SPI) can be performed at low magnetic field, but requires long acquisition times and is often impractical especially for high resolution images [3,6].
High-field NMR instruments, e.g. 7 T or 9.4 T, have been used for imaging in various disciplines, including the geological and medical fields.A commonly used pulse sequence is rapid acquisition with relaxation enhancement with inversion recovery (RARE-Inv-Rec), which uses flip angles of 90°and 180°for an echo recovery [7].This pulse sequence works efficiently with high proton density areas, such as in those found in soft tissue or porous media with pores on the scale of millimeters.Fluids trapped in these samples have long T 2 values (ms), allowing for a long repetition time of the applied pulse sequences [8].With small pores (spanning nm to μm), which are commonly found in bone or geological samples, the T 2 values are much shorter, (tens of μs), and the nuclear spins under study become mostly relaxed before commonly used pulse sequences are fully applied [8].Other pulse sequences must be employed in tight porous media that allow for efficient FID acquisition, along with short repetition times.
Zero Echo Time (ZTE) is a pulse sequence that has been used for image collection of tight porous media.This pulse sequence uses a smaller excitation flip angle (5°), constant gradients and radial encoding techniques [1,3,4,9,10,11].When compared to the more common pulse sequences that use 90°-flip angles, the 5°-flip angle used in ZTE allows for a much shorter pulse time, leading to less relaxation of the magnetization before the receiver is switched on.Another aspect of the traditional imaging pulse sequences that causes them to take longer for full application is the gradient switching process.Additionally, gradient switching can be very taxing on the instrument hardware and cause excessive heating.To avoid gradient switching, ZTE uses nonselective radio frequency (RF) excitation while the gradients are already acting [12,13].The strength of the gradient gradually increases over time, so that spatial encoding is still possible, but the encoding happens with zero delay.The ramping gradients also reduce eddy current effects within the magnet [9].This technique has been used in multiple fields such as imaging tissues with T 2 value on the μs scale in animal bones [2] and human teeth [10,11].The ZTE pulse sequence has also been applied to image porous media such as dolomite, reef reservoirs, and sandstone [1,3,4,11].
Subsequent to the collection of images of a porous medium, pore-cluster analysis (PCA) is frequently used to evaluate properties in geological samples, such as porosity, surface-to-volume ratio, and mass transport such as single-phase (groundwater) or multiphase flow (oil/gas) [15][16][17].Pore clusters are defined by connected pore voxels that are in full contact with each other.The criterion commonly referred in connectivity analysis is that vertex and edges are not considered true connections [15].Pore-cluster analysis allows the examination of pore connectivity through the 3D space and can be used in simulation exercises such as random walks [17].
The primary goal of the work described herein is to test the ability of ZTE, combined with PCA, to resolve both pores and fractures in rock cores with a variety of pore distributions [1,3,4,14], such as sandstone and carbonate lithologies.We initially use a porous medium consisting of 4-mm glass beads submerged in water to test the efficiency and accuracy of a PCA protocol implemented through an in-house developed Wolfram Mathematica 12.0 notebook.We provide evidence that when ZTE combined with PCA is applied to tight porous media, connectivity can be determined, and analyzed producing efficient characterization of tight porous media.Resolution limits are shown to affect the results of this analysis.

Materials:
Solid glass beads were purchased from Fisher Scientific.Berea Sandstone and Indiana Limestone (IL) were purchased from Kocurek Industries Inc. Madison Limestone (ML) was collected from the Rock Springs Uplift (Wyoming) experimental well [18].

Fracturing Procedure :
To analyze fractured rock, select carbonate samples were subject to a fracturing procedure in our laboratory.The process was carried out with a bench vise.The sample was placed horizontally, and the stress was gradually increased until the first crack was observed.Epoxy was immediately applied on the cracks exterior in order to avoid sample disaggregation (crumbling).

Samples:
To test the RARE and ZTE pulse sequences in porous media, various NMR samples were prepared.All samples were loaded in 15 mL Falcon tubes with the dimensions of O.D. 17 mm and length of 120 mm.These tubes were used to reduce the evaporation or loss of water saturation during NMR data collection.
All core plugs were saturated with demineralized water by placing the sample in a vacuum cell for 24 hrs.Sample 1 : 20.866g of 4.0 mm glass beads immersed in 8.0 ml of water.The sample dimensions for the rock cores core are described in Table 1 Fractured Madison Limestone 8

TD-NMR:
In Time-Domain NMR, the relaxation information, whether longitudinal (T 1 ) or transverse relaxation (T 2 ), can be subject to deconvolution to obtain relaxation distribution.T 2 is affected by three underlying relaxation mechanisms: bulk, surface and diffusion relaxation [20,21], as follows: In aqueous solution, the bulk relaxation time usually ranges from 2,000 to 3,000 ms.Its contribution to the harmonic average, therefore, becomes negligible in comparison.Diffusion relaxation, caused by the Brownian movement of the proton spins in an inhomogeneous magnetic field, can be neglected in the presence of small gradients and short echo spacing during measurements [22,23].Therefore, the dominant relaxation within the pore space is surface relaxation, which is a function of the surface-to-volume ratio (S /V ) and surface relaxivity (ρ 2 ).Therefore, T 2 can be approximated in most cases as [23,24]: ) where d is the pore diameter and F s is the pore shape factor.For example, F s for cylindrical and spherical pores equals 2 and 3, respectively.The converted T 2 distribution can be used to estimate the pore-size distribution based on Eq. ( 2).In this study, T 2 measurements for Berea sandstone, Indiana Limestone, and Madison Limestone samples were performed at a resonance frequency of 5.7 MHz and echo spacing of 0.5 ms using a Bruker LF110 instrument at room temperature.The instrument runs a 110-mm diameter, absolute-proton probe with long receiver dead time, but high accuracy.The Carr-Purcell-Meiboom-Grill (CPMG) pulse sequence was used during measurements.The number of scans and the recycle delay in this study were 64 and 4s, respectively.Conversion of relaxation signal into a continuous distribution of relaxation components (T 2 ) was performed using the CONTIN application provided by the equipment vendor.Subsequently, the pore diameter distribution (PDD) is estimated through the transform of T 2 distribution using Eq. ( 2).As surface relaxivity relates to mineralogy and is dependent on the approach used for its determination, this surface property varies among different rock types and even among samples of the same rock type [23,[25][26][27] (Table 2). ) and a matrix of 512 for a resolution of 24.60 x 24.60 x 200 μm 3 was used.The total acquisition time for 32 scan averages, number of scans (NS), was 9 hours and 10 minutes.This method was also tested on the sandstone core, but visualization of the water distribution in the pore space was not successful [5].

ZTE MRI:
Acquisition parameters were: bandwidth = 250,000 Hz, acquisition time = 0.512 ms, repetition time = 2.048 ms, excitation pulse angle = 5.0 degrees, missing data = 1.7 dwells, minimum spooling = 2.0 cycles/pixel, FOV = (2.0 cm) (2.0cm) (2.0 cm) to (4.0 cm) (4.0 cm) (4.0 cm).Dummy scans used were 2,000, resulting in a duration of dummy scans of 4.096 s.With an imaging matrix of 256 x 256 x 256 the achieved nominal resolution was 78.12 x 78.12 x 78.12 μm 3 .Oversampling factor was set to 8. With NS =320, the total acquisition time was 7 hours and 7 minutes per one 3D image.Other combinations of the FOV and matrix were also tested to take account of possible aliasing artifacts.

ZTE Image Analysis:
Principal Component Analysis (PCA) of the 3D ZTE MRI data was performed using Wolfram Research Mathematica Version 12.0.

Image Processing and Analysis:
All images were reconstructed from 2D-slice datasets into 3D volumes, processed and analyzed using a workflow coded in an in-house Mathematica notebook.Figure 2 depicts the workflow developed.This workflow includes: 1. MRI-slice selection: In this step, slices outside the region of interest are visually selected out.This allows maximizing useful information to render 3D volumes.2. 3D rendering: The remaining slides from step 1 are used to render a 3D image.An automated algorithm in Mathematica is used to this end.3. Cropping: Cartesian and cylindrical cropping methods are employed to remove unwanted regions containing elevated levels of noise, particularly from volumes with essentially no signal.
4. Filtering: Sharpening, smoothing and Gaussian filters are applied to prepare the 3D images for binarization and thereby clustering analysis.5. Binarization and segmentation: A binary image is produced through the selection of an adequate threshold value.This step can be conducted automatically through available library functions in Mathematica.However, in this work, we explored the threshold selection process by comparison with alternate sources of data, e.g.gravimetric porosity determination.More on this is discussed in the Results and Discussion section.6. Cluster analysis.In this step, morphological components (different pore systems that are isolated from each other) are clustered and the overall porosity is calculated.The binary image is dependent on choosing a threshold value that effectively represents the pore structure without diluting the signal with low pixel values.When the pore space is saturated with aqueous phase, T 2 values in a pore diameter distribution (PDD) are directly associated with the number of protons within the pore space.In turn, the values of T 2 relate to pore sizes [28] (Eq.( 2)). Figure 3 (A-C) illustrates the T 2 distribution and PDD of three core plugs with two surface relaxivities.The larger surface relaxivity is used to set the upper bound pore diameter to compare with the nominal resolution in MRI analyses in this section.This provides the upper resolution limit needed to visually resolve pores.
As Figure 3A shows how longer T 2 values correspond to larger pore sizes and vice versa.The pore diameter of deionized water-saturated Berea Sandstone estimated by its T 2 distribution exhibits a trimodal distribution ranging from 3 μm to 100 μm (Figure 3A), suggesting that three populations of pore sizes can be observed.
In order, a dominant peak (roughly centered around 400 ms) is associated with the largest pores (30 μm-100 μm).An intermediate peak (centered at roughly 150 ms) represents the mesopores, at a length scale from 10 μm to 30 μm.Finally a smaller peak (centered at roughly 40 ms) resides in the micropore region ranging from 3 μm to 10 μm.The connectivity among pores in a rock sample can be characterized by the connection among multiple peaks in PDD [29].The existence of less defined or somewhat overlapping peaks in the Berea sandstone PDD likely reflects NMR diffusion coupling, which would imply well-developed connections among pores in this rock, most likely through pore-throat connections.
Figure 3B shows the PDD of Indiana Limestone with two bound values of surface relaxivity.In contrast to Berea Sandstone, the PDD of Indiana limestone contains more widespread multiple peaks.The rightmost peak, roughly centered at 1000 ms, represents the macropores, with pore diameters between 40 μm and 200 μm.The leftmost peak, located at approximately 50, indicates the pore diameter occurring in the range from 0.6 μm to 2 μm.The mesopores, ranging from 2 μm to 40 μm, occupy a relatively large portion of the pore volume in the Indiana Limestone.Furthermore, the multiple peaks above 2 μm confirm the high porosity of Indiana Limestone as the amplitude of the PDD peak is directly related to the porosity of the rock samples.The continuous PDD peak connection suggests the multiple well-connected pores of Indiana limestone.
Figure 3C shows a bimodal PDD for the Madison Limestone sample.More specifically, the predominant narrow (right-most) peak, at roughly 350 ms, represents the diameter of the majority of intergranular pores.This PDD is approximately 80.17 μm and a small portion of micropores with diameter less than 3 μm are also present and are shown by the most left peak.Additionally, the disconnected bimodal distribution in Figure 3C suggests a rather disconnected pore system in this sample compared to Berea Sandstone and Indiana Limestone (Figures 3A and 3B).The narrow peaks displayed by this sample suggest the presence of a more homogeneous pore space.It is worth mentioning that vugs (or dissolution pores) cannot be detected in the PDD plots, since the surface relaxation rate in the vugs is small, leading to signals too weak to be detected.Moreover, the porous medium of the Madison Limestone sample is dominated by intergranular porosity, and the vugs are randomly placed and surrounded by intergranular pores.Therefore, the signal from the vugs is further weakened by diffusion coupling between them [23,30].

MRI Glass beads:
The first sample, a 4.0-mm glass bead pack, was used to test the Mathematica notebook workflow described in the materials and methods section.The high connectivity and simplicity of this porous medium makes it ideal for testing the Mathematica notebook for accuracy in porosity calculations and PCA.The larger pore space between individual beads allows for water molecules to yield T 2 values longer than those they would have in a rock sample due to less confinement.These conditions prompted the use of the RARE-Inv-Rec pulse sequence to reproduce the structural features of the simple bead sample.Figure 4 shows the results of image collection with RARE-Inv-Rec (Figure 4A) and the application of the Mathematica notebook to the glass bead sample (workflow steps 1-6, Figure 4B-4E).The set of MRI images was uploaded into the notebook and some slices in the set were selected.The selected 2D slices were used to construct a 3D rendering of the sample.The 3D rendering was then taken through a series of cropping stages including Cartesian and cylindrical cropping (Figure 4B and Figure 4C).Once the 3D rendering was finalized through the cropping stages, a set of filters was applied.These filters include color convert from RBG to greyscale, sharpen, and a Gaussian filter (Figure 4D).As shown by the color scale to the right of the image in Figure 4D, the higher the proton density in the sample, the more yellow the color in the figure will be.This also translates to pixel values between 0 and 1, which are later used in the binarization process.The bright yellow pixels n Figure 4D represent the water that fills the pore space around the beads, whereas the beads appear as empty spaces.Lower proton density can also be observed near the surface of the beads.The high connectivity of the fluid expected from this simple is also clearly detected.Once the filters are applied to the cropped 3D image, the 3D image is then transformed to its binary form (Figure 4E), which allows for further analysis to be applied to the image via PCA.The binary transformation was applied using threshold values between 0.05 and 1.0 in 0.05 increments.This process permits a comprehensive analysis of the sample.Figure 6 shows the analysis of the glass-bead sample at threshold 0.45. Figure 6A presents the applied pore cluster analysis to this sample, resulting in a singular pore system, with all the pores being connected.The histogram in Figure 6B shows that there likely is one main pore system and, with the mode consisting of one voxel probably generated by imaging noise and processing artifacts.The singular pore system surrounds the glass beads.These results indicate that the PCA protocol developed for our studies faithfully reflects the physical reality of the system under consideration.Core Analysis: MRI image sets of Berea sandstone, Indiana Limestone, and Madison Limestone, with pore spaces tighter than other cores previously reported in the literature [1,3,14], were collected using the ZTE pulse sequence.These images were processed and analyzed using the Mathematica notebook described previously.Figure 7 compares the three different core samples in their 3D reconstructions with filters applied.The yellow hue indicates the presence of water in the pore spaces.It is clear from Figure 7 that all three samples are saturated with water, as seen by the bright yellow hue throughout the cores.The Berea sandstone appears more homogeneous with even dispersal of both yellow (water) and red (solid material or space devoid of water) (Figure 7A and 7D).The Indiana Limestone core shows that higher water saturation occurs on the top and bottom of the core (Figure 7B and 7E).The Madison Limestone has larger and more disconnected pore spaces, exemplified by the more concentrated yellow hue wrapping through the core in Figure 7C and 7F.
Figure 7: Side and front views of the 3D reconstruction of (A) and (D) Berea Sandstone, (B) and (E) Indiana Limestone, and (C) and (F) and Madison Limestone.
The Berea Sandstone sample was further analyzed using PCA.From the calculations described earlier, a plot comparing porosity and number of clusters is generated.In the porosity and number of clusters graph as functions of threshold values (Figure 8A), a cluster peak (threshold 0.40) is observed, and no apparent cluster valleys are detected.The cluster peak is observed when the corresponding threshold value diminishes lower voxel values enough to disconnect the sample.At this cluster peak the sample is the most disconnected.However, a calculated porosity of 0.1946, similar to the measured gravimetric value of 0.1800 (Table 3) is found at this threshold value.Ideally, the Berea Sandstone sample should be completely connected, however, it is possible that the resolution achieved with MRI cannot account for the small (tens of microns) throat connections between the pores that dominate this system.This assessment is supported by the TD-NMR results.
In Figure 3(A-C), the black dotted lines represent the resolution of the ZTE MRI technique achieved for image collection.As is clear from this figure, not all the pore-size distributions can be resolved in the collected images, therefore, the connectivity deduced through MRI will not always be a true representation.Based on the TD-NMR results, MRI should allow resolving some of the Berea Sandstone pores, but not the pore throats.If we consider the Berea Sandstone sample exhibiting total connectivity (only one cluster) (Figure 8A), as is true for threshold values 0.05-0.25,then the porosity would be 1.0.However, the gravimetric porosity of this sample is 0.180 (Table 3).Since the calculated porosity is close to the gravimetric porosity a threshold value 0.40, the Berea sample was analyzed at the threshold where the peak is observed.The binary 3D reconstruction of sandstone was transformed to the PCA form (Figure 8B).The PCA separates morphological components based on their connectivity.The sandstone sample appears homogenous with many, relatively, small clusters as seen in Figure 8B.This is supported by the histogram in Figure 8C, where there are many clusters with a range of sizes.Most of the clusters are in the range of 1-10 voxels in size, which is relatively small.Based on these results, Berea Sandstone consists predominately of smaller pores (low pixel count), with pore throats that cannot be resolved with the current MRI method.Although the chosen threshold produces the gravimetric porosity value, it does not produce proper connectivity of the sample, due to resolution limitations, as indicated by the TD-NMR results.Cluster analysis was also performed on the non-fractured Indiana limestone sample.Figure 9 shows a porosity and number of clustervs .threshold plot.The valleys are present at the threshold values of 0.60 and 0.80, indicating where the sample is the most connected.In this analysis, a porosity of 0.1477 was calculated at a threshold value of 0.60.This porosity is similar to the measured gravimetric porosity of 0.173 with a deviation of 14.8% (Table 3).The cluster valley at the threshold value 0.60 (Figure 10A) is indicative of the entire pore system, including small pores, and vugs (dissolution pores).Whereas, the secondary valley (threshold 0.80, Figure 10B) represents higher proton density areas, such as vugs, which do not include the smaller pores.As seen in Figure 10B, the vugs have been isolated, not only from the smaller pores, but also from one another, as the threshold values for the bulk of the small pores have been put outside of the range of interest.These vugs are the connection between groups of smaller pores.Because the vugs have larger water volumes and tend to be space apart through the sample, compared to the small pores and throats in the Berea sample, they fall within the threshold limits and are discernable from one another.Higher connectivity is observed in spite of some of the pores and pore throats being below the resolution limits, as seen by Figure 3B, leading to the conclusion that this sample is dominated by the vug system.Figure 3B also shows that a majority of the porosity is contained in the larger vugs, therefore confirming that this sample is dominated by the vugs.Therefore, this sample is reliably characterized with this workflow, showing relatively higher connectivity when compared to other samples.3).As seen for the non-fractured Indiana Limestone, the first cluster valley (Figure 11) at threshold 0.35, is indicative of small pores and vugs.These vugs connect the small pore systems as shown in Figure 12A.The second valley at threshold of 0.55 (Figure 11) is a depiction of the vug system.When the small pore systems are ignored at a higher threshold value of 0.55, the vugs can be seen distinctly as indicated in Figure 12B.The histograms in Figure 12 show that there are fewer clusters when compared to the Berea sandstone, meaning that this is a more highly connected pore system.There are a greater number of higher proton density distributions in the Madison limestone when compared to the Indiana limestone, which can be correlated to a higher number of larger pores being observed in the Madison limestone when compared to the Indiana limestone.We attribute this result to the structural differences between Indiana and Madison limestone.The imaging and analysis processes previously described were also applied to fractured Indiana and Madison limestone samples, with the goal of observing both fracture and pore systems simultaneously.As previously seen in Figure 9 (non-fractured Indiana limestone), the cluster valley for this sample are observed at threshold values for the non-fractured Indiana limestone are 0.60, and 0.80.The fractured Indiana limestone sample shows similar characteristics to the non-fractured Indiana limestone of small and large vugs as shown in Figure 13.There is a new valley at the threshold 0.30, Figure 14A, because the sample is now more connected than before (the small pores and vugs are more accessible due to the fracture) and intense features in fractures and their vicinity exhibit less signal relaxation.This is exemplified by the threshold value of 0.60, which remains the same for both the fractured and non-fractured Indiana limestone samples, but the number of clusters reduces in the fractured Indiana Limestone (at threshold 0.60) changing from 6 cluster in Figure 10A to 2 clusters in Figure 14B.A cluster valley has also appeared at the threshold value of 0.70, suggesting that there is a new feature that has a proton density that is not a part of the small pore system or the vug system.This new proton density is assigned to the fracture that was manually created in the sample.Figure 14C clearly shows the fracture through the sample (green and gold feature).This fracture has made the vugs more accessible in the Indiana Limestone as shown in Figure 14D (threshold 0.80).There are more vugs present in the fractured sample when compared the non-fractured sample.This result is consistent with the fracture connecting the sample more so than before.Based on this analysis, the pore space in the sample has become more accessible to water, now that a fracture has been introduced.The pores, whether small or vug like, are more connected and therefore could allow better transport through the sample.35 is constant between the non-fractured and the fractured sample, meaning that the high connectivity that was previously observed held constant between the samples.The up shift from threshold 0.55 to 0.60 is significant because there is now more water in the fracture, which appears to be a part of the same threshold range and the vugs.The threshold value 0.60 (Figure 15B) indicates that this sample is more connected than its counterpart (threshold 0.55) in Figure 12B.Evidence of this is provided by the histograms associated with the images, since Figure 12B shows 9 clusters while Figure 15B has 7 clusters.The cluster sizes are also different.The biggest cluster in Figure 15B is significantly larger than the largest cluster in Figure 12B, by a four fold.This can be interpreted to mean that the fracture falls in the same threshold value range as the vugs, and connects not only the entire sample, but also specifically the vugs.This is anticipated, because water is closer to bulk behavior given similar feature sizes of vugs and fractures in this case.Porosities calculated following the workflow and used for the cluster analysis are shown in Table 3.These porosities are comparable to the measured gravimetric porosities.Berea Sandstone has a percent error of 8.1%, while Indiana and Madison limestone have percent errors of 14.8% and 4.2%.Our research shows a practical application of the ZTE MRI pulse sequence in the characterization of sandstone and carbonate samples.The in-house developed PCA is useful, showing a new method of data comparisonvia simultaneous comparison of porosity and the number of clusters present at a given threshold.The glass beads sample confirmed the accuracy of the workflow in an in-house coded Mathematica notebook.This notebook was then applied to other porous samples such as Berea Sandstone and various limestone samples.Connectivity patterns throughout the samples were deduced from the PCA and could correlate various features in the sample with specific threshold values.This leads to an in depth comparison between non-fractured and fractured samples, where each feature could be separately visualized and analyzed.Unfortunately, resolution limitations do not allow for full connectivity characterization as many small pores and pore throats are below the resolution achieved through MRI, as confirmed by TD-NMR NMR measurements.In general, the PDD obtained from the T 2 distribution can be considered a complement to MRI in pore structure characterization.This is because the resolution limitations of MRI do now allow visualization of the smaller pore sizes.In Figure 3A-3C, the black dotted line indicates the resolution of the MRI technique presented herein.From the results, it becomes clear that not all the pore-size distributions can be resolved by the MRI pulse sequence used in this work, therefore the connectivity deduced will not always be a complete representation of the connected pore space.Results also indicate that some of the Berea Sandstone pores should be visible in the MRI images, but not the pore throats.Additionally, most of the large pores in the Indiana Limestone should also be discernable, with a larger volume of the small pores missing.Lastly, the Madison Limestone sample shows that the bulk of its porosity lies in the large pores, which should be resolved with MRI.TD-NMR confirms that the samples are saturated and have a large range of pore sizes and that MRI is resolving pores as opposed to noise or artifacts.
In view of the connection between the signal intensity and feature, threshold choices represent scale sweeps, namely the threshold selection serves to identify different pore groups in the samples.

Figure 2 .
Figure 2. Mathematica Notebook workflow Results and Discussion: TD-NMR Pore Diameter Distributions:

Figure 4 .
Figure 4. Workflow example.(A) Slice selector; (B) Cartesian cropping; (C) cylindrical cropping; (D) filtered image; (E) and binary form.With the aim of selecting the optimum threshold value for the glass-bead sample, various calculations were performed on the 3D rendered image, and information regarding porosity, number of clusters, and size of clusters was derived.These calculations were carried out for each of the considered thresholds and are plotted in Figure5.For this sample, a valley in the number of clusters vs threshold plot is observed between threshold values 0.40 and 0.45.The porosity of the glass beads calculated with Mathematica was found to be 0.378, at threshold value 0.45 located in the valley.On the other hand, Nakashima et al .determined a gravimetric porosity of 0.370 for a glass-bead sample[16].These results indicate that a threshold value of 0.45 should be optimum for the glass-bead sample.A glass-bead sample is known to have only one pore system, and therefore should consist of a low number of clusters in the PCA (the valley).The results of this experiment indicate that the glass-bead sample can be successfully used as a model to test the Mathematica code.

Figure 5 .
Figure 5. Large Beads Porosity and Number of Clusters vs.threshold value

Figure 6 :
Figure 6: PCA of the glass bead sample.(A) PCA image; (B) Histogram of pore size distribution; and (C) photograph of the glass beads in a 15 mL falcon tube, where the blue line indicates the FOV.

Figure 8 .
Figure 8. PCA analysis of Berea Sandstone.(A) Porosity and number of cluster graph vs. threshold value plot, (B) PCA image, (C) Histogram of cluster-size distribution, and (D) photograph of sample (blue line indicates the FOV).

Figure 9 :
Figure 9: Indiana Limestone.(A) Porosity and number of clusters vs. threshold, and (B) photograph of sample (blue line indicates the FOV).

Figure 10 :
Figure 10: Indiana limestone sample cluster analysis at different threshold values: (A) 0.60, and (B) 0.80 that correspond to the valleys found in Figure 9.

Figure 11
Figure11shows the PCA analysis of the non-fractured Madison Limestone sample.The valleys are observed at threshold values of 0.35 and 0.55.The gravimetric porosity for this sample is close to the MRI porosity calculated at a threshold value of 0.55 (Table3).As seen for the non-fractured Indiana Limestone, the first cluster valley (Figure11) at threshold 0.35, is indicative of small pores and vugs.These vugs connect the small pore systems as shown in Figure12A.The second valley at threshold of 0.55 (Figure11) is a depiction of the vug system.When the small pore systems are ignored at a higher threshold value of 0.55, the vugs can be seen distinctly as indicated in Figure12B.The histograms in Figure12show that there are fewer clusters when compared to the Berea sandstone, meaning that this is a more highly connected pore system.There are a greater number of higher proton density distributions in the Madison limestone when compared to the Indiana limestone, which can be correlated to a higher number of larger pores being observed in the Madison limestone when compared to the Indiana limestone.We attribute this result to the structural differences between Indiana and Madison limestone.

Figure 11 :
Figure 11: Madison Limestone.(A) Porosity and number of clusters vs. threshold, and (B) photograph of the sample where the blue line indicates the FOV.

Figure 12 :
Figure 12: Madison limestone cluster analysis at the different threshold values: (A) 0.35, and (B) 0.55 that correspond to cluster valleys in Figure11.

Figure 13 :
Figure 13: Fractured Limestone.Porosity and Number of Clustersvs .threshold for: (A) Indiana Limestone, and (B) Madison Limestone.(C) photograph of the fractured Indiana limestone where the blue line is the FOV.

Figure 14 :
Figure 14: Fractured Indiana Limestone cluster analysis at threshold values: (A) 0.30, (B) 0.60, (C) 0.70, and (D) 0.80 that correspond to cluster valleys found in Figure 14.The fractured Madison Limestone does not behave significantly different from its counterpart shown in Figure 12.The non-fractured Madison Limestone sample has cluster valley thresholds of 0.35, and 0.55, compared to the fractured sample with cluster valley thresholds at 0.35 and 0.60 in Figure 13B.The threshold value

Figure 15 :
Figure 15: Fractured Madison Limestone cluster analysis at threshold values: (A) 0.35 (A), and (B) 0.60 that correspond to the cluster valleys found in Figure 14.

Table 1 .
. Description of the rock core sample dimensions Table1.Description of the rock core sample dimensions Table

Table 2 .
Surface relaxivities for sandstone and limestone samples Table2.Surface relaxivities for sandstone and limestone MRI was carried out using a Bruker Avance III 300 MHz wide-bore NMR (8.9 cm) with ParaVision 5.1 software, used for acquisition and reconstruction.The probe used is the Bruker Micro 2.5 MICWB40-1H 040/025 quadrature transmission radio frequency coils that can generate a gradient strength up to 1.5 T/m at 60 A.
RARE MRI:Acquisition parameters were: echo time of 12 ms, effective echo time of 36 ms, rare factor of 8, repetition time of 3665.46 ms, effective spectra bandwidth of 37,878.8Hz, duration of read diphase of 1.3 ms, duration of 2D phase gradient of 1 ms, repetition spoiler duration of 2.0 ms with 30.0%strength, excitation length of 2.7 ms with 2,000 Hz bandwidth, refocusing pulse (180 degrees) length= 1.71 ms with 2,000 Hz bandwidth, number of dummy scans = 2.A field of view (FOV) of 1.26 cm x 1.26 cm with 80 slices (slice thickness = 200 μm

Table 3 .
Comparison porosities of the cluster calculations performed on each sample.