Longitudinal and transverse nuclear magnetic resonance (NMR) relaxation signatures in porous rock were simulated on the microscale to examine and quantify how physical hydrologic parameters, such as rock-surface properties and pore sizes, affect longitudinal and transverse NMR signals of real, complex media. Parameters studied were: magnetic field strength, rock susceptibility, pore coupling, and surface reactivity. Using the finite element method (FEM), simulations of the spatial- and time-dependent magnetization evolution in arbitrary pore geometries, diffusion regimes, and heterogeneous distributions of rock surface properties, i.e., surface relaxivity, were compiled using an adapted generic diffusion model coupled with magnetic gradient field calculations. The numerical simulations were validated using analytical solutions that are available for simple pore geometries. We observed a pore-size-dependent ratio of transverse T2 and longitudinal T1 relaxation times, and thus a pore-size-related and rock-susceptibility-dependent effective transverse surface relaxivity was deduced. This can be used to improve estimates of pore sizes and thus of permeability from transverse NMR relaxometry measurements. Simulations of connected pore systems showed significant influences of interpore coupling at hydrologically relevant pore sizes, e.g., fine sands. Depending on the dominant diffusion regime, the typically heterogeneous distribution of surface relaxivities in rocks and sediments, i.e., geological noise, can lead to a significant underestimation of derived pore sizes and thus of permeability.

To advance the fundamental understanding necessary to translate measured nuclear magnetic resonance (NMR) relaxometry data into pore-space properties, NMR signals were simulated at the microscale using a finite element analysis. Subjects studied were a pore-size- and frequency-dependent ratio of transverse and longitudinal relaxation times, coupling, and heterogeneity of surface relaxivity.

In petrophysical applications of NMR, the measured relaxation signals originate from the fluid-filled pore space. Hence, in water-saturated rocks or sediments, the water content directly corresponds to the initial amplitude of the recorded NMR relaxation signal. The relaxation rate (longitudinal/transverse relaxation time ratio T1/T2) depends on the pore size and physicochemical properties of the rock–fluid interface (surface relaxivity) and on the concentration of paramagnetic ions or molecules (e.g., dissolved O2) in the fluid phases (bulk relaxivity) (e.g., Bryar et al., 2000; Keating and Knight, 2007). Therefore, NMR can provide estimates of hydraulic parameters, i.e., pore size distribution or permeability. In this respect, there is a large variety of empirical or semiempirical relations (Kenyon et al., 1988; Dunn, 2002; Sen et al., 1990). Other relations frequently used as permeability estimators from NMR are based on geometrical considerations such as the Kozeny–Carman model (Carman, 1956, Ch. 1). These permeability estimators apply the analytical relations between pore size and NMR relaxation behavior derived by Brownstein and Tarr (1979) for single unconnected pores. These relations, however, are limited to simple, idealized pore models such as cylinders, spheres, or planar shapes. Furthermore, the influence of internal gradients on the transverse relaxation processes is usually not accounted for. As Pape et al. (2009) demonstrated for low-porosity reservoir rocks, however, considering the effects of internal gradients can significantly improve the estimation of pore sizes, and thus of permeability, from T2 relaxation times.

To better understand the fundamental quantitative correlations between complex porous media and their respective NMR parameters such as relaxation times or surface relaxivity, it is necessary to revert to numerical simulations. This work aimed at improving the basic understanding of these parameters based on pore-scale simulations, thereby advancing the interpretation of NMR data by avoiding restrictive approximated interpretation schemes, e.g., for deriving pore size distributions, connectivity, inner surfaces, or permeability.

Nguyen and Mardon (1995) initially proposed using a FEM to simulate NMR processes at the pore scale and to eventually couple these equations with transport equations (e.g., Navier–Stokes) in a porous medium to investigate and quantify the links between NMR and permeability. Due to the limited computing capacities at that time, however, their investigations where restricted to small, simple, two-dimensional pore geometries and respective finite element meshes. In this study, we developed a multiphysics approach using a FEM to numerically simulate NMR T1 and T2 relaxation data at the microscale to study how physical and hydrologic parameters affect NMR relaxometry data. The studied parameters that influence NMR signals are: internal field gradients, pore connectivities, and various diffusion regimes. The COMSOL Multiphysics programming environment (COMSOL, Burlington, MA) was applied to simulate the arbitrary model geometries, thereby allowing a multiphysics three-dimensional modeling of the NMR behavior of the porous medium model along with the electromagnetic properties of the model. Nuclear magnetic resonance signatures can be simulated in both time and eigenvalue domains, which allows NMR relaxation times and amplitudes to be directly derived without the need for prior processing of the time-domain data, e.g., by inverse Laplace transformation. This approach thus reduces computational costs and avoids ambiguities when estimating parameters by inverse modeling.

Theory of NMR Relaxation in Porous Media

The phenomenon of NMR occurs when the nuclei of certain atoms immersed in a static magnetic field B0 are then exposed to a second oscillating magnetic field. In hydrologic and petrophysical practice, usually the nuclear magnetization of H1 in the pore fluid (e.g., water or oil) is relevant (Dunn, 2002; Kenyon, 1997). The relaxation of observable NMR signals in pore fluids results from the interaction of nuclear spins with their surrounding environment, i.e., the bulk fluid and pore walls.

Bloch–Torrey Equations

The evolution of magnetization during excitation and the subsequent NMR relaxation of the energized nuclei in free fluid is well described with the phenomenological Bloch–Torrey equations (Bloch, 1946; Torrey, 1956). They describe the x′, y′, and z components of the magnetization vector m as a function of time in the rotational frame of reference—precessing along the z axis with the angular Larmor frequency ωL = 2πγB0—according to 
where m = m(r,t) denotes the magnetization density vector in the rotational frame of reference and γ is the gyromagnetic ratio of the protons. The diffusion tensor D simplifies to D = ‖D‖ in the case of isotropic diffusion being assumed in the presented simulations. External magnetic fields in addition to B0 = (0,0,B0)T = constant are denoted by Bex = [Bx,By,G(r,t)r + Δω(r)/γ]T, where B1 = Bx + jBy is the complex oscillating radio frequency (rf) field of the NMR excitation pulse and Δω denotes the frequency offset between B1 and ω0. The global external gradient field vector of the static magnetic field is given by Gk(r,t) = [Bk(r,t) − B0k]/rk with k = x,y,z.
These relations represent the fundamental equations of motion of the magnetic moment vector m by considering diffusion plus an additional sink term governing the attenuation of the system to equilibrium. The relaxation constants T2 and T1 reflect the homogeneous decay of the macroscopic nuclear magnetization (spin-spin or transverse relaxation time constant, T2) and energy exchange with the surrounding environment (spin-lattice or longitudinal relaxation time constant, T1) returning the system to equilibrium magnetization density m0(t = 0) = m0uz. After termination of the rf excitation field, the external magnetic field vector simplifies to Bex = [0,0,G(r,t)r]T and Eq. [1] can be written as 
where the complex transverse magnetization mmx + jmy describes the free induction decay (FID) of the transverse magnetization in the observational xy plane that, in practice, is mainly dominated by signal dephasing of the total magnetization M due to field inhomogeneities (local gradients) in the static field B0. To eliminate the complex dephasing term of the transverse magnetization in Eq. [2a], and thus refocusing M, a second excitation pulse shifted by 180°, i.e., jγrGm, can be applied after a time t = τ stimulating a Hahn echo (Hahn, 1950). The repetition of such a procedure results in the well-known CPMG echo train (Carr and Purcell, 1954; Meiboom and Gill, 1958), which is commonly used in T2 relaxation experiments. The total transverse magnetization experiences further irreversible degradation due to dephasing, however, as the molecules diffuse through the inhomogeneous magnetic field during the time interval τ. To take into account such diffusive motion of the excited spins in the presence of a magnetic gradient field during the time interval τ, an additional time-dependent sink term (irreversible) has to be incorporated into Eq. [2a], yielding 
for transverse relaxation with a diffusion gradient and FID (e.g., Dunn, 2002).

Magnetization Evolution in a Restricted Pore Space

Considering NMR relaxation in a confined system, e.g., a saturated porous medium, the relaxation rate of excited nuclei is further accelerated due to surface interactions at rock–fluid interfaces. To explain the NMR relaxation behavior in an environment where diffusion is restricted, e.g., a single isolated pore, Brownstein and Tarr (1979) introduced Neumann boundary conditions to Eq. [1] at the fluid interfaces, e.g., pore walls.

At the (active) pore surface—i.e., within a thin volume layer of thickness d in the range of a few Angstroms near the pore surface (Fig. 1 )—the excited nuclei experience an accelerated surface relaxation rate 1/τs. The specific mechanism for such relaxation depends strongly on the magnetic dipolar interactions of the spins in the fluid with magnetic dipoles within this layer or at the surface of the fluid–rock interface (Torrey, 1956; Foley et al., 1996). Such interactions are very short ranged and are commonly assumed to be in the range of a few Angstroms. The strength of the NMR relaxation at the surface is then characterized by the surface relaxivity ρs defined as ρsd(1/τs − 1/Tb), where Tb is the relaxation time of the bulk fluid and τs is the surface relaxation time within the surface fluid layer of thickness d.

The dominant mechanism governing the surface relaxation processes in porous rocks can be attributed to the hyperfine interactions of the fluid nuclei with paramagnetic centers at the rock surfaces, e.g., most commonly Fe3+ and Mn2+. Thereby, the transversal and longitudinal surface relaxation rates (1/τs1,2) related to fluid nuclei aligned with paramagnetic ions at the rock surface is characterized by different multipolar relaxation mechanisms such that ρs2s1 ≥ 1 with median experimental values of 1.59 (e.g., Korringa et al., 1962; Kleinberg et al., 1993). In this study, however, we have focused mainly on transverse relaxation processes enhanced by diffusion in inhomogeneous magnetic fields that arise from susceptibility contrasts between the pore fluid and rock matrix (internal gradients) and as such are not related to surface relaxation. Therefore, without loss of generality, we use ρs1 = ρs2 as a common intrinsic rock property input parameter on the ith fraction of active NMR rock surface Si, e.g., ρSi = 10 μm/s: 
\begin{eqnarray*}&&\mathrm{{\rho}}_{S}m(\mathrm{\mathbf{r}},t){+}D{\hat{n}}_{i}\mathrm{{\nabla}}m(\mathrm{\mathbf{r}},t)\mathrm{{\vert}}_{S_{i}}{=}0{\ }(\mathrm{active})\\&&D{\hat{n}}_{i}\mathrm{{\nabla}}m(\mathrm{\mathbf{r}},t)\mathrm{{\vert}}_{S_{i}}{=}0{\ }(\mathrm{passive})\end{eqnarray*}
where i is the surface outward normal at Si. The surface relaxivity ρs is then a measure of the loss of magnetization and reflects the molecular interactions at the fluid boundary S. In this respect, the gas–fluid interface in a partially saturated system would have to be treated as a passive boundary. In hydrologic considerations, the corresponding fluid phase would be considered an active volume (water), whereas the gas phase (air) would be a noncontributing passive volume, i.e., m0(t = 0) = 0. Note, that in case of an NMR active second gas or fluid phase, magnetization (spin) density m0,2 ≠ 0 and diffusion coefficient D2 would have to be set accordingly for such phase.

Summing up the introduced NMR diffusion Eq. [2b] and [3] together with the boundary conditions as given in Eq. [4] allows a comprehensive description of the evolution of magnetization in a porous medium for time-domain modeling. The formulation is valid for T1 relaxation processes and can also be used to calculate the FID and transverse relaxation behavior in the presence of internal and external magnetic gradients.

Considering an experimentally observable scale, the spatially distributed magnetization densities m integrate to a total effective macroscopic magnetization: 
where V is the total volume of the investigated pore space and M0 is the effective magnetization in equilibrium.
In practice, e.g., to reduce computation time and to directly derive relaxation constants, it can be convenient to solve the above differential equations not in the time domain but rather in terms of their eigenvalues by separating variables according to 
Substituting m with Eq. [7], the orthogonal eigenfunctions En(r) and eigenvalues λn ≡ 1/Tn basically satisfy 
\[{-}\left\{\frac{1}{T_{\mathrm{b}}}{-}D\left[\mathrm{{\gamma}\mathbf{G}}(\mathrm{\mathbf{r}})t\right]^{2}\right\}E_{n}{+}\mathrm{{\nabla}}D\mathrm{{\nabla}}E_{n}{\ }(T_{2})\]
\begin{eqnarray*}&&{-}\left\{j\mathrm{{\gamma}\mathbf{r}G}(\mathrm{\mathbf{r}}){+}D\left[\mathrm{{\gamma}\mathbf{G}}(\mathrm{\mathbf{r}})t\right]^{2}{+}\frac{1}{T_{\mathrm{b}}}\right\}E_{n}{+}\mathrm{{\nabla}}D\mathrm{{\nabla}}E_{n}{=}\\&&{-}\mathrm{{\lambda}}_{n}E_{n}{\ }(\mathrm{FID})\end{eqnarray*}
for the longitudinal and transverse NMR relaxation behavior. The boundary conditions for T1, T2, and FID are given by 
for active and passive boundaries, respectively. For simple geometries, e.g., spherical (Fig. 1), cylindrical, or planar shapes, Eq. [8a] can be solved analytically. From such considerations, the basic relations between NMR T1 relaxation times and the geometrical properties of isolated pores, e.g., the pore radius, were derived by Brownstein and Tarr (1979) and are commonly used for the interpretation of NMR relaxometry measurements. They described the relaxation behavior of magnetization in surface- and diffusion-limited regimes. In the case of a surface-limited (fast) diffusion regime and assuming ds/(vds) ≪ 1, a pore is characterized by a single eigenvalue λ, i.e., the NMR-labeled molecules diffuse much faster to the pore wall than the time it takes for the excited spins to decay to the equilibrium state; hence, 
for longitudinal relaxation and 
when considering transverse relaxation times T2app measured in CPMG, where DG is the irreversible diffusion gradient relaxation rate due to magnetic gradients, and s/v represents the surface/volume ratio of the pore. Furthermore, by assuming a constant gradient field G(r,t) = constant and linear time steps t = nΔt with a constant echo interval (echo time) Δt = 2τ, Eq. [8b] can analogously be solved analytically. Thus, the transverse relaxation behavior in the presence of a linear (external) magnetic gradient field can be derived. Considering nonlinear (internal) magnetic gradients, Eq. [8b] and [8c] can only be solved numerically.

The analytical solutions of the above equations, however, are only valid under the following conditions: (i) relaxation of the nuclei begins and ends within the same pore, i.e., pores act independently; (ii) surface relaxivity is constant; and (iii) pore geometry is known and constant throughout the investigated sample. Therefore, to assess the limitations of such approximated solutions and to improve the quantitative understanding of T1 and T2 NMR relaxation signatures of complex porous systems (e.g., connected pores, variations of surface properties due to different rock minerals, the impact of internal gradients on derived T2–based pore size distributions and permeability, partial saturation, etc.) the above NMR diffusion equations have to be solved numerically.

Calculation of the Primary Magnetic Fields (Internal Gradients)

It is well known that internal magnetic field gradients in saturated rocks originate from susceptibility contrasts between the pore fluid and the minerals in the rock matrix (e.g., Kleinberg et al., 1994; Coates et al., 1999; Godefroy et al., 2001). To model the increased decay rate observed in the transverse NMR relaxometry, i.e., CPMG echo trains and FID, due to the diffusive motion of the NMR-labeled molecules in a nonlinear magnetic gradient field (internal gradient), the spatial distribution of field gradients within the simulated fluid–matrix system can be calculated according to Ampere's law (Nabighian, 1987): 
and outer boundary conditions 
at the model limits and with inner continuity requirements 
at interfaces S1,2 between inner model domains such as the pore–matrix boundary. The global Larmor frequency f0 = γB0 = γμ0μrH0 is defined by the static external primary magnetic field H0. The gradient in the direction of the equilibrium magnetization B0 = Bz at the point r is given by G(r) = ∂/∂zB(r). The terms A, H, and B denote the local magnetic potentials, fields, and flux densities, respectively; μr(r) = 1 + κ(r) is the spatial distribution of the magnetic permeability (susceptibility κ) in the system. In this study, we set μrfluid = 1 and assumed mean values for μrrock throughout the rock matrix and thus μrrock > μrfluid/gas. Dissolved ferro- or paramagnetic ions such as Fe3+ or Mn2+, however, can significantly increase the magnetic permeability of the pore fluid. Actually, μr can also exhibit significant variations within a rock sample, e.g., due to individual grains or rock fractions with a high Fe content. Typically, susceptibility values for sediments and rocks range from κ ≤ 10−4 SI (e.g., clean sandstone or silica) up to κ ≥ 10−2 SI (e.g., volcanic rocks).

Numerical Realization

As the description of the evolution of magnetization in the pore fluid is very similar to classical differential equations for diffusion processes, approaches frequently used in fluid dynamic sciences can be used to solve Eq. [1–3] with boundary conditions according to Eq. [4] in the time domain, e.g., random walk strategies, lattice Boltzmann methods (LBMs), or FEMs.

We used COMSOL Multiphysics, a commercial program based on finite elements that includes simulation toolboxes for general convection and diffusion processes in arbitrary three-dimensional geometries. Optionally, these models can be linked to other physical simulations, such as corresponding magnetic field calculations (internal gradients), fluid flow, or heat transfer. In the FEM, the geometric model is divided into mesh elements (Fig. 2 ) and the above differential equations are solved for magnetization densities mi(t) at each mesh element i of the model (Zienkiewicz, 2005).

In contrast to random walk methods or LBMs, the FEM approach permits the NMR diffusion equations to be solved in their eigenvalue domain, thus significantly reducing the computational cost. Consequently, arbitrary pore geometries can be simulated and a direct derivation of NMR relaxation times (eigenvalues) can be derived directly without the need for consecutive inverse modeling in the time domain, e.g., by inverse Laplace transformation.

In this study, we considered only the evolution of magnetization densities of the NMR-labeled H2O molecules in the saturated pores, i.e., initially m(x,y,z,t = 0) = M0/V. The simulations were compiled using a medium performance dual core personal computer with 2Gb of RAM. Depending on the model design, the computational time was generally in the range of several seconds up to a few minutes.

Results and Discussion

Model Validation Using Analytical Solutions

To validate the simulations, we compared the numerical compilations with analytical solutions that are available for isolated pores with simple geometries. Simulation results for an isolated spherical pore in surface-dominated (one dominant decay time per eigenvalue) and diffusion-dominated (multiple decay times per eigenvalue) regimes using 4166 mesh elements are in a very good agreement with the analytical results with deviations of <0.1% (Fig. 3 ). In comparison, Nguyen and Mardon (1995) used an axial symmetry with eight nodes in their initial models.

To validate our simulations with respect to transverse NMR relaxation in a gradient field as stated in Eq. [3] and [8b], we simulated an unrestricted relaxation environment in free water with an intrinsic bulk relaxation time Tb = 3 s. The simulations used constant gradients G(r,t) = G = constant throughout the model, with a zero surface sink term at the model boundaries, i.e., ρs = 0. In Fig. 4 , the square of the magnetic field gradient is plotted vs. effective transverse relaxation times at different CPMG echo spacings for analytical (lines) and numerical (symbols) calculations. The results are in good agreement. As DG relaxation becomes the dominant term in Eq. [3] for increasing echo times 2τ and field gradients G, the slope of each curve becomes proportional to the diffusion constant of the fluid.

Frequency Dependence of Transverse Relaxation Time Due to Internal Gradients

In NMR applications, the primary static magnetic field B0—usually generated by strong permanent or electro-magnets—defines the resonance frequency (Larmor frequency fl) of the investigated nuclei, e.g., 1H isotopes. Typically, Larmor frequencies used in petrophysical NMR relaxometry range from kilohertz (e.g., the Earth's magnetic field: B0 ∼ 30−50 μT) up to a few 10s of megahertz (e.g., MRI: B0 ≥ 0.25 T). The resulting internal magnetic field gradients in rocks stem from susceptibility contrasts between the pore fluids and rock matrix or are due to minerals with different susceptibilities. Such internal gradients naturally increase as the field strength of the magnetic field increases, and thus an enhanced observed CPMG relaxation rate 1/T2 will be observed at higher Larmor frequencies. The dependence of the transverse NMR relaxation time on the Larmor frequency as a function of susceptibility contrast between pore water and rock matrix in a saturated porous rock were simulated for the model of a single capillary pore (radius R = 1 μm, ρs = 10 μm/s) for a frequency range of fl = 103 to 107 Hz (Fig. 5 ).

The simulations were compiled for a susceptibility range of κ = 10−4 SI (e.g., clean sandstone) up to κ = 0.1 SI (e.g., volcanic rocks). The simulated NMR relaxation times show the strong increase in relaxation rates with increasing field strength of B0 and thus reflect the impact of internal gradients on the resulting CPMG relaxation measurements in petrophysical NMR. Note that at very low B0 fields, e.g., with Larmor frequencies in the kilohertz range or less, the influence of the gradient fields becomes less dominant and eventually T2 draws near to the respective T1 relaxation times, even for the high κ contrasts utilized in this model. (In low magnetic fields and negligible internal gradients, T2T1). Therefore, NMR relaxometry experiments at lower magnetic fields, i.e., low-field NMR, yield a more accurate estimate of petrophysical properties from T2 relaxation times, with a tradeoff, however, toward lower signal amplitudes as the spin magnetization scales with magnetic field intensity. On the basis of such numerical case studies, the dependency of T2 on both rock susceptibility and Larmor frequency can be used to quantify and compensate for such interfering influences of high susceptibility contrasts when deriving pore size distributions from NMR relaxometry.

Enhanced Transverse Relaxation Time as a Function of Pore Size in the Presence of Internal Gradients

For Larmor frequencies fl ≤ 4 MHz—which is generally considered low-field NMR—the diffusion effects for T2 can be minimized by using a sufficiently small echo spacing τ in a CPMG experiment. Then T2 can be used as a T1 proxy to derive pore size distributions. In this respect, a constant T2/T1 ratio for all pore sizes is assumed (e.g., Kleinberg et al., 1993).

In the pore space, however, internal magnetic field gradients, and thus the DG, are at a maximum near the pore–rock interfaces, i.e., where a high contrast in κ occurs. In Fig. 6 , the spatial distribution of DG within a saturated capillary pore for a matrix susceptibility κmatrix = 10−2 SI, echo spacing τ = 0.5 ms, and Larmor frequency fl = 2 MHz is plotted. The relaxation times generally decrease toward the pore walls (dark blue), which significantly affects the resulting effective transverse relaxation time of this pore, i.e., T2 = 0.5 s (no gradient) and T2 = 5.5 × 10−2 s (with internal gradients). Following the above considerations in a small pore, a larger fraction of its total NMR-labeled volume is affected by the gradient field, and thus one can assume a higher impact of the diffusion gradient on T2 NMR signatures that originate from the smaller pores.

To study such a dependence, we simulated longitudinal (T1) and transverse (T2) NMR relaxation signatures of saturated capillary pores in the presence of internal gradients. The simulations presented in Fig. 7 were calculated for capillary pore geometries using a global Larmor frequency of 2 MHz and (frequently used) CPMG echo times τ = 0.5 ms. The T2/T1 ratios generally show a strong dependence on the pore radius that becomes more distinct at low surface relaxivities ρs at the walls. The T2/T1 ratios are more strongly affected at low ρs because of the long relaxation times for low ρs. In contrast, the short relaxation times at high ρs cannot be reduced much further by the diffusion gradient. Actually, as surface relaxivity usually reflects paramagnetic substances at the pore wall, e.g., hematite, an additional increase in ρs would be expected with increasing internal gradients in a rock sample due to an increased surface relaxation rate 1/τs. Even so, the results plotted in Fig. 7 support a strong relation between pore size and the T2/T1 ratio.

As a consequence of an increase in the rock's susceptibility, the resulting effective transverse relaxation times decline. As depicted in Fig. 8 , such a decrease in T2 strongly depends on the respective pore size. The simulations show a significant impact on the resulting relaxation rates at κ ∼ 10−2 SI for the large pores (R = 100 μm). A significant decrease in effective transverse relaxation for the small pore sizes (R = 1 μm) occurs already at κ ∼ 10−3 SI. In practice, it is commonly assumed that ρs does not vary significantly throughout the medium, i.e., ρs = constant for all signal-contributing pores, although there is no direct experimental justification. From Eq. [9a] and [9b] it follows that the surface relaxivity for a single pore is given by 
\[\mathrm{{\rho}}_{s_{k}}{=}\frac{T_{\mathrm{b}}{-}T_{k}}{T_{\mathrm{b}}T_{k}}\frac{v}{s}{\ }\mathrm{with}{\,}k{=}1,2\]
where v/s is the (active) volume to surface ratio of the pore. For a pore space consisting of several pore sizes, i.e., integrating over the entire pore space, Eq. [13] extends to 
where k is a weighted relaxation time, Ii and Ti are the relative signal contribution and respective relaxation time of the ith pore of the system (note that this relation is strictly only valid in a surface-limited diffusion regime), and V = Σvi and S = Σsi are the total volume and surface of the system, respectively.

For the simulations shown in Fig. 8 to 10910 , a simple capillary pore with radius R was used, thus v/s = R/2. Following from Eq. [14] with a known volume/surface ratio of the pore model, we can derive a pore-size-specific enhanced transverse surface relaxivity parameter ρ2s that increases with increasing matrix susceptibility (Fig. 9). The parameter ρ2s(κ, R, fl) is then also a function of the matrix susceptibility and Larmor frequency as well as the pore size and reflects the enhanced transverse relaxation rate of the pores due to the internal magnetic field gradients in the porous medium.

As illustrated in Fig. 10, on using a constant surface relaxivity value (diamonds), i.e., ρ2s = 10 μm/s, to estimate the respective pore sizes, the pore radii will be significantly underestimated when the susceptibility of the matrix increases beyond a pore-size specific threshold. This susceptibility threshold decreases with decreasing pore radii, and thus the size of small pores will be more strongly underestimated than large pores. In contrast, when using an enhanced surface relaxivity model that takes into account the rock susceptibility, the pore radii are estimated correctly (dots). The slight increase at high susceptibilities is due to the beginning of the transition from a surface-limited to a diffusion-limited regime of the relaxation processes, i.e., ρ2sR/D > 1. Generally, this susceptibility threshold will shift to smaller values for NMR experiments at higher Larmor frequencies and vice versa. These considerations concur with observations from laboratory studies on low-porosity reservoir rocks (Pape et al., 2009). In practice, to estimate the pore-size-dependent effective surface relaxivity by inverse modeling, the susceptibility and porosity (e.g., NMR porosity) of the rock sample need to be independently measured. In this context, developing an enhanced inverse modeling based on the above simulations is still an objective of ongoing research.

Pore Coupling

Common quantitative interpretation of NMR relaxometry data in porous media uses the classical relation of NMR relaxation times and pore sizes for single isolated pores as derived by Brownstein and Tarr (1979). Because the pores are usually connected in porous rocks, the concept of separated individual pores is limited for the investigation of coupling between pores, i.e., the diffusive motion of NMR-labeled molecules between pores. Due to the diffusion constant of water, i.e., the diffusion length, however, pore coupling effects become more significant in NMR relaxometry in rocks with small pore sizes, e.g., in low-porosity reservoir rocks, shales, or clay.

A simple model of a capillary pore (R = 10 μm) that is connected to four smaller pores (R = 5 μm) by thin channels (pore throats) with widths 10−2 μm is shown in Fig. 11A . The total volume of the four small pores is equal to the volume of the large pore; thus we obtain equal amplitude fractions for each pore size. Simulations of longitudinal relaxation (T1) in terms of eigenvalues were performed for different throat lengths and widths of this coupled system. Figure 11B shows the respective T1 relaxation times and normalized amplitudes for throat lengths ranging from 0.1 to 80 μm. The arrows indicate the corresponding direction of change with decreasing throat length.

The coupled system is characterized by two relaxation times that each decrease as coupling increases (red arrows). For high throat lengths (weak coupling), no significant diffusive exchange between pores takes place and the pores can be considered to be isolated, with relaxation times of 0.25 and 0.5 s. The same relative amplitudes of 0.5 for each pore size reflect the equal pore volumes of large and small pores. When decreasing the throat length and thus increasing the coupling strength, however, the fully coupled system will eventually be characterized by a single effective relaxation constant and a respective normalized signal amplitude of 1. Because molecular exchange by diffusion takes place between the pores, the pore sizes estimated by the classical approach fail to describe this coupled system. Note that as the dominant relaxation time of the system converges to an effective value, the lower part of the spectrum decreases along with its loss in signal strength. These implications, relevant for the NMR relaxation behavior of coupled pore systems at the pore scale, agree with the results for porous medium models obtained by McCall et al. (1991).

Heterogeneity of Surface Relaxivities

Porous rocks and sediments usually are composed of fractions of different mineral grains such as quartz, clay, or hematite. Naturally such minerals can have different properties with respect to their paramagnetic ion content at the inner surfaces and thus show highly different mineral-specific surface relaxivities, i.e., ρsquartz ≪ ρshematite.

In practice, to estimate pore sizes from NMR relaxometry relaxation data, effective mean values of ρs are derived with uncertainties due to measurement errors. Additionally, there will also be deviations caused by the mineral grain distributions within the rock matrix. Such natural noise usually cannot be accounted for, and information with respect to its influence and significance on the derivation of pore parameters is not well known.

To assess such influences for fast- and slow-diffusion NMR relaxation regimes, we set up a simple model box (Fig. 12 ) of 144 spherical mineral grains (grain radii R = 1 and 10 μm for fast and slow diffusion, respectively), which were assigned low (ρsa = 1 μm/s) and high (fast diffusion models: ρsb = 50 μm/s; slow diffusion models: ρsb = 100 μm/s) surface relaxivities, representing clean (e.g., quartz) and impure (e.g., hematite) mineral grains. These two surface properties were assigned randomly to the grains in the models. Figures 12A to 12D illustrate examples of random realizations (fast-diffusion model) of transverse magnetization density at the first CPMG echo (τ = 0.5 ms) for different fractions (0, 25, 75, and 100%) of high-relaxivity grains. Simulated data (fast-diffusion regime) of 50 random realizations for different fractions (0–100%) of high-relaxivity grains in the model are plotted in Fig. 13A (time domain) and 13B (eigenvalue domain). Dashed lines represent the limits for 0 and 100% high-relaxivity grain contents. Figures 14A to 14C illustrate three excerpts for high-relaxivity grain contents of 6, 20, and 60%. Here, a bimodal eigenvalue distribution at low ρsb grain contents, i.e., isolated occurrences of ρsb grains, converges to a weighted mean monomodal relaxation behavior with increasing high-relaxivity grain contents. Similar simulations compiled for diffusion-limited relaxation, i.e., grain sizes are upscaled by a factor of 10 and ρsb = 100 μm/s, are plotted in Fig. 15 . Here, at 0% high-relaxivity grain content, the relaxation is purely surface limited. Upon introducing high-relaxivity grain fractions to the model, the relaxation behavior becomes increasingly diffusion limited, and thus the system becomes characterized by an ensemble of relaxation times (e.g., Fig. 16 , bottom). Upon increasing the fraction of high-relaxivity grains, the characteristic eigenvalues λ of the system increase along with an increased scattering of mean relaxation times and subsequently derived surface relaxivities for the 50 model realizations (indices 1–50). As expected, the scattering appears generally stronger in the surface-limited diffusion environment. If the model's total volume and surface area are known, Eq. [14] can be used to estimate the surface relaxivities and thus derive pore sizes from the simulated relaxation data (Fig. 17 and 18 ). The effective surface relaxivity of the system increases exponentially with increasing content of ρsb grain surfaces.

The derived effective pore sizes for the fast-diffusion model, i.e., volume/surface ratios, were generally found to be in very good agreement of <1% with the model. In the case of intermediate to slow diffusion, however, the derived mean surface relaxivities were typically underestimated due to the resulting faster effective relaxation rates. Consequently, the deduced pore sizes were underestimated, exhibiting a maximum deviation (∼15%) from the model at 20% high-susceptibility grain fraction. This maximum deviation can be attributed to a stronger impact from clustering (e.g., compare Fig. 12B and 12C), which then promotes the development of two distinct slow and fast relaxation regimes within the model. Furthermore, coupling between pore fractions of the model will also contribute, however, being of only minor significance in the presented models.

Conclusion and Outlook

Pore-scale simulations of longitudinal and transverse NMR relaxation signatures, taking into account the enhanced T2 relaxation rates that arise in the presence of internal magnetic gradients in porous rocks and sediments, were implemented in a multiphysics finite elements program. The simulations can be compiled for arbitrary pore geometries with spatial distributions of surface and magnetic properties of the fluid and rock, and thus allow the fundamental NMR relaxation behavior to be studied in complex pore systems, e.g., derived from computed tomography. The finite element approach permits the NMR diffusion equations to be solved in their eigenvalue domain and thus the NMR relaxation times (eigenvalues) to be directly derived. The computational cost can thereby be significantly reduced and simulations can be compiled in two and also in three dimensions using a standard personal computer.

Simulations of transverse NMR relaxometry data support previous observations (e.g., Pape et al., 2009) indicating that the use of T2 distributions as a permeability estimator tends to overestimate the contribution of smaller pores and thus to underestimate the rock permeability. To compensate for this drawback, the concept of an enhanced transverse surface relaxivity parameter that depends on pore size and rock susceptibility was introduced.

Modeling of the transverse relaxation behavior reflects the frequency dependence of internal gradients in the rocks, which, in principle, could be used to derive effective surface relaxivities by inverse modeling or enhancing existing (semi)empirical models. Therefore, it seems promising to apply such an approach to data obtained from frequency-dependent NMR relaxometry data, e.g., utilizing transverse NMR data recorded at two different Larmor frequencies, to better estimate pore-size distributions and permeability. In this context, from the observed shifts of the transverse relaxation times at different frequencies, it seems promising to quantify frequency- and pore-size-related effective surface relaxivities and thus derive more reliable permeability estimates, in particular, for low-porosity reservoir rocks, shales, or clay. Direct simulations of connected pores at the pore scale confirm an averaging of decay times along with an overall shift to higher decay rates related to the coupling strengths within the pore system. Simulated NMR signatures of random grain mixtures indicated that, depending on the dominant diffusion regime, the commonly heterogeneous distribution of surface relaxivities in rocks and sediments, i.e., geological noise, can lead to a significant underestimation of derived pore sizes and thus of permeability.

On the basis of such pore-scale simulations, we aim to implement an NMR interpretation scheme for T2 relaxometry data that will take into account pore-size-dependent transverse relaxation rates due to susceptibility contrasts between the rock and pore fluid. Using the presented finite element approach, simulations can be readily compiled for multiphase systems and thus will also be used in further studies to correlate the NMR properties and hydrologic parameters of partially saturated rocks, e.g., to derive the relative permeability and to study the transport properties of film flow. In a next step to further investigate and quantify the relationship between NMR and the transport properties of porous media, we will jointly simulate NMR and fluid flow (e.g., Navier–Stokes).

This study was conducted within the framework of the Transregional Collaborative Research Centre 32 (SFB TR 32). We thank the German Research Foundation (DFG) for funding.

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