Model Synthetic Samples for Validation of NMR Signal Simulations

Simulations of nuclear magnetic resonance (NMR) signal from fluids contained in porous media (such as rock cores) need to account for both enhanced surface relaxation and the presence of internal magnetic field gradients due to magnetic susceptibility contrast between the rock matrix and the contained fluid phase. Such simulations are typically focussed on the extraction of the NMR T2 relaxation distribution which can be related to pore size and indirectly to system permeability. Discrepancies between such NMR signal simulations on digital rock cores and associated experimental measurements are however frequently reported; these are generally attributed to spatial variations in rock matric composition resulting in heterogeneously distributed NMR surface relaxivities (ρ) and internal magnetic field gradients. To this end, a range of synthetic sediments composed of variable mixtures of quartz and garnet sands were studied. These two constituents were selected for the following reasons: they have different densities allowing for ready phase differentiation in 3D μCT images of samples to use as simulation lattices and they have distinctly different ρ and magnetic susceptibility values which allow for a rigorous test of NMR simulations. Here these 3D simulations were used to calculate the distribution of internal magnetic field gradients in the range of samples, these data were then compared against corresponding NMR experimental measurements. Agreement was reasonably good with the largest discrepancy being the simulation predicting weak internal gradients (in the vicinity of the quartz sand for mixed samples) which were not detected experimentally. The suite of 3D μCT images and associated experimental NMR measurements are all publicly available for the development and validation of NMR simulation efforts. Model porous media samples (composed of quartz and garnet sands with distinctly different densities, magnetic susceptibilities and surface relaxivities) were used to systematically assess NMR signal modelling of saturating fluids. Computation of 3D internal magnetic field map from digital rocks, followed by random walk simulations of the NMR signal of the saturating fluids. Good agreement between the simulated and experimental internal magnetic field gradient probability distributions. Model porous media samples (composed of quartz and garnet sands with distinctly different densities, magnetic susceptibilities and surface relaxivities) were used to systematically assess NMR signal modelling of saturating fluids. Computation of 3D internal magnetic field map from digital rocks, followed by random walk simulations of the NMR signal of the saturating fluids. Good agreement between the simulated and experimental internal magnetic field gradient probability distributions.

1 Introduction and Background 1 H nuclear magnetic resonance (NMR) is a non-invasive measurement technique applied to a wide range of opaque samples. Here the focus is on its application to fluid in rock systems, which enable the prediction of a range of petrophysical properties (e.g. porosity, pore size, permeability, and tortuosity) necessary for reservoir characterization (Arns et al. 2005(Arns et al. , 2007Freedman and Heaton 2004;Kenyon 1997;Straley et al. 1995;Yang et al. 2019). NMR T 2 signal relaxation measurements are commonly performed to this end in both the laboratory and during well logging in the field. Interpretation of such acquired NMR signal and its extrapolation to quantitative petrophysical property prediction is however far from trivial, not least due to the presence of additional internal magnetic field gradients arising from magnetic susceptibility differences between the fluid and the containing solid matrix. In order to improve the quantitative understanding of these complex relationships, detailed pore-scale 3D simulations of the NMR signal evolution, mostly based on molecular random walks, are a valuable tool.
NMR T 2 relaxation time measurements are performed using the Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence or variations thereof (Hürlimann 1998;Sun and Dunn 2002;Zhang et al. 2016). The measured T 2 relaxation time constant can be described as the effective sum of three contributions as follows (Mitchell et al. 2010): where D 0 refers to the self-diffusion coefficient of the contained fluid; is the gyromagnetic ratio of relevant nuclear spin ( 1 H); g is any internal field gradient and t e is the time interval between the refocusing r.f. pulses in the CPMG train. The first contribution is the intrinsic relaxation characteristic of the relevant fluid, T 2,bulk . This is generally either well known or readily measured; note however that rock matrix dissolution can result in the in-situ and time-dependant introduction of dissolved paramagnetic impurities (Connolly et al. 2017) that result in a significant reduction in T 2,bulk . The second contribution, S V , is a function of interaction between the fluid molecules and the rock surface. Knowledge of the surface relaxivity, , allows for determination of the surface to volume ratio, S V , of the pore space containing the fluid and thus an estimate of the pore size. This is only tractable if the third contribution, which results from diffusion of the fluid molecules through internal magnetic field gradients, is negligible. This is generally assumed via the use of weak magnetic fields (g scales with the magnetic field used (e.g. Lewis and Seland 2016)) which however also result in a reduction in signal-to-noise ratios (SNR) and/or the use of the minimum possible value of t e (Mitchell et al. 2010). An inherent assumption in Eq. 1 is that the T 2 value for the fluid in a rock core system can be described by single values of ρ and g, respectively. In reality, rocks present complex mineralogy with heterogeneous distributions of ferromagnetic and paramagnetic content (e.g. illite and chlorite) which challenge these assumptions (Elsayed et al. 2021). Consequently, the estimation of petrophysical parameters based on measured T 2 distributions can be corrupted. Initial attempts at simulating and measuring internal magnetic field gradients in porous media, focussed on model pore space geometries-these included both packed spheres (Audoly et al. 2003;Cho et al. 2012;Dunn 2002;Majumdar and Gore 1988;Song 2003) and periodic arrays of cylinders (Sen and Axelrod 1999) for which analytical solutions are readily available. One common approach to simulate internal magnetic field gradients is the dipolar sum method (Brown 1961;Drain 1962;Glasel and Lee 1974), in which the magnetic field is approximated from a summation of individual dipole fields usually associated with each voxel in the simulation lattice. Extrapolation of this method to geometrically much more complex digital rocks using 3D μCT images of the internal pore space of sandstone rock samples was performed by Arns et al. (2005) resulting in reasonable estimates of system permeability. However, rock surface composition heterogeneity was identified as being a major challenge to the accuracy of these permeability estimates (Arns et al. (2007)). Similarly, Benavides et al. (2017) showed that simulated NMR T 2 distributions with the assumption of a uniform surface relaxivity (and thus surface composition homogeneity) were inconsistent with corresponding experimentally measured T 2 distributions. Extrapolation to magnetic field gradients in 3D digital rocks (Arns et al. 2011) reproduced the expected trends in diffusion-relaxation behavior for a variety of fluids but no direct comparison against independently acquired experimental data was made. Connolly et al. (2019) used both finite element methods (FEM) and the dipolar sum method to predict the internal magnetic field gradients in six different sandstone digital rocksvery similar probability distributions were produced by these two simulation methods. However, when compared to experimental measurements, the simulations failed to predict the stronger magnetic field gradients observed experimentally (an example for the case of Idaho Grey sandstone is shown in Fig. 1). This was ultimately attributed to both the assumption of a single surface relaxivity and a single value of magnetic susceptibility, and hence the assumption of homogeneous rock composition.
Given this consistent lack of agreement between simulation and experiment, the NMR simulation methodology is tested here against model randomly packed porous media characterized by readily quantifiable values of surface relaxivity and magnetic susceptibility for the constituent particles. The resultant synthetic samples consist of mixtures of variable amounts of quartz and garnet sands respectively, with grain sizes sieved to be in the range from 100 to 150 μm. These model components are selected as they provide sufficient electron density contrast to allow for unambiguous solid-phase differentiation in μCT images of the samples. They also provide significantly different values of ρ and magnetic susceptibility such that the local NMR response will distinctly depend on local porous medium composition. Consequently, we are able to compare NMR-based simulations of the internal magnetic field gradients in these samples in detail with independent experimental measurements and hence test the veracity of the simulation methodology.

Material and Sample Preparation
The synthetic porous media samples used in this study were formulated from predominately quartz sand (featuring iron content < 1 wt%) and garnet sand. The quartz sand was sourced by Fisher Scientific whilst the garnet was sourced from GMA Garnet Group. The garnet datasheet indicated that it was predominately (> 96 wt%; iron-rich) almandines in composition. Both sands were sieved to a narrow size distribution between 100 and 150 µm so as to provide comparable pore size distributions. The densities of the quartz and garnet sands were gravimetrically measured to be 2.64 and 4.14 g/cm 3 respectively. Synthetic samples were prepared by vigorously mixing and packing the two sand components according to the bulk volume ratios reported in Table 1. The samples were then saturated with brine solution (97% de-ionized water and 3% sodium chloride, sourced from Sigma Aldrich) prior to NMR measurements. Post NMR measurements, the samples were set with resin to allow for minimal disturbance of their packing or microstructure during transport prior to 3D μCT imaging of their respective pore space micro-structures.

NMR Measurements
A comparatively low magnetic field Magritek Rock Core Analyser featuring a 53 mm inner diameter radio frequency (r.f.) coil tuned to the 1 H (which was exclusively detected from the saturation fluid used) resonance of 2 MHz was used to perform all NMR measurements on the range of samples at ambient condition (~ 25℃, 1 atm). T 2 measurements were performed using a standard CPMG (Carr-Purcell-Meiboom-Gill) pulse sequence as schematically shown in Fig. 2. A total of 32 logarithmically spaced values of t e from 100 μs (the minimum possible) to 20 ms were used for each sample-this variation of t e allows for determination of the internal magnetic field gradients in the sample via appropriate application of Eq. 1 ) and as detailed below. Each train of echoes was acquired for a total of 3500 ms with signal averaging for measurements on each sample performed eight times.

Digital Rock
3D digital rocks of the various samples were generated using μCT images acquired using a Zeiss VersaXRM-500 X-ray µCT scanner. Images were acquired using a voxel resolution of 3.93 µm over a total volume of 44.36 mm 3 (900 3 voxels). These 3D images were filtered and gated to identify the three phases present-pore space, quartz sand and garnet sand using visualization software, Avizo by Thermofisher Scientific. Phase allocation was unambiguous on account of the large electron density contrast between the three phases present-sample 3D images (a subsample of 450 3 voxels (5.55 mm 3 ) are shown in Fig. 3 for sample IDs 2-4 (Table 1). Results presented below were found to be independent of which 450 3 voxel subsample was selected.

Simulations
A detailed description of the 3D simulations in digital rocks as used here, of both internal magnetic field gradients and subsequently the NMR signal evolution (based on random walks), is contained in Connolly et al. (2019). A summary is provided here for completeness:  Table 1 The initial simulation step involves quantification of the internal magnetic field within the pore space of the discretized image of the sample. This requires the magnetic susceptibility of the two sands used as input parameters (the determination of which is discussed below), more specifically the magnetic susceptibility contrast, , between the saturation fluid and the respective sand. Calculation of the internal magnetic field is performed via the dipolar sum method, in which each solid voxel in the digital rock is effectively equated to a uniformly magnetized sphere. When the sample is subjected to an external static magnetic field, 0 the in the z-direction, the spherical grain voxels of radius, R , emanate a magnetic dipole moment, , that can be described in Gaussian units (Song 2003) as follows: The contribution to the magnetic field from a solid voxel i to an arbitrary pore space voxel can then be calculated using: where r is the distance between the source solid voxel i and the arbitrary pore space voxel, θ is the angle between them relative to the z-direction and z = 0 + .
The final magnetic field within the pore space (shown here for voxel, k) can then be determined by summing up the contributions from all solid voxels (a total number of N solid voxels in the remainder of the digital rock) as follows: where i z i − k is the magnetic field contribution of the i th solid voxel to this k th voxel positioned in the pore space. With the resultant computed 3D magnetic field map within the pore space, the internal gradient, , can then be calculated using Eq. (5), and then, readily converted into internal gradient probability distributions.
NMR T 2 relaxation was then simulated via a random walk, whereby 'walkers' that are representative of the diffusing water molecules that provide the NMR signal, were randomly assigned locations in the pore space of the digital rock. These walkers were then allowed to diffuse through the pore space within the digital rock via Brownian motion with the diffusion distance, l , after each time step, t s , in a random direction calculated using Eq. (6).
Entry into the solid material as a result of these diffusive steps was prohibited, however, time in these circumstances was still increased by t s .
The resultant evolution in the NMR signal, M(t), can then be described by three contributions (which are the three right-hand side terms of Eq. (1)): 1. The first contribution is the simplest; signal attenuation or magnetization decay M(t) M 0 due to transverse spin relaxation of the bulk fluid, T 2,bulk , as per Eq. (10): where M 0 is the NMR signal in the absence of any signal relaxation. 2. The second contribution due to surface relaxation was implemented with a "killing probability", , that the walker can lose its magnetization following a collision with the solid surface of the porous media. The killing probability is determined by the surface relaxivity parameter, , (as shown in Eq. (1)) via Eq. (11) (as derived in Bergman et al. 1995): 3. Signal attenuation due to diffusion in the internal magnetic field gradients (as simulated) can be included in the evolution of M(t) as follows: is the accumulated phase shift experienced by a walker after t n = n × t s , where B z,j is the local (calculated) magnetic field corresponding to the position of the walker in the simulation lattice and 0 is the external applied magnetic fields. This simulation approach has been validated in the work of Connolly et al. (2019) where each of these three contributions were rigorously tested using model simulation lattices and hence known analytical solutions for the evolution of the NMR signal. Mitchell et al. (2010) showed that the probability distribution of NMR T 2 relaxation times for the fluid in a saturated porous medium, H T 2 , can be indirectly accessed via Eq. 13, where P(g) is the probability distribution of internal magnetic field gradients in the sample pore space.

Internal Magnetic Field Gradient Analysis
(8) y j+1 = y j + l sin sin , (9) z j+1 = z j + l cos , where t = n × t e is the n th echo in the acquired CPMG train as shown in Fig. 2 and t e is the duration between echoes. Via 2D inversion, Eq. 13 allows for the extraction of the 2D probability space of H T 2 versus P(g) . The effective applicability of Eq. 13 is limited to the so-called asymptotic short-time (ST) relaxation regime. In this regime, the diffusion length ( l e ∼ √ D 0 t e ) is considerably shorter than the characteristic pore length scale, l s . For our measurements and simulations, the largest echo time was 20 ms; for the self-diffusion coefficient of the brine solutions (2.2 × 10 -9 m 2 /s) this corresponds to a length scale of ~ 6 microns which is considerably smaller than the pore sizes which will be of the order of 100 microns.
Solving for internal magnetic field gradient distributions, P(g) , as well as T 2 distributions, H T 2 , are generally ill-conditioned problems. The solutions from simple least square regression are typically unmeaningful due to the presence of experimental noise (Connolly et al. 2017). To mitigate this issue, a Tikhonov regularisation technique (Tikhonov 1963) is applied, where a penalty function, weighted by a smoothing quantity (the magnitude of which is regulated by a regularisation parameter) is implemented during the 2D data inversion process. It is critical to select an optimum regularisation parameter, this can be achieved via the generalized cross-validation method (Hollingsworth and Johns 2003).

Results and Discussions
The estimated values of ρ and v are initially detailed for the two sand components. These are required input parameters for the 3D NMR simulations of the various samples using the digital rocks as simulation lattices. The results of these simulations (probability distributions of the internal magnetic field gradients) are subsequently compared against independent experimental measurements.

Magnetic Susceptibility
The magnetic susceptibilities of quartz and garnet sand were measured using the vibrating sample magnetometer (VSM) of a Quantum Design Physical Property Measurement System (PPMS-9) (Mohammed et al. 2020). The vibration frequency was 40 Hz and the amplitude was 2 mm. The magnetic polarization, M of each sample was measured as a function of external magnetic field strength, H (ranging from − 5 to 5 T), and can be described by Eq. (14) for diamagnetic and paramagnetic materials: where v is the required magnetic susceptibility.
The magnetic polarization of quartz and garnet sand samples were measured at the full range of available magnetic field (− 5 T to 5 T and then back to − 5 T). No hysteresis was observed, suggesting that the samples contain negligible amounts of ferrimagnetic or ferromagnetic materials. The magnetic field of interest at which all NMR measurements were performed is 0.05 T. Hence the value of magnetic susceptibility for both the qutz sand and = v garnet sand was determined from magnetic field data ranging from -0.1 to 0.1 T. This is shown in Fig. 4; the magnetic susceptibility values thus determined were 5.9 × 10 -6 and 2.5 × 10 -3 (SI) for the quartz sand and garnet sand respectively. These data can be verified via a semi-empirical approach as detailed by Lewis and Seland (2016). Equation (15) relates the estimated mean pore size, a, of the sample to the magnetic susceptibility contrast . and the mean internal magnetic field gradient, G.
We estimate G for each pure sand sample below from the experimental measurements (as detailed in Sect. 3.3: Internal Magnetic Field Gradients). The mean pore size was estimated based on the known packing size and porosity; this was estimated to be ~ 50 microns. Via this independent approach, the effective magnetic susceptibilities (SI) of quartz sand and garnet sand were estimated to be 7.2 × 10 -6 and 4.1 × 10 -3 , which given the crude assumptions involved are consistent with the magnetometer measurements detailed above.

Surface Relaxivity
Determination of the surface relaxivity value for the two sand components was performed via simulation of the CPMG echo train for the minimum echo time, t e (100 μs). These data are then regularised to produce T 2 distributions as a function of the surface relaxivity, ρ for each pure sand sample. The value of ρ which minimizes the difference between these simulated distributions of T 2 and the experimentally determined distributions of T 2 is thus determined. The optimum fits for the two sands are shown in Fig. 5. In both cases, agreement is excellent and resulted in ρ values of 12.5 µm/s and 98.5 µm/s for the quartz and garnet sands respectively. The more distinct peaks for the simulation data in Fig. 5b reflect the inherently better SNR for the simulations relative to the experimental data. The above assessment assumes no contribution from the background gradients at the minimum echo time (100 μs). To test this, we simulated the T 2 distributions for a sample of 25% quartz sand and 75% garnet sand (sample 4 in Table 1) for the minimum echo time with and without the effect of internal magnetic fields included. The resultant T 2 distributions are shown in Fig. 6. The simulation result shows that the impact of internal gradient is effectively negligible at this minimum value of t e and produces excellent agreement with the independently measured experimental T 2 distribution.
For validation purposes, we estimated the S/V ratio of the pore space from the digital rock images for the two sands using the maximal ball algorithm (Al-Kharusi and Blunt 2007; Silin and Patzek 2006;Silin et al. 2003) to extract the pore and throat network from the 3D digital rock model of each sample. With the pores being approximated as spheres, this allowed ready estimation of the mean pore space S/V ratio. By adjusting ρ, the experimentally measured values of T 2 for these samples were matched to the relevant mean pore space S/V values producing estimates of ρ of 22 µm/s and 80 µm/s for the quartz and garnet sands respectively. Given the geometrical assumptions used, these values are consistent with the measurements reported above. Fig. 6 T 2 distributions of the saturation fluid measured sample 4 (Q25G75) as predicted with and without the consideration of internal magnetic field gradients and as measured experimentally. All distributions correspond to measurements or calculations where the minimum echo time was used 1 3

Internal Magnetic Field Gradients
In this work the internal magnetic field gradient probability distributions were determined via three methods: • The first simulation-based method was to calculate these distributions directly based on the 3D magnetic field maps generated via the dipolar method (as detailed above) for the pore space of the various digital rocks of the samples. This is hereafter referred to as DipSum_IGD. • The second simulation-based method was to perform random walk simulations in these 3D magnetic field maps for the various digital rocks, thus producing CPMG signal decay data series for the range of echo time considered. The resultant 2D data set was inverted as per Eq. 13 to produce the 2D distribution of internal magnetic field gradients against T 2 . By collapsing the T 2 dimension, the 1D distribution of internal magnetic field gradients was determined. This is hereafter referred to as RW_IGD. • The third experimental method focussed on the acquired CPMG experimental data as a function of echo time. This was also inverted using Eq. 13 to produce the distribution of magnetic field gradients in one of the dimensions. This is hereafter referred to as Exp_IGD. Figure 7 shows sample 2D slices from each of the 3D digital rocks for each of the five samples (as per Table 1). In each case, a segmented image is shown in which the two constituent sands are readily identified followed by the equivalent image in which the computed internal magnetic field gradients using the dipolar sum method are shown in the pore space. As the garnet sand content increases, stronger magnetic gradients are clearly evident across the sample as would be expected. Local stronger gradients are clearly evident in the immediate vicinity of the garnet sand particles with values up to 2 T/m evident (reflecting the greater magnetic susceptibility of 5.1 × 10 -3 for garnet sand versus 7.2 × 10 -6 for the quartz sand).
The probability distributions of internal magnetic field gradients as predicted by the three methods (DipSum_IGD, RW_IGD & Exp_IGD as detailed above) are shown for each sample in Fig. 8. There is broad good agreement across the samples for the independent determinations of these distributions between the simulation-based and experimental methods. In the case of 100% quartz sand (Fig. 8a), two discrepancies are evident. A small population of stronger gradients (> 0.5 T/m) is evident in the experimental data that is not replicated in the simulation data. This we believe reflects a very small amount of heterogeneously distributed ferric content in the quartz sand, which was absent in the case of the simulations which assumed a homogeneous composition. In terms of the population at weaker gradients (< 0.1 T/m), the agreement regards the modal value is excellent; however, the experimental distribution is broader than either simulation. These weak gradients result in minimal change in the NMR signal as a function of echo time (t e ), thus these data are more heavily influenced by noise which is present in the experimental data to a much greater extent than in the simulation data. Regularisation will produce broader peaks for noisier data. In the case of mixed sand samples ( Fig. 8b-d) agreement is excellent in terms of the stronger magnetic field gradient values. Where a significantly weaker magnetic field gradient peak is present in the simulation data(as in Fig. 8b) this is largely not captured by the experimental data; this peak is associated with NMR signal from the immediate porescale vicinity of the quartz sand particles. The origin of this discrepancy potentially lies in the ability of the experimental data (in the presence of experimental noise) to cater for three orders of magnitude variation in internal magnetic field gradient strength. The evolution in NMR signal with echo time is dominated by the much stronger magnetic field gradient peak, with any variations of this signal due to the much weaker magnetic field gradient peak being more readily compromised by experimental noise.
An alternative explanation for the discrepancies might be diffusive signal averaging where the molecules diffuse during the duration of the CPMG train between regions dominated by either quartz or garnet sand. However, such diffusion effects would also be present in the RW_IGD simulation data, and thus this potential contribution to the discrepancy can be discounted. Nonetheless, agreement between the simulations and experimental data is comparatively good, particularly relative to the literature with significantly stronger gradients being simulated unlike in Fig. 1 above. Collectively this data (range of digital rocks with clear phase identification coupled with experimental CPMG data as a function of echo time) is openly available from figshare (https:// doi. org/ 10. 6084/ m9. figsh are. 17161 730. v1) Fig. 8 Internal gradient distributions obtained from the three methods (labelled as DipSum_ IGD, RW_IGD & Exp_IGD) for a sample 1 (Q100), b sample 2 (Q75G25), c sample 3 (Q50G50), d sample 4 (Q25G75) and e sample 5 (G100) and is promoted here as a suitable model porous media system for systematic testing of NMR signal response simulation codes for porous media.

Conclusions
Model porous media samples consisting of random packings of quartz and garnet sands are studied as an effective means of systematically assessing and validating NMR signal modeling for fluids contained in porous media. These two constituent sands present a significant difference in electron density (allowing for good phase differentiation in μCT images to generate digital rocks) as well as both surface relaxivity and magnetic susceptibility, thus allowing for a range of NMR responses to be probed. Comparisons of the simulated and experimental distributions of magnetic field gradients were generally in good agreement with contributions from stronger local magnetic fields being adequately captured by the simulations. The discrepancy between the simulated and experimental results is significantly reduced in comparison with that shown in the literature (Arns et al. 2011;Chen et al. 2005;Connolly et al. 2019), this serves to suggest that the much greater compositional and microstructural heterogenoeity of the real rock samples considered in these works is responsible for the greater discrepancy observed and thus not the implemented simulation codes and their required assumptions. Discrepancies existed when significant weaker magnetic field gradients were present in combination with stronger magnetic field gradientsthese were not adequately captured by the effective detection range of the experimental data, which was attributed to sensitivity to very minor NMR signal perturbation by comparatively weak internal magnetic field gradients as a function of echo time in the presence of comparable experimental noise. All data for these model systems are publicly available and are recommended as suitable model systems for NMR simulation code validation and comparison. permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.