Numerical analysis of nuclear magnetic resonance relaxation–diffusion responses of sedimentary rock

The nuclear magnetic resonance (NMR) relaxation–diffusion response of porous reservoir rock is frequently used, e.g. in oil field applications, to extract characteristic length scales of pore space or information about saturating fluids. External gradients are typically applied to encode for diffusion. In reservoir rocks, field inhomogeneities due to internal gradients can even at low fields be strong enough to interfere with this encoding. Furthermore, the encoding for diffusion coefficients of fluids takes a finite amount of time, during which diffusing fluid molecules can experience restricted diffusion. Both effects can combine to make the interpretation of the diffusion dimension of a relaxation–diffusion measurement difficult. We use x-ray-CT images of porous rock samples to define the solid and fluid phases of reservoir rock and simulate the full experimental pulse sequence, taking into account the static applied field, external gradients and internal gradients as a function of susceptibility of each component, and surface and bulk relaxation properties of fluids and fluid–fluid and fluid–solid interfaces. We carry out simulations of NMR relaxation–diffusion measurements, while explicitly tracking the time-dependent diffusion coefficient in each fluid as well as associated local gradients. This allows us to quantify the influence of restricted diffusion and internal gradients for common choices of experimental parameters.


Introduction
Nuclear magnetic resonance (NMR) methods offer a non-destructive way of inferring geometric characteristics and fluid properties of the spatial distributions of proton-containing fluids in the pore space of porous materials. A particular experiment, the measurement of diffusion-relaxation (T 2 -D) correlations, has been introduced relatively recently to the petroleum industry [1,2]. It uses the fact that the diffusion coefficients of some fluids, such as oil and water, are typically quite different to separate the relaxation responses of the different fluids using the additional diffusion dimension. This approach combines the application of T 2 relaxation analysis [3]- [8] with an encoding for the time-dependent diffusion coefficient of saturating fluids [2], [9]- [12], thus combining the advantages of T 2 relaxation measurements with diffusion measurements [13]- [15].
The T 2 -D experiment consists of two distinct pulse sequences in succession, a diffusion encoding sequence followed by a Carr-Purcell-Meiboom-Gill (CPMG) sequence used for the detection of the signal while also recording the T 2 relaxation decay (see figure 1). A short time spacing in the CPMG train recording T 2 relaxation is chosen to suppress diffusion effects by frequent refocusing of the transverse relaxation. In addition, the diffusion encoding sequence only includes short time intervals during which T 2 relaxation is present, resulting in a separable kernel describing the physics of the problem. We consider, in particular, the pulsed field gradient-stimulated echo (PFGSTE) encoding sequence [16], which has advantages in keeping the kernel separable even for longer diffusion intervals by lessening the effect of relaxation during the diffusion encoding interval. It is useful for systems in which the longitudinal relaxation rate T 1 is larger than the transversal relaxation rate T 2 . The diffusion kernel k 1 of the experiment for a homogeneous background gradient g b and a homogeneous applied gradient g a , assuming free diffusion, is given by [10,15,16]  Tanner NMR-pulsed field gradient-stimulated spin-echo sequence for an applied variable field gradient g a over a constant time interval δ in the presence of constant background gradient(s) g b , see (1)- (4).
where |g a | = g a , |g b | = g b , t d = + 2τ is the total diffusion encoding interval time from the first π 2 pulse to the first echo, γ is the magnetogyric ratio, τ is the spacing between the first two π 2 pulses, δ is the length of the gradient pulse, and δ 1 and δ 2 are the short times preceding and following the gradient pulse with τ = δ 1 + δ + δ 2 [10]. is the storage time between the second and third π 2 pulse, where only T 1 relaxation is present. The detection kernel k 2 of the sequence is a standard CPMG train described by with t E being the CPMG echo spacing and k noting the kth echo, resulting in the combined kernel K = k 1 k 2 . This formulation allows one to invert for diffusion coefficient and relaxation times using a two-dimensional (2D) inverse Laplace transform [17,18]. Here, we use the L-curve method combined with a non-negative least squares (NNLS) algorithm solving an amplituderegularized inversion problem [19,20]. It is well known that in porous rocks both restricted diffusion and internal gradients play an important role. Restricted diffusion leads to a reduction in attenuation in an applied gradient field by preventing the diffusing spins from testing the full variability of the inhomogeneous magnetic field [21]- [25]. Susceptibility differences between rock and fluid components generate internal magnetic fields increasing field inhomogeneity [8,15], [26]- [29], contributing to the loss of signal attributed to diffusion. These effects lead to the measurement of an apparent rather than absolute diffusion coefficient spectrum [15,30]. The effects of internal gradients might be accounted for by explicitly encoding for them [20,31] or by experimentally suppressing the internal gradients (e.g. [10]), which is not always possible or practical. We see from equation (1) that the balance between internal and external gradients can be changed either by adjusting the timing of the individual terms by influencing the coefficients C a , C b or C ab or by varying the strength of the applied gradients or background gradients, the latter e.g. by varying the static field strength. It should be noted that for internal gradients, g b cannot be assumed to be homogeneous and is a local quantity. This, in combination with restricted diffusion, results in various subtleties in the interpretation of T 2 -D relaxation-diffusion measurements. Significant ones include the fact that the diffusion in 4 the internal field might involve a shorter time interval than the diffusion encoding dimension, which includes a T 1 relaxation interval and thus acts in general on a different scale if working with a device providing a permanent gradient. Likewise, the presence of restricted diffusion can easily invalidate the assumption of the Gaussian phase approximation often used in the interpretation of T 2 -D maps.
Very high internal gradients with distribution tails exceeding 1000 G cm −1 have been reported for reservoir rocks at low fields [26]. It is clear that an accurate analytical treatment of the T 2 -D NMR response in porous rocks is not easy [30,32]. To aid quantitative treatment, it is advantageous to model the NMR response numerically. This allows one to control physical and experimental parameters explicitly, while attempting to match experiment as closely as possible. NMR random walk simulations go as far back as Carr and Purcell in 1954 [33]. Newer treatments include the effects of relaxation [34]- [37], diffusion [38] and internal gradients [39,40] in reservoir rocks. In earlier work, the surface relaxivity typically was treated as a free parameter to match experiments and numerical simulation, leading to much larger surface relaxivity in numerical simulations compared to experiments. Recently, it has been shown that a consistent set of parameters with a constant surface relaxivity can be used to match experiments and numerical simulations at different field strengths by explicitly modelling the internal magnetic field effect in complex geometries [40].
In this work, we extend the approach of Arns et al [40] for the T 2 relaxation response to T 2 -D responses. We use high-resolution x-ray-CT images [37,41] to define the distribution of the solid and fluid phases of the reservoir rock, and simulate the full experimental pulse sequence, taking into account the static applied field, external gradients, internal gradients as a function of the susceptibility of each component, and surface and bulk relaxation properties of fluids and fluid-fluid and fluid-solid interfaces. We directly compute the restricted diffusion coefficient spectrum and the apparent diffusion coefficient spectrum, and quantify the change to the diffusion spectrum under the influence of internal gradients for Berea and Bentheimer sandstone.
This paper is organized as follows. In the next section, we introduce the samples used in this study, followed by a description of the modelling of sub-micrometer structure. We then summarize the calculation of physical properties. Finally, we present the results followed by a discussion and conclusion section.

Image acquisition and processing
In this study, we consider two sandstone samples, Bentheimer sandstone and Berea sandstone, which are often used as benchmark rocks in the petroleum industry. Both samples were 5 mm cylindrical plugs and were imaged on the ANU x-ray-CT facility at about 3 µm resolution. In addition, x-ray powder diffraction (XRD) analysis was performed to determine the mineral content of the samples. Slices through the x-ray-density fields of the samples are depicted in the left column of figure 2. To simulate different fluid saturations inside the pore space, we use the covering radius map [42] as an approximation for the long-time equilibrium distribution, which might be experimentally achieved by a slow imbibition process and/or fluid redistribution by diffusion (middle and left columns of figure 2). Alternatively, direct imaging of fluid distributions might be used [43]- [45]. , clay region (dark grey) and pore space (black and white). The pore space is partitioned into two fluids using a morphological approximation to fluid distributions, white being the non-wetting fluid. Wetting fluid saturations are ≈50% (middle) and ≈25% (right).

Internal magnetic field calculation
The internal magnetic field is calculated in a dipole approximation [46] by assigning each voxel in the image an effective isotropic magnetic susceptibility and convoluting the dipole field with the susceptibility field. The dipole field B dip for r a (a is the radius of the dipole and r is the distance from the dipole center) is given as for unit volume. Here, µ 0 = 4π × 10 −7 N A −2 is the magnetic permeability of vacuum and m is the magnetic dipole moment for a unit volume V = 3 , with being the lattice spacing, Inside the sphere (r a), we have If a porous sample is inserted in a static applied field of strength H , the susceptibility field χ (r) results in an induced internal magnetic field B int . For sufficiently small susceptibilities, this field can be calculated using the convolution operation Given the material distributions from x-ray-CT analysis and fluid saturation simulations (see figure 3(b)), and the constituents of the samples from XRD analysis, we can calculate the internal fields at a desired field strength (9). The relevant magnetic susceptibilities are summarized in table 1. We calculate the effective susceptibilities of the clay regions assuming a brine fraction of 50% and using a volume-weighted average for χ clay , The resulting intrinsic internal gradient distributions at image resolution can be calculated using assuming that the applied static magnetic field B 0 is oriented in the z-direction and |B 0 | |B int |.   (9) for a sphere discretized on a Cartesian regular grid with spacing = 1 15 a for a given susceptibility contrast and a static magnetic field in excellent agreement with theory (6). (b) Sketch of the material distributions and material property assignments inside the porous partially saturated rock. Here, χ i stands for the susceptibility of the ith mineral, which in general can be a tensor.  The maximal internal gradient calculated for Berea sandstone is about 90 G cm −1 , consistent with experimentally derived maximum gradients for Berea 400/500 of about 75 G cm −1 [14]. We emphasize that we used an average susceptibility for the clay regions. Larger-scale susceptibility contrast and heterogeneity would result if the clay minerals were distributed differently with potentially higher maximum gradients. No diffusion-weighted average is taken here. Both Bentheimer and Berea sandstone have some clay fraction, which is modelled as a brine-saturated effective medium with 50% porosity. For Berea, this effective medium is paramagnetic, while for Bentheimer, the clay regions are diamagnetic, leading to much stronger susceptibility contrast in the Berea sandstone. In addition, we note that the gradient distributions for pure fluids are quite similar in shape due to the similar susceptibilities. Distinct partitions into low-and high-gradient regions result from the partitioning of pore space into wetting (brine) and non-wetting (dodecane) fluid. The higher range of the internal gradient spectrum is of the same order as tool gradients.

Nuclear magnetic resonance (NMR) response simulation
The spin relaxation of a saturated porous system is simulated by using a lattice random walk method [34,48]. Initially, the walkers are placed randomly in the 3D pore space. At each time 9 step i, the walkers are moved from their initial position to a neighboring site and the clock of the walker is advanced by τ i = 2 /(6D 0 ), where D 0 is the bulk diffusion constant of the relevant fluid, reflecting Brownian dynamics. We treat each random walk as the movement of a spin packet with initial strength M w (t = 0). At each time step i of length τ i , the strength of the walker is reduced by the survival probability S i , with the strength of the walker at time t = i τ i given by Here, for bulk relaxation and S s = 1 − ν for surface relaxation. For steps within the same fluid, S s = 1. The killing probability ν is related to the surface relaxivity ρ via where A is a correction factor of order 1 (here, we take A = 3/2) accounting for details of the random walk implementation [34,48]. This leads to So far we have discussed the classical approach of using bulk relaxation and surface relaxation processes. In addition, we model in detail the phase accumulation of the random walker in the internal magnetic field for the CPMG sequence. The evolution of the phase for a particular spin in relation to the Lamor frequency at the starting position ω 0 = γ B z (0) is given by [33] φ where τ j is a variable time step of the random walk controlled by the local diffusion coefficient and lattice spacing. This problem can still be conveniently solved with the same random walk techniques. Note that the Gaussian phase approximation is not invoked. The π and π 2 pulses of the acquisition sequence are considered instantaneous and either flip the sign of the accumulated phase or stop/start the accumulation of φ D . All random walk calculations presented in this paper were carried out using 10 6 random walks per parameter set on a 1000 3 voxel central domain of the x-ray-CT images with a grid resolution of = 1 10 a ≈ 290 nm (Bentheimer) and = 1 20 a ≈ 140 nm (Berea), where a is the voxel size. A 4 s CPMG echo train was acquired for each setting.

Restricted diffusion and diffusion-averaged internal gradients
In porous media, the diffusing spins encounter fluid-solid or fluid-fluid boundaries, which in our simulation are considered impenetrable-we assume that no magnetization transfer across fluid-fluid interfaces takes place. Such phase boundaries lead to a time-dependent diffusion coefficient, while at the same time leading to a better averaging of local fields. To reach an analytical expression for the magnetization decay, typically a homogeneous linear background gradient (1) and the Gaussian phase approximation and/or the concept of effective gradients are used [10]. In [32] a parabolic gradient is considered. Here, we directly use (15) and calculate the effective gradient of a particular random path in the following way, with τ = t N and where the sum extends over all terms for which r j = r 0 . In the limit of no diffusion or severely restricted diffusion, this effective gradient approaches zero. For partially saturated rock, we expect stronger restricted diffusion effects and better averaging of internal gradients for the fluid with the higher diffusion coefficient. The relevant time intervals where the magnetization vector precesses in the transversal plane are the encoding time δ for diffusion in the external gradient and the encoding time τ for internal gradients. While in experiments eddy currents and other artefacts need to be considered leading to finite but small values of δ 1 and δ 2 with δ 1,2 δ, numerically one can use δ 1 = δ 2 = 0, and therefore τ = δ. We show in figure 6 for Bentheimer sandstone and in figure 7 for Berea sandstone the distributions of g 2 eff averaged over the time intervals τ = 5 ms, = 80 ms and = 320 ms at different saturations for the case of vanishing external gradients. We see that for the fully brine-saturated case, the diffusional averaging of the internal gradients is more effective than for the dodecane-saturated samples. Furthermore, as seen for distributions of the local internal fields, the internal fields in the non-wetting dodecane phase are significantly lower. The shape of the distribution functions of |g eff | is markedly smoother than for |g b | distributions. Diffusional averaging is evident in both the fully saturated and the partially saturated samples. The effective gradient distributions of Berea sandstone are in good agreement with measured values for a Berea sandstone with higher susceptibility contrast (Berea 100, figure 4 in [26]).
To further illustrate our point about diffusional averaging, we calculate directly from the image the diffusion coefficient distributions where r 0 and r t denote the start and end positions of a random walk and t is time. We consider the same time intervals as above (τ = 5 ms, = 80 ms and = 320 ms), depicted in figure 8 for Bentheimer sandstone and in figure 9 for Berea sandstone. For the fully saturated samples, the diffusion coefficient distributions are centered on the bulk diffusion coefficient of brine (D 0b = 2.3 × 10 −5 cm 2 s −1 ) and dodecane (D 0d = 8 × 10 −6 cm 2 s −1 ). There is evidence for restricted diffusion, particularly for longer diffusion times. For the partially saturated samples, we see a significant shift in the diffusion coefficient distribution of the wetting brine phase, while the dodecane diffusion coefficient distribution essentially stays unchanged for the partial saturation range considered. At these partial saturations, there is little contrast between brine and dodecane in terms of diffusion coefficients. However, the product of Dg 2 b will be different because of the higher internal gradients in the water phase at partial saturations. For a 1D NMRderived diffusion spectrum, the diffusion spectrum would be weighted by T 2 relaxation. We now consider the full T 2 -D NMR simulation in the presence of restricted diffusion and internal gradients.

Diffusion-relaxation analysis
We show in figures 10 and 11 the NMR T 2 -D relaxation diffusion responses for the classical PGSE sequence using τ = 5 ms, ∈ {80, 320} ms and a CPMG spacing of t E = 0.4 ms. Using (2)-(4), we can calculate the values of C a , C b and C a,b . For our particular choice, we have C a = C b and C a,b = 2C a , e.g. the prefactors in (1) are of similar order. We see that for both sandstones the peaks in the 100% brine-saturated samples shift to smaller diffusion coefficients for the longer diffusion intervals , implying that restricted diffusion has a much stronger influence on the signal than internal gradients for the chosen at 2 MHz static field. We remark that we considered a sequence with constant δ, varying the strength of the applied gradient g a . A sequence with varying δ, e.g. using a permanent gradient, would give rise to stronger internal gradient effects.

Conclusions
We have presented a numerical technique to simulate and analyse T 2 -D diffusion spectra on the basis of x-ray-CT images, avoiding common approximations, that is able to deal with nonlinear internal gradients. The derived effective gradients are in good agreement with the literature. For the contrast of susceptibilities considered at low fields, the internal gradients did not impact strongly on the diffusion coefficient spectrum. It was assumed that the clay regions exhibit homogeneous susceptibility.  [49]. Note that for crude oils a distribution of diffusion coefficients and bulk relaxation times present in the complex fluid would cause it to follow the sloped line.
We have shown that the internal gradients are much larger in the wetting phase for typical susceptibility contrasts and can expect that at higher fields the difference in the internal gradient distributions would provide contrast to aid in the derivation of wetting properties and fluid configurations, as the variation in experimental parameters allows one to control this contrast mechanism. This is in addition to surface relaxivities, which contribute to an incoherent T 2 relaxation mechanism.
In this work, we considered dodecane as a non-wetting fluid, which has a well-defined diffusion coefficient. For many crude oils, a distribution and correlation of T 2 and D values make the fluids more distinguishable. We conclude that the numerical simulations are in general agreement with experimental techniques and support the interpretation of NMR experiments involving 2D inverse Laplace transforms.