\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[greek,ngerman,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Analysis of ZTE MRI Application to Sandstone and Carbonate}



\author[1]{Kyle Covington}%
\author[1]{Alexander Goroncy}%
\author[1]{Teresa Lehmann}%
\author[2]{Zuhao Kou}%
\author[3]{Heng Wang}%
\author[2,3]{Vladimir Alvarado}%
\affil[1]{University of Wyoming College of Arts and Sciences}%
\affil[2]{University of Wyoming College of Engineering and Applied Science}%
\affil[3]{University of Wyoming}%


\vspace{-1em}



  
  \date{June 16, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Techniques such as micro-tomography (\selectlanguage{greek}µ\selectlanguage{english}CT) and nuclear magnetic resonance
(NMR) methods have been used to characterize porous media for decades.
Magnetic Resonance Imaging (MRI) offers direct visualization of pore
architecture through a vast number of pulse sequences. In this work, we
tested the MRI pulse sequence Zero Echo Time (ZTE) in the study of
sandstone and carbonate. ZTE has been used previously to image highly
porous geological material with relative success. In this research, ZTE
was used to study porous systems tighter than those previously reported.
We show that ZTE can resolve not only pore systems, but also a
combination of fractures and pores simultaneously. Additionally, Pore
Cluster Analysis (PCA), combined with ZTE, can be used to analyze
connectivity of relatively large volume core plugs. By combining
Time-Domain NMR pore-size analysis and PCA, we show that careful
selection of resolution is necessary to understand transport in porous
media.%
\end{abstract}%



\sloppy


\textbf{Introduction:}

For many years, researchers have employed various experimental
techniques to characterize hydrocarbon-bearing 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 \selectlanguage{greek}µ\selectlanguage{english}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 (\selectlanguage{greek}µ\selectlanguage{english}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.

\selectlanguage{greek}µ\selectlanguage{english}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, \selectlanguage{greek}µ\selectlanguage{english}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 \selectlanguage{greek}µ\selectlanguage{english}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 \selectlanguage{greek}µ\selectlanguage{english}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\textsubscript{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\textsubscript{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\selectlanguage{ngerman}° 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\textsubscript{2} values (ms), allowing for a long repetition time of
the applied pulse sequences {[}8{]}. With small pores (spanning nm to
\selectlanguage{greek}µ\selectlanguage{english}m), which are commonly found in bone or geological samples, the
T\textsubscript{2} values are much shorter, (tens of \selectlanguage{greek}µ\selectlanguage{english}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\selectlanguage{ngerman}°), 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\textsubscript{2} value
on the \selectlanguage{greek}µ\selectlanguage{english}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-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.

\textbf{Materials and Methods:}

\emph{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{]}.

\emph{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).

\emph{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. \textbf{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.\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
Table 1. Description of the rock core sample dimensions & Table 1.
Description of the rock core sample dimensions & Table 1. Description of
the rock core sample dimensions & Table 1. Description of the rock core
sample dimensions &\tabularnewline
\midrule
\endhead
\textbf{Sample} & \textbf{Description} & \textbf{Outer Diameter (mm)} &
\textbf{Length (mm)} & \textbf{Length (mm)}\tabularnewline
2 & Berea Sandstone & 15 & 70 & 70\tabularnewline
3 & Indiana Limestone & 8 & 43 & 43\tabularnewline
4 & Madison Limestone & 8 & 43 & 43\tabularnewline
5 & Fractured Indiana Limestone & 8 & 23 & 23\tabularnewline
6 & Fractured Madison Limestone & 8 & 23 & 23\tabularnewline
\bottomrule
\end{longtable}

\emph{TD-NMR:}

In Time-Domain NMR, the relaxation information, whether longitudinal
(T\textsubscript{1}) or transverse relaxation (T\textsubscript{2}), can
be subject to deconvolution to obtain relaxation distribution.
T\textsubscript{2} is affected by three underlying relaxation
mechanisms: bulk, surface and diffusion relaxation {[}20,21{]}, as
follows:

\(\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2bulk}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2surface}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2diffusion}}}\)(1)

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 (\emph{S} /\emph{V} ) and
surface relaxivity (\(\mathrm{\rho}_{\mathrm{2}})\). Therefore, T\textsubscript{2}
can be approximated in most cases as {[}23,24{]}:

\(\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2}}}\mathrm{\approx}\frac{\mathrm{1}}{\mathrm{T}_{\mathrm{2surface}}}\mathrm{=}\mathrm{\rho}_{\mathrm{2}}\frac{\mathrm{S}}{\mathrm{V}}\mathrm{=2}\mathrm{\rho}_{\mathrm{2}}\frac{\mathrm{F}_{\mathrm{s}}}{\mathrm{d}}\)(2)

where \emph{d} is the pore diameter and \emph{F\textsubscript{s}} is the
pore shape factor. For example, \emph{F\textsubscript{s}} for
cylindrical and spherical pores equals 2 and 3, respectively. The
converted T\textsubscript{2} distribution can be used to estimate the
pore-size distribution based on Eq. (2). In this study,
T\textsubscript{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\textsubscript{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\textsubscript{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-27{]} (Table 2).\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Table 2. Surface relaxivities for sandstone and limestone samples &
Table 2. Surface relaxivities for sandstone and limestone samples &
Table 2. Surface relaxivities for sandstone and limestone samples &
Table 2. Surface relaxivities for sandstone and limestone
samples\tabularnewline
\midrule
\endhead
Parameter & Berea sandstone & Indiana limestone & Madison
Limestone\tabularnewline
\(\mathrm{\rho}_{\mathrm{2}}\) (\selectlanguage{greek}μ\selectlanguage{english}m/s) & 16-26 & 2.6-16 & 5-38.7\tabularnewline
\bottomrule
\end{longtable}

\emph{MRI:}

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.

\emph{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.8 Hz, 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 \selectlanguage{greek}µ\selectlanguage{english}m) and a matrix of 512
for a resolution of 24.60 x 24.60 x 200 \selectlanguage{greek}µ\selectlanguage{english}m\textsuperscript{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{]}.

\emph{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 \selectlanguage{greek}µ\selectlanguage{english}m\textsuperscript{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.

\emph{ZTE Image Analysis:}

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

\emph{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:

\begin{enumerate}
\tightlist
\item
  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.
\item
  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.
\item
  Cropping: Cartesian and cylindrical cropping methods are employed to
  remove unwanted regions containing elevated levels of noise,
  particularly from volumes with essentially no signal.
\item
  Filtering: Sharpening, smoothing and Gaussian filters are applied to
  prepare the 3D images for binarization and thereby clustering
  analysis.
\item
  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.
\item
  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.
\end{enumerate}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 2. Mathematica Notebook workflow

\textbf{Results and Discussion:}

\emph{TD-NMR}

\emph{Pore Diameter Distributions:}

When the pore space is saturated with aqueous phase,
T\textsubscript{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\textsubscript{2} relate to pore sizes {[}28{]}
(Eq. (2)). Figure 3 (A-C) illustrates the T\textsubscript{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\textsubscript{2} values correspond to
larger pore sizes and vice versa. The pore diameter of deionized
water-saturated Berea Sandstone estimated by its
T\textsubscript{2}distribution exhibits a trimodal distribution ranging
from 3 \selectlanguage{greek}μ\selectlanguage{english}m to 100 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}m-100 \selectlanguage{greek}μ\selectlanguage{english}m). An
intermediate peak (centered at roughly 150 ms) represents the mesopores,
at a length scale from 10 \selectlanguage{greek}μ\selectlanguage{english}m to 30 \selectlanguage{greek}μ\selectlanguage{english}m. Finally a smaller peak (centered
at roughly 40 ms) resides in the micropore region ranging from 3 \selectlanguage{greek}μ\selectlanguage{english}m to
10 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}m and 200 \selectlanguage{greek}μ\selectlanguage{english}m. The leftmost peak, located at
approximately 50, indicates the pore diameter occurring in the range
from 0.6 \selectlanguage{greek}μ\selectlanguage{english}m to 2 \selectlanguage{greek}μ\selectlanguage{english}m. The mesopores, ranging from 2 \selectlanguage{greek}μ\selectlanguage{english}m to 40 \selectlanguage{greek}μ\selectlanguage{english}m, occupy a
relatively large portion of the pore volume in the Indiana Limestone.
Furthermore, the multiple peaks above 2 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}m and a small portion of micropores with
diameter less than 3 \selectlanguage{greek}μ\selectlanguage{english}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{]}.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

Figure 3. T\textsubscript{2} distribution and PDD analysis of (A) Berea
sandstone, (B) Indiana limestone, and (C) Madison Limestone sample.

\emph{MRI}

\emph{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\textsubscript{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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

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 Figure 5. For this sample,
a valley in the number of clusters \emph{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 \emph{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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Figure 5. Large Beads Porosity and Number of Clusters
\emph{vs.}threshold value

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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

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.

\emph{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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

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

Cluster analysis was also performed on the non-fractured Indiana
limestone sample. Figure 9 shows a porosity and number of
cluster\emph{vs} . 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

Figure 9: Indiana Limestone. (A) Porosity and number of clusters vs.
threshold, and (B) photograph of sample (blue line indicates the FOV).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

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 shows 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 (Table 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}

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.

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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}

Figure 13: Fractured Limestone. Porosity and Number of Clusters\emph{vs}
. threshold for: (A) Indiana Limestone, and (B) Madison Limestone. (C)
photograph of the fractured Indiana limestone where the blue line is the
FOV.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image14/image14}

\end{center}
\end{figure}

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 of 0.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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image15/image15}

\end{center}
\end{figure}

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.

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\%.

Table 3. Comparison porosities of the cluster calculations performed on
each sample.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\textbf{Core} & \textbf{MRI Porosity} & \textbf{Gravimetric
Porosity}\tabularnewline
\midrule
\endhead
\textbf{Berea} & 0.1946 & 0.1800\tabularnewline
\textbf{Indiana} & 0.1477 & 0.1731\tabularnewline
\textbf{Madison} & 0.1245 & 0.1300\tabularnewline
\textbf{Indiana Fracture} & 0.1506 & -\tabularnewline
\textbf{Madison Fracture} & 0.1813 & -\tabularnewline
\bottomrule
\end{longtable}

Conclusions:

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
comparison\emph{via} 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\textsubscript{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.

Acknowledgements:

This work was supported as part of the Center for Mechanistic Control of
Water-Hydrocarbon-Rock Interactions in Unconventional and Tight Oil
Formations (CMC-UF), an Energy Frontier Research Center funded by the
U.S. Department of Energy, Office of Science under DOE (BES) Award
DE-SC0019165.

References:

\begin{enumerate}
\tightlist
\item
  Fheed A, Krzy\selectlanguage{polish}ż\selectlanguage{english}ak A, \selectlanguage{polish}Ś\selectlanguage{english}wierczewska A. Exploring a carbonate reef
  reservoir -- nuclear magnetic resonance and computed microtomography
  confronted with narrow channel and fracture porosity. \emph{J Appl
  Geophys} . 2018;151:343-358.
\item
  Giovannetti G, Guerrini A, Salvadori, P. Magnetic Resonance
  Spectroscopy and Imaging for the Study of Fossils. \emph{Mag Reson
  Imaging} . 2016;34:730-742.
\item
  Weglarz WP, Krzyzak A, Stefaniuk M. ZTE imaging of tight sandstone
  rocks at 9.4T-Comparison with standard NMR analysis at
  0.05T.\emph{Magn Reson Imaging} . 2016;34:492--495.
\item
  Weglarz WP, Krzy\selectlanguage{polish}ż\selectlanguage{english}ak A, Machowski G, Stefaniuk M. ZTE MRI in high
  magnetic field as a time effective 3D imaging technique for monitoring
  water ingress in porous rocks at sub-millimetre resolution. \emph{Mag
  Reson Imaging} . 2018; 47: 54-59.
\item
  Ratzsch KF, Friedrich C, Wilhelm M. Low-field rheo-NMR: A novel
  combination of NMR relaxometry with high-end shear rheology. \emph{J
  Rheology} . 2017;61:905-917.
\item
  Beyea SD, Balcom BJ, Prado PJ, Cross AR, Kennedy CB, Armstrong RL,
  Bremner TW. Relaxation time mapping of short T\textsubscript{2} nuclei
  with single-point imaging (SPI) methods. \emph{J Magn Reson} .
  1998;135:156--64.
\item
  Rinck P. \emph{Magnetic Resonance in Medicine} . \emph{The Basic
  Textbook of the European Magnetic Resonance
  Forum}(11\textsuperscript{th} edition). Berlin: Blackwell Science,
  2017.
\item
  Grunewald EB, Knight R. The effect of pore size and magnetic
  susceptibility on the surface NMR relaxation parameter
  T\textsubscript{2}*. \emph{Near Surf Geophys} . 2011;9:169-178.
\item
  Weiger M, Pruessmann KP. MRI with Zero Echo Time. \emph{eMagRes} .
  2012;1:311--322.
\item
  Weiger M, Pruessmann KP, Bracher AK, Kohler S, Lehmann V, Wolfram U,
  Hennel F, Rasche V. High-resolution ZTE imaging of human
  teeth.\emph{NMR Biomed} . 2012;25:1144-1151.
\item
  Weiger M, Brunner DO, Dietrich BE, M\selectlanguage{ngerman}üller CF, Pruessmann KP. ZTE
  imaging in humans. \emph{Mag Reson Med} . 2013;70:328-332.
\item
  Schieban K, Weiger M, Hennel F, Boss A, Pruessmann KP. ZTE imaging
  with enhanced flip angle using modulated excitation. \emph{Mag Reson
  Med} . 2015;74:684-693.
\item
  Weiger M, Brunner DO, Tabbert M, Pavan M, Schmid T, Pruessmann KP.
  Exploring the bandwidth limits of ZTE imaging: Spatial response, out
  of band signals, and noise propagation. \emph{Mag Reson Med} .
  2015;74:1236-1247.
\item
  Nikolaev MY, Kazak AV. Liquid saturation evaluation in organic-rich
  unconventional reservoirs: A comprehensive review. \emph{Earth Sci
  Reviews} . 2019;194:327-349.
\item
  Nakashima Y, Watanabe Y.~Estimate of transport properties of porous
  media by microfocus X-ray computed tomography and random walk
  simulation.~\emph{Water Resour Res} . 2002;38:1272.
\item
  Nakashima Y, Nakano T, Nakamura K, Uesugi K, Tsuchiyama A, Ikeda S.
  Three-dimensional diffusion of non-sorbing species in porous
  sandstone: computer simulation based on X-ray microtomography using
  synchrotron radiation. \emph{J Contam Hydrol} . 2004;74: 253-264.
\item
  Nakashima Y, Kamiya S. Mathematica Programs for the Analysis of
  Three-Dimensional Pore Connectivity and Anisotropic Tortuosity of
  Porous Rocks using X-ray Computed Tomography Image Data. \emph{J Nucl
  Sci Technology} . 2007;44:1233-1247.
\item
  Wang H, Alvarado V, McLaughlin JF, Bagdonas DA, Kaszuba JP, Campbell
  E, Grana D. Low-field nuclear magnetic resonance characterization of
  carbonate and sandstone reservoirs from Rock Spring Uplift of
  Wyoming.\emph{J Geophys Res: Solid Earth} . 2018;123:7444--7460.
\item
  Wang H, Alvarado V. Ionic strength-dependent pre-asymptoic diffusion
  coefficient distribution in porous media-determination through the
  pulsed field gradient technique. \emph{J Nat Gas Science and
  Engineering} . 2018;49:250--9.
\item
  Zhao PQ, Wang ZL, Sun ZC, Cai JC, Wang L. Investigation on the pore
  structure and multifractal characteristics of tight oil reservoirs
  using NMR measurements: Permian Lucaogou Formation in Jimusaer Sag,
  Junggar Basin. \emph{Mar Petrol Geol} . 2017;86:1067-1081.
\item
  Xu HJ, Fan YR, Hu FL, Li CX, Yu J, Liu ZC, Wang FY. Characterization
  of pore throat size distribution in tight sandstones with Nuclear
  Magnetic Resonance and High-Pressure Mercury Intrusion.\emph{Energies}
  . 2019;12:1528.
\item
  Daigle H, Thomas B, Rowe H, Nieto, M. Nuclear magnetic resonance
  characterization of shallow marine sediments from the Nankai Trough,
  integrated ocean drilling program expedition 333. \emph{J Geophys Res:
  Solid Earth} . 2014;119:2631--2650.
\item
  Wang H, Alvarado V, McLaughlin JF, Bagdonas DA, Kaszuba JP, Campbell
  E, Grana D. Low-field nuclear magnetic resonance characterization of
  carbonate and sandstone reservoirs from Rock Spring Uplift of
  Wyoming.\emph{J Geophys Res: Solid Earth} . 2018;123:7444--7460.
\item
  Lyu CH, Ning ZF, Wang Q, Chen MQ. Application of NMR T2 to pore size
  distribution and movable fluid distribution in tight
  sandstones.\emph{Energ Fuel} . 2018;32:1395-1405.
\item
  Horlimann MD. Restricted diffusion in sedimentary rocks: Determination
  of surface-area-to-volume ratio and surface relaxivity. \emph{J Mag
  Reson A} . 1994;111:169-178.
\item
  Fleury M. NMR surface relaxivity determination using NMR apparent
  diffusion curves and BET measurements. Paper presented at the
  International Symposium of the Society of Core Analysts, Calgary,
  Canada, September10-12, 2007.
\item
  Carneiro G, Souza A, Boyd A, et al. Evaluating pore space connectivity
  by NMR diffusive coupling. Paper presented at the SPWLA 55th Annual
  Logging Symposium, Abu Dhabi, United Arab Emirates, May 18-22, 2014.
\item
  Lai J, Wang GW, Wang ZY, Chen J, Pang XJ, Wang SC, Zhou ZL, He ZB, Qin
  ZQ, Fan XQ. A review on pore structure characterization in tight
  sandstones. \emph{Earth Sci Review} . 2018;177:436-457.
\item
  Yao YB, Liu DM, Che Y, Tang DZ, Tang SH, Huang WH. Petrophysical
  characterization of coals by low-field nuclear magnetic resonance
  (NMR). \emph{Fuel} . 2010;89:1371-1380.
\item
  Godefroy S, Korb JP, Fleury M, Bryant RG. Surface nuclear magnetic
  relaxation and dynamics of water and oil in macroporous
  media.\emph{Phys Rev E} . 2001;64:021605.
\end{enumerate}

\selectlanguage{english}
\FloatBarrier
\end{document}

