Freely available online through the author-supported open access option.

Abstract

Nuclear magnetic resonance (NMR) spectroscopy has become an increasingly important technique in recent years for the characterization of porous materials. In particular, its importance in petroleum exploration has been enhanced by new NMR well-logging techniques for the characterization of both rocks and natural fluids. Such advanced techniques have been accepted as one of the most valuable logging services. There are several areas currently being developed in the study of pore structure and rock heterogeneity including the physics of internal magnetic fields, the two-dimensional NMR paradigm, heterogeneity decomposition, and the mathematics of Laplace inversion. Although the methods have been developed for rocks, they are generally applicable to partially saturated porous media such as the vadose zone.

Nuclear magnetic resonance has become an increasingly important technique for characterization of porous materials, in particular for petroleum exploration. Such advanced techniques have been accepted as one of the most valuable logging services. This development, although originally for rocks, could be applicable to partially saturated formations such as the vadose zone.

Nuclear magnetic resonance spectroscopy is considered one of the most powerful analytical methods for determining the molecular structure, chemical composition, and physical properties of a wide range of materials including liquids, solids, and biological tissues. Such techniques, however, have been primarily limited to laboratory use and relatively pure materials until recently. With the development of inside-out NMR systems such as a well-logging tool (Kleinberg, 1995) and NMR-MOUSE (Eidmann et al., 1996), and applications in petroleum exploration, the use of NMR to characterize porous materials has been growing rapidly. Because porous media are often filled with fluids (e.g., water) in their natural state, such as underground rock formations, sand packs, and cements, measurements of the spatial structure of the fluid distribution sheds light on the solid granular structure and thus the properties of such materials. Recent NMR development has demonstrated the use of relaxation, diffusion, and imaging techniques to quantify pore structure, fluid properties, and rock heterogeneity. There has been accelerated development and application of NMR techniques in many areas of research and commercial applications for porous media. Here, I outline a few aspects of the technical development trends in which I am most involved and discuss their potential applications in porous media such as soils. I focus on the essential physics involved in the new techniques, which is generally applicable. This is not, on the other hand, a thorough review of all aspects of the field nor a compendium of all relevant references. A more thorough and in-depth review would be very valuable.

Background: Relaxation, Diffusion, Internal Field

When fluids saturate a porous rock, the solid granular structure is reflected in the spatial distribution of the fluids. Furthermore, the fluid–rock interactions at the granular surfaces are important for the rock properties and the fluid flow inside the rock. In natural environments, rocks are often filled with water, which provides strong NMR signals. Also, the NMR property of water is affected by the rock–fluid interactions including surface relaxation, restricted diffusion, and dispersion. These effects can be readily measured by modern NMR technology.

The idea of using spin relaxation as a probe of pore space was first proposed by Senturia and Robinson (1970). In the early 1980s, it was understood that spin relaxation can be used to determine the surface/volume ratio (Cohen and Mendelson, 1982; de Gennes, 1982). It has been shown that different relaxation rates reflect different pools of water in the microstructure of cement and rocks (Kleinberg, 1995; Halperin et al., 1989). Since the early 1980s, NMR has been used in commercial well-logging services for oilfield exploration. In the early 1990s, pulsed NMR was introduced into the oil and gas industry and NMR has since become increasingly important for the evaluation of subsurface earth formations.

For the fluids inside the pore space, molecular diffusion is a key element linking the NMR behavior of the fluid to the properties of the porous material. For simple pore geometries, such as a spherical pore, the theory of molecular diffusion is well understood. The differential equation governing the evolution of the longitudinal or the transverse components of nuclear spin magnetization carried by diffusing molecules is the well-known Torrey–Bloch equation (Torrey, 1956): 
\[\frac{{\partial}}{{\partial}t}m(r,t){=}D\mathrm{{\nabla}}^{2}m(r,t){-}\mathrm{{\mu}}m(r,t)\]
[1]
where D is the bulk diffusion constant, μ is the bulk spin relaxation rate, m is the magnetization deviation from its equilibrium, r is a spatial coordinate, and t is time. The detected signal m(t) is an integral of m over the detected sample volume. A solution can be found in the generic form: 
\[m(r,t){=}{{\sum}_{n{=}0}^{{\infty}}}a_{n}\mathrm{{\phi}}_{n(r)}\mathrm{exp}\left(\frac{{-}t}{T_{n}}\right)\]
[2]
where ϕn(r) and Tn are eigenfunctions and eigenvalues of the Helmholtz equation. For simple pore geometries, such eigenfunctions can be calculated analytically (Brownstein and Tarr, 1979). For a complex geometry, numerical solutions can be obtained (de Hoop and Prange, 2007). These equations are the basis for understanding many relaxation and diffusion phenomena of NMR in porous media. Understanding of the eigenfunctions is very important for application to complex pore geometry such as interacting pores as an example of the complex diffusion dynamics (Song et al., 2008b).

In addition to relaxation and diffusion behavior, another important feature of fluids in porous media is the presence of a significantly inhomogeneous magnetic field inside the porous structure. This is due to the magnetic susceptibility difference between the solid grains and the interstitial fluids. When subject to a strong magnetic field such as an NMR magnet, the grains are magnetized. These magnetized grains in turn produce a local magnetic field inside the pore space. It has been shown that such a magnetic field exhibits a correlation length scale that is closely related to the structural correlation length reflected in the structure factor (Audoly et al., 2003). The effects of such an internal field can be incorporated into the diffusion equations (Eq. [1] and [2]) and has been used to determine pore structure (Song et al., 2000).

Internal Field Correlation

Nuclear magnetic resonance diffusion measurements are usually performed by measuring the decoherence of the diffusing spins' magnetization under the influence of applied magnetic field gradients. Any medium, however, for which a difference in magnetic susceptibility Δχ exists between the solid and liquid components will exhibit a spatially varying internal magnetic field when placed in a homogeneous external field B0. Such background gradients can be estimated as ΔχB0/d, where d is the pore size. For instance, this background gradient has been estimated to be up a few hundred gauss per centimeter at high magnetic field strengths (Hürlimann, 1998). Several techniques have been developed to reduce the effect of such internal field gradients on measurements of the diffusion constant (Cotts et al., 1989; Sun et al., 2003; Galvosas et al., 2004). Such internal fields often cause artifacts and a loss of signal, in particular for imaging experiments (e.g., Hall et al., 1997; Votrubová et al., 2000).

While it is not entirely surprising that the internal field can be a fingerprint of the underlying pore geometry, in fact, quantitative connections exist between the pore structure and the magnetic field structure. For example, an interesting scaling relationship was reported for NMR at different scales (Borgia et al., 1995). A numerical study of the internal magnetic field within the pore space of a random, dense pack of magnetized spheres (Audoly et al., 2003) has demonstrated the near equivalence of two paired correlation functions: that of the mass distribution of the grains, and that of the internal magnetic field in the pore space. Intuitively, because the dipolar magnetic field decays monotonically, the pore size is the only length that determines the magnetic field variation due to susceptibility difference.

Recently, the NMR technique (Cho and Song, 2008) has been developed to directly measure the two-point field correlation function in a randomly packed sphere samples. The qualitative idea is the following. To detect the field correlation function, we used nuclear spins to probe the magnetic field difference between two positions by the following method. First, a pulsed field gradient NMR method was used to select spins by their translational diffusion displacement. For this experiment, the effect of the internal field was nullified. Then, a similar experiment was performed but with the internal field effect activated through pulse sequence design. The field difference between the two positions causes a signal decay that is related to the magnetic field correlation at that displacement. The method was found to be robust and measurements on different grain sizes, diffusion, encoding times, and displacement resolutions were highly consistent. This technique represents a novel method to characterize the structure factor of a porous material, which is a fundamental structural property, as it is the key for solid crystals via measurements of x-ray and neutron scattering.

For porous samples with larger grains, it is possible to use magnetic resonance imaging (MRI) to directly image the spatial distribution of grains and fluid. A technique was developed to measure the internal field gradient structure in such samples through direct imaging (Cho et al., 2009). The technique uses a conventional stimulated echo with three 90° radio frequency pulses and a spin-warp imaging is performed at the end of the sequence. The imaging is performed as the delay time (td) between the second and third pulses are varied. The signal as a function of td shows a decay behavior similar to spin-lattice relaxation, T1. Because of the presence of the internal field gradients (gint), however, the actual decay time constant (Teff) contains a contribution from gint: 
\[\frac{1}{T_{\mathrm{eff}}}{=}\frac{1}{T_{1}}{+}\mathrm{{\gamma}}^{2}t_{e}^{2}Dg_{\mathrm{int}}^{2}\]
[3]
where γ is the gyromagnetic ratio, and te is the echo time. Since T1 can be measured separately, this technique can obtain an image of gint2. The measurement on a capillary pack showed excellent comparison with theoretical calculated internal field gradients.
In addition, Cho et al. (2009) also used this method to image the direction of the field gradient. They applied a weak external gradient during the encoding period so that the total gradient (g) was the sum of the applied (gext) and the internal gradient: g = gext + gint. This technique directly measures the squares of the total gradient, thus 
\[g^{2}{=}g_{\mathrm{int}}^{2}{+}g_{\mathrm{int}}^{2}{+}g_{\mathrm{int}}{\cdot}g_{\mathrm{ext}}\]
[4]
Because the gradient is a vector, the dot here signifies the vector dot product. Thus, by performing experiments with different applied gradient directions, Cho et al. (2009) was able to obtain a full image of the gradient direction.

For three-dimensional porous media such as rocks, extensive work has been reported on the measurement of the pore size distribution using the internal field effect. The internal field correlation causes an effective decay and the decay rate is closely related to the diffusion time across the pores. As a result, a measurement of the decay rate distribution is equivalent to the pore size distribution. This method is called decay due to diffusion in internal field (DDIF) and it has been applied to both sandstones and carbonate rocks. For sandstones with a relatively uniform grain size distribution, the DDIF pore size distribution shows a unimodal shape (Fig. 1 ) consistent with a simple pore size distribution.

Mercury porosimetry is commonly used to characterize pore throat size in a rock, as shown in Fig. 1. On the other hand, DDIF measures both the pore body and the pore throat. The overlay of the two data immediately identifies which part of the pore space is the pore body and which is the throat, thus obtaining a model of the pore space. In the case of Berea sandstone, it is clear from Fig. 1 that the pore space consists of large cavities of about 85 μm and they are connected via 15-μm channels or throats. Such comparison is vital to understanding fluid flow in such porous media. Such experiments are the basis for a more detailed analysis of the pore shape (Chen and Song, 2002).

Experiments using DDIF have also been applied to carbonate rocks with complex pore geometry such as the presence of microporosity, which is intimately connected to the macropore network. The presence of micropores significantly affects the behavior of electrical conductivity and its interpretation; however, it is not easily identified by other methods, including NMR relaxation measurements. This is because the water molecules diffuse in and out of the micropore regions into the macropore regions with the lifetime of the spin relaxation. As a result, the relaxation time reflects an average surface/volume ratio instead of identifying the two pore systems with different relaxation times. The DDIF method is, on the other hand, directly based on diffusion dynamics and so does not suffer from such an averaging effect and it has been shown to readily identify the presence of micro-, meso- and macropores in many carbonate rocks (Song et al., 2000, 2002a; Kenyon et al., 2002). It has also been used to characterize the evolution of a pore network and the growth of microporosity due to diagenesis, the relationship between pore body and throat, and pore filling and drainage processes. This method has also been applied to biological porous media to study their structure and diffusion dynamics (Sigmund et al., 2008; Lisitza et al., 2007).

Two-phase flow (air and water) has been studied by a combination of MRI and DDIF techniques (Chen et al., 2003). Due to its limited resolution, MRI does not address water saturation at the pore scale. The DDIF data show drastically different pore fillings in different partial saturation processes, such as cocurrent and countercurrent imbibition. For example, in countercurrent imbibition, water was found to only fully fill the small pores while leaving significant air space in the large pores. Such partial saturation is very relevant both in oil reservoirs with the presence of oil and water and also in the vadose zone, where water and air occupy the pore space. In particular, water and air saturations change continuously from the top of the soil to the depth of the water table.

Two-Dimensional Nuclear Magnetic Resonance Spectroscopy

Spin relaxation measurements have been the key technique for characterizing porous media for many years. The principle of such measurement is based on the theory of surface relaxation (T1 and T2 are longitudinal and transverse relaxation times) and the measured relaxation rate reflects the surface/volume ratio of the sample (Cohen and Mendelson, 1982): 
\[\frac{1}{T_{1}}{=}\frac{1}{T_{1b}}{+}\mathrm{{\rho}}_{1}\frac{S}{V_{\mathrm{p}}}\]
[5]
and 
\[\frac{1}{T_{2}}{=}\frac{1}{T_{2b}}{+}\mathrm{{\rho}}_{2}\frac{S}{V_{\mathrm{p}}}\]
[6]
where T1b and T2b are the bulk relaxation times of the fluid, ρ1 and ρ2 are the surface relaxivity for the two relaxation processes, S is the rock surface area, and Vp is the pore volume. A rock consisting of multiple pore sizes (thus different surface/volume ratios, SVRs) would exhibit a distribution of relaxation spectrum corresponding to these different SVRs and the pore sizes can be characterized by T2 or T1 distributions. For example, it is well-known that the signal at T2 < 33 ms corresponds to the capillary bound water in small pores (Morriss et al., 1993; Kenyon et al., 1995; Kenyon, 1992).

This type of measurement works well in single-phase fluid-saturated rocks. When saturation is more complex, however, such as air and water in the vadose zone, oil and water in a petroleum reservoir, or water and organics in contaminated soils, the interpretation of the relaxation spectrum can be very difficult. For example, the T2 of oils and organic components can vary due to their molecular weight and structures so that they overlap with that of the water signal. In this case, the T2 distribution alone is no longer an effective means to differentiate them. Also, mud invasion further complicates the identification of different fluids at near-bore regions.

To further enhance the resolution of different fluids, two-dimensional NMR has been developed in recent years to combine diffusion and relaxation measurements to obtain the joint distribution of the diffusion constant and relaxation, i.e., DT2 map. It is conceptually different than the one-dimensional experiment in that a matrix of signals is acquired as a function of two independent variables. An example of the two variables could be that one is sensitive to relaxation and the other sensitive to diffusion. A specially designed two-dimensional Laplace inversion algorithm (Song et al., 2002b) was used to obtain the two-dimensional spectrum of diffusion and relaxation. Many variants of this two-dimensional NMR concept have been developed recently to study polymer mixtures, liquid crystals, crude oils, pore structures, cements, and other materials (Song et al., 2008a). Laboratory and field experiments are being performed routinely. An example of the two-dimensional DT2 map is shown in Fig. 2 , illustrating the ability to separate water and oil phases in a porous rock (Hürlimann et al., 2002). Many studies have reported application of this general method for a variety of subjects, for example, petroleum fluids (Hürlimann et al., 2009), food science (Hills, 2006; Hürlimann et al., 2006), cements (McDonald et al., 2005), chemical exchange (Lee et al., 1993), liquid crystals (Hubbard et al., 2005), and porous media (Washburn and Callaghan, 2006; Song et al., 2008b).

Imaging and Heteogeneity Analysis

Nuclear magnetic resonance relaxation and diffusion methods are important to characterize pore-scale properties of rocks and fluids. Magnetic resonance imaging, on the other hand, allows the imaging of the interior of the rocks to obtain the spatial distribution across a much larger scale. Figure 3 shows cross-sectional images of several rock samples obtained on a 2T MRI system (Pomerantz et al., 2008). The images demonstrate a large range of heterogeneity behavior such as uniform and heterogeneous at the millimeter scale and the centimeter scale.

Although visual inspection of the images provides an intuitive understanding of the heterogeneity (such as the magnitude and its length scales), however, it is highly desirable to discuss heterogeneity quantitatively to make quantitative comparisons with different samples. A new method has been recently reported to quantitatively interpret the statistics of heterogeneity in terms of the magnitude and the length scale of the heterogeneity composition. The basic idea is to use the two-point correlation function to describe the statistics of the heterogeneity: 
\[\mathrm{{\gamma}}(h){=}\frac{1}{2N(h)}{{\sum}_{i{=}1}^{N(h)}}\left[z(x_{i}){-}z(x_{i}{-}h)\right]^{2}\]
[7]
where h is the lag separating the two points in the pair, z(x) is the quantity of interest measured at position x, and the sum runs across all position pairs with distance h; N(h) is the number of pairs with separation h. This function γ is called the variogram in geostatistics. If the sample heterogeneity is of a single length scale, a model function such as the exponential function can be used to describe the decay of the correlation. For samples with multi-length-scale heterogeneity, we expect γ to be a sum of several length scales: 
\[\mathrm{{\gamma}}(h){=}{{\sum}_{i{=}1}}C_{i}\mathrm{{\gamma}}_{\mathrm{model}}(h,a_{i})\]
[8]
where γmodel(h,a) is the intrinsic variogram with a length scale of a, and the full variogram of the sample is the superposition of a series of variograms with different length scales. The coefficient Ci is the weighting factor for each length scale and is called the heterogeneity spectrum, and it is the goal of this heterogeneity characterization.

The above definition of the heterogeneity analysis assumes infinite resolution of γ(h) and unlimited sample size to characterize full heterogeneity scales. Experimental variograms will always suffer, however, from finite resolution and finite sample size. As a result, any experimental variogram will have reduced sensitivity for heterogeneity scales much shorter than the resolution and much larger than the sample size. To interpret an experimental variogram, such a sensitivity function was developed in recent studies to allow inversion of the measured variogram to obtain the heterogeneity spectrum, as shown in Fig. 3. The heterogeneity spectra quantitatively show the amount of heterogeneity in absolute units and as a function of the length scale. Further details of the theory and experimental analysis can be found in Pomerantz et al. (2008, 2009).

The need to visualize spatial characteristics is ubiquitous in studies of natural materials. The heterogeneity spectrum method can be applied to any imaging technique and most importantly allows direct comparisons of different techniques and different measurement scales (Pomerantz et al., 2009).

Inversion: Resolution and Uncertainty

Typically, NMR relaxation and diffusion are manifested as decaying signals. Data analysis often involves Laplace inversion to obtain a spectrum (or distribution) of relaxation times or diffusion constants. The inversion is ill conditioned in the sense that small noise in the data can cause large changes in the spectrum. Similar ill-conditioned problems are of great importance in many other fields of science and engineering, for example, geophysics and astronomy. One of the critical issues regarding such algorithms is the resolution of the resulting spectra. For example, given a spectrum with a single peak, one may ask the question whether the spectral width is determined by the finite signal/noise ratio of the data or does it reflect the true width of the underlying phenomenon? A second related question can be: how far away must two peaks be to be identified as two independent contributions?

The mathematical framework for such inversion can be described generally as the following equation: 
\[\mathbf{\mathrm{M}}{=}\mathbf{\mathrm{KF}}\]
[9]
where M is a vector of measurement data, F is the discretized relaxation spectrum, and K is the kernel matrix. The goal of the inversion is to obtain F from M. The kernel function is known for a specific experiment. For example, for T2 relaxation, Kij = exp(−τi/T2j), where τi is the decay time for the ith data point and T2j is the jth component of T2. Many of the properties of the inversion are determined by the structure of the kernel matrix, which is governed by the functional form of the kernel and the choice of the τ values.
One of the properties of the K matrix is its singular value decomposition (SVD). The SVD is a method to decompose a non-square matrix into intrinsic singular values and vectors in a fashion similar to the eigenvalues and eigenvectors of a square matrix. The result of such decomposition is that the spectrum F can be written as a sum of several singular vectors. Due to the ill-conditioned nature of Laplace inversion, however, the series of singular values used for the decomposition needs to be truncated at the smallest singular value determined by the data noise. The corresponding singular vector, in fact, provides the highest resolution in the T2 domain (Song et al., 2005). Using this procedure, we obtained an estimate of the spectral resolution using 
\[\mathrm{{\Delta}{\lambda}}{\sim}\frac{\mathrm{{\lambda}}_{\mathrm{max}}{-}\mathrm{{\lambda}}_{\mathrm{min}}}{1.2\mathrm{log}\left[M(0)/\mathrm{{\sigma}}\right]}\]
[10]
where λ = log10(T2), λmax and λmin correspond to the maximum and minimum of the T2 range, M(0) is the signal at zero time, and σ is the noise variance of the data. This result can be used to determine whether an observed spectral width is the true signature of the data or is due to limited signal/noise of the data so that it is not reflective of the underlying signal.
Another issue with Laplace inversion is the difficulty in determining the uncertainty of the inversion results due to the nonlinear nature of the inversion. Two general methods have been discussed in the recent literature based on the understanding of the noise characteristics of the data. The essence of the ill-conditioned problem is that many solutions fit the experimental data within the allowed noise of the data. Statistical analysis of this ensemble of solutions allows the determination of the uncertainty of the inversion results, such as quantities (I) that can be evaluated from the spectrum: 
\[I{=}{{\sum}_{j}}w(T_{2j})\mathrm{F}(T_{2j})\]
[11]
where w(T2j) is the an arbitrary function of T2. For example, to obtain porosity, w = 1; for capillary bound fluid, w = 1 for T2 < 33 ms and w = 0 for T2 > 33ms.

The first method uses a non-negative least squares (NNLS) algorithm to identify the solution ensemble and to obtain the maximum and minimum of such a quantity (I) within the ensemble (Parker and Song, 2005; Song, 2007). The solution ensemble is characterized by the least-square misfit (χ2) to be less than Aσ2, where the numerical factor A determines the confidence level of the ensemble. Using NNLS, a minimum misfit is calculated as a function of I [X2(I)]. In general, the curve of X2 vs. I is concave and the intersection with X2 = Aσ2 defines the boundary of acceptable solutions. These two intersection points thus determine the maximum and minimum of the allowed values of I. This technique considers all solutions constrained only by the data and much more than just the regularized solutions. It does not directly find all the solutions, however, but evaluates the boundary of such an ensemble and thus can be calculated fast.

The second technique directly evaluates many solutions to achieve a representation of the entire solution ensemble (Prange and Song, 2009). It uses a Monte Carlo algorithm to produce a random solution sample (e.g., a test T2 spectrum) and then uses the χ22 to evaluate the probability density function of such a solution sample. An efficient sampling algorithm was developed in order to speed up the Monte Carlo method and 10,000 samples can be obtained within a few seconds. Such a Monte Carlo method allows a direct visualization of the extent of the solution variations as well as quantification of any inversion results (such as I) and their uncertainty. Compared with the NNLS method, the Monte Carlo method determines the probability distribution instead of only the boundaries.

Even though these inversion techniques are discussed in term of Laplace inversion, their application is quite general for many other ill-conditioned problems such as imaging of shallow groundwater using the magnetic resonance sounding technique.

Conclusions

Magnetic resonance techniques have seen considerable development in recent years primarily due to demand in the oil and gas industry for formation evaluation and reservoir characterization.These advances can be straightforwardly transferred to the application of shallow groundwater storage and management. Nuclear magnetic resonance spectroscopy is sensitive to water saturation and its quantitative partition in different pools (clay-bound, capillary-bound, and movable water) and is thus uniquely suited for understanding water movement in aquifers and well productivity. Thus NMR should become one of the routine geophysical tools for water prospecting and monitoring.