NMR Relaxation Modelling in Porous Media with Dual-Scale-Resolved Internal Magnetic Fields

A pore-scale forward modelling approach for NMR relaxation responses of sandstones incorporating their dual-scale nature is presented. The approach utilises X-ray micro-CT images to capture inter-granular porosity and scanning electron microscopy images to reconstruct clay regions via a resolved clay micro-structure model. A key to calculating the NMR response with resolved clay micro-structure is the development of a dual-scale internal magnetic field calculation. This is achieved by a separation of near- and far-field effects in a dipole approximation of the internal field with periodic clay micro-structures, the latter of which take local clay pocket porosity into account. Tri-linear interpolation of the micro-CT image before calculation of the internal magnetic field further reduces errors in the transition regions between coarse- and fine-scale structure, with final discretisation level matching the fine-scale clay micro-structure model across the whole domain. The method is validated against direct calculations of model media at full resolution and applied to Bentheimer sandstone. Measured and simulated NMR T2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_2$$\end{document} relaxation responses, including relaxation time distribution shape, are in excellent agreement and distributions of internal magnetic field gradients at the highest spatial resolution as well as diffusion-averaged effective gradients are reported.


Introduction
Due to advances in micro-structure characterisation and computational techniques, the modelling of pore-scale transport phenomena is increasingly important in ground water engineering and geostorage applications. Increasing the accuracy of rock representations would significantly improve the quality of forward simulations, which are increasingly important, e.g. in the context of carbon capture and storage. Reservoir rocks including sandstone (Øren and Bakke 2002) and carbonate rocks (Biswal et al. 2007) typically exhibit dual-scale micro-structures. Capturing multiple scales of porosity in rocks, like vugs, macro-, and micro-porosity, is a long-standing problem: while nano-computed tomography (nano-CT) and focused ion beam scanning electron microscopy (FIB-SEM) offer exceptional resolution of micro-structures, this comes at the expense of a very narrow field of view (FoV), typically missing larger-scale heterogeneity of the pore space entirely (Bai et al. 2013;Welch et al. 2017). X-ray micro-computed tomography (micro-CT) is capable of capturing heterogeneity at core scale but is limited to micron resolution (Arns et al. 2003). The problem can be partially overcome by the introduction of so-called superresolution techniques (Chang et al. 2004;Glasner et al. 2009;Wang et al. 2018b), which increase the effective resolution by a factor 4 to 8 with the help of auxiliary high-resolution information for training. This approach is insufficient to accurately represent the microstructure of clay regions when applied at X-ray micro-CT imaging resolution. Effective medium theories (Bruggeman 1935;Kirkpatrick 1971) describe the macroscopic behaviour of a micro-structure at larger scales, avoiding the need for explicitly resolved spatial information. A combination of the latter with direct imaging in the context of digital core analysis was proposed in the past (Arns et al. 2005). "Macroscopic" pore-scale spatial information is directly captured by micro-CT imaging of an appropriate resolution, and an effective medium representation is used for the finest elements of the system (micro-porosity, clays, etc.), including in recent upscaling approaches (Jiang and Arns 2021). However, this strategy has a significant limitation-every effective element or region needs to be large enough to possess the effective property of interest. If this condition is not, fulfilled problematic boundary conditions between regions of resolved and non-resolved porosity may arise. Another convenient approach for representing multiple spatial scales of rocks is pore-network modelling (Øren and Pinczewski 1994;Stig and Øren 1997). Due to strong abstractions of the local geometry and the need to incorporate mineralogy, it is less suited for the predictive modelling of NMR responses. NMR transverse relaxation time ( T 2 ) distributions are frequently used in reservoir characterisation to estimate pore size distributions and permeability (Kleinberg et al. 1994). However, T 2 responses are sensitive to internal magnetic field inhomogeneities induced by susceptibility contrast between solids and fluids in rocks. Molecular diffusion in the presence of strong internal magnetic field gradients leads to additional dephasing, which enhances transverse relaxation, potentially invalidating the condition of 'fast diffusion' (Brownstein and Tarr 1979). Accounting for relaxation enhancement due to internal magnetic field gradients in pore-scale simulations requires the accurate calculation of these fields.
Over the past decades, two main numerical methods for deriving the internal magnetic field distribution in rocks have been established. The first method derives the magnetic field distribution by calculating the magnetic potential using the finite element method (FEM) to solve the Maxwell equations in zero-current condition. For instance, the internal magnetic field and internal gradients are simulated using thin sections of Berea sandstone (2D model) for multiple constituent phases (Chen et al. 2005). The 3D internal magnetic field of a single pore is calculated in Fordham and Mitchell (2018). More recently the 3D internal magnetic field was computed for digital rock samples (Connolly et al. 2019), assigning a mean magnetic susceptibility across all solid components.
The second method to calculate the internal magnetic field is based on a linear superposition of dipolar fields (also known as dipole approximation). Internal field distributions of random sphere packs were calculated (Drain 1962;Majumdar and Gore 1988;Audoly et al. 2003). The approach has been investigated for a periodic cylinder pack system, which represented the simplest pore geometry for visualisation of internal field distribution (Sen and Axelrod 1999). The application of the dipole approximation to 3D micro-CT images of rocks including for multiple mineral phases was demonstrated in the past (Arns 2004;Arns et al. 2007Arns et al. , 2011. However, the magnetic field distribution of the clay phase was calculated assuming a continuous effective medium with constant volume-weighted susceptibility for the entire clay region due to a lack of image resolution. Utilising directly a high-resolution micro-structure of clay regions poses a significant computational challenge. Such a micro-structure could be chosen to be periodic in a way that it still covers regional micro-structure patches (Arns et al. 2009). The latter resulted in local porosity variations and, therefore, average susceptibility variation at voxel scale. Calculating highresolution internal magnetic fields then must consider near-field effects at intermediate scales. We note that resolving the micro-structure of kaolinite requires resolution of the order of 20 nm (Cui et al. 2021). A direct approach using, e.g. an n 3 c voxel sub-domain, n c = 1000 , of a micro-CT image at 2 μm resolution would result in fields of size 100,000 3 , which is impractical.
Most NMR response simulations are based on random walk algorithms, where each 'walker' is representative of a nuclear spin package (Mendelson 1990;Bergman et al. 1995). 2D anisotropic pore geometries involving different pore shapes were considered by Wilkinson et al. (1991). Comparisons of simulation and experiment for the mono-sized grain consolidation model led to good agreement with experiments (Straley and Schwartz 1996). However, internal magnetic field gradients were not included in these earlier studies. The dephasing spins are considered in a uniform gradient in a porous material (Valckenborg et al. 2002) and for ferromagnetic grains (Valckenborg et al. 2003). The NMR response is simulated with constant gradient based on a 3D model that includes arbitrary bi-modal pore size distributions by compacted spheres (Toumelin et al. 2003). Consider now NMR relaxation response simulations incorporating a general internal magnetic field. Random walk simulations on X-ray micro-CT images of heterogeneous sedimentary rock samples comparing different definitions of pore size were presented in Arns (2004). D-T 2 responses of Bentheimer and Berea sandstone were modelled based on micro-CT images including internal field effects (Arns et al. 2011). Recent work simulating T 2 responses led to partial agreement with experiments; short relaxation times are not recovered (Connolly et al. 2019).
Another method for NMR relaxation response simulation is based on the FEM to solve the Bloch-Torrey equations for the evolution of the magnetisation. The magnetisation decay due to diffusion in the internal magnetic field gradient in a single pore of complex shape was explored in Fordham and Mitchell (2018). A recent study (Mitchell et al. 2019) involving intra-and inter-particle porosity compares random walks and FEM on simplified models. Due to different magnitudes of length scales in micro-and macro-pores, microporous grains were treated as effective medium.
The increasing sophistication of NMR simulators for porous media responses opens new applications regarding the extraction of intrinsic physical properties of complex materials when a high-resolution 3D physical representation is available, e.g. by micro-CT imaging. Recently, a comprehensive and efficient Bayesian optimisation framework to extract multiple intrinsic physical parameters simultaneously was introduced . They demonstrated the framework on Bentheimer sandstone, extracting quartz surface relaxivity as well as effective diffusion and transverse relaxation time of unresolved clay regions.
The latter suggests that dynamic changes in surface properties, e.g. as result of wettability alteration (Shikhov et al. 2018(Shikhov et al. , 2019, may be understood spatially resolved and inform modelling of two-phase flow (Wang et al. 2018a). For such an approach, it would be desirable to replace effective clay parameters with clay structure and surface properties. A first step in this direction has been presented by Cui et al. (2021), where a dual-scale model of Bentheimer sandstone described the heterogeneity of the rock including micro-porosity in clay regions. However, a multi-stage approach was used for the simulation of the NMR responses.
In this work, we present a dual-scale NMR solver which supports multiple clay microstructures and accounts for the resultant internal magnetic fields at the highest resolution. The paper is organised as follows: we first introduce the materials including a simplified three-layer model, kaolinite model, and the representation of Bentheimer sandstone. This is followed by the dual-scale internal field calculation utilising the dipole approximation, based on a decomposition of the field into high-and low-resolution components and accounting for near-and far-field effects. We then introduce the dual-scale random walk approach to solve for NMR responses. Section 3 illustrates the accuracy of the simulation approach utilising the model material, and applications to Bentheimer sandstone with excellent match to experimental measurements. Finally, we report on the distributions of spatial internal magnetic field gradients and diffusion-averaged effective gradients for different regions of the sample.

Materials and Methods
The study includes dual-resolution samples which consist of two categories. One is the coarse-resolution micro-CT image (around 2 to 4 μm resolution). The other is the fineresolution artificially reconstructed 3D model (down to about 10 nm resolution). There can be N different fine-scale models with different porosity and micro-structure representing coarse-scale (unresolved) phases in the low-resolution image. The internal fields of these domains will be named as coarse field and fine field in the later discussion, with transition regions giving cause to near-field effects. The simple three-layer model is first introduced here to better explain the dual-scale internal field calculation methodology. We then transfer to complex rock geometry.

Three-Layer Model
The coarse-resolution sample ( Fig. 1a) contains three parallel layers, assumed, respectively, as pore (top), unresolved "clay" phase (middle) which will be replaced by resolved micro-structure, and solid (bottom). Image size is 20 × 20 × 40 cubic voxel with resolution of 3 μm/voxel. The 3D micro-structure is constructed by a level cut through a Gaussian Random Field (GRF). For such fields, the resultant two-point correlation function can be calculated analytically and fitted to measured correlation functions, e.g. from 2D highresolution SEM micrographs. We employ the field-field correlation function (Cahn 1965;Roberts 1997;Arns et al. 2004): with cut-off scale r c = 0.7355, correlation length = 0.9095 and domain size d = 7.6434 in voxel units. Actual values are given for concreteness but have no particular meaning other than defining the model structure. With the known g(r), realisations of GRF Z( ) are numerically generated using the fast Fourier transform (FFT) method utilising a one-level cut with level cut parameter = 0.17. The fine-scale model structure is generated with periodic boundary conditions as 400 3 voxel domain with the model domain covering 8 3 voxels of the coarse-resolution image (Fig. 1b). To generate a fully resolved model, we map the coarse-resolution sample to a conforming subgrid. Each coarse-scale voxel corresponds to n 3 s fine-scale voxels. Here, we chose n s = 3 μm∕0.06 μm = 50 . This results in a 1000 × 1000 × 2000 voxel micro-structure at 60 nm resolution, where the unresolved middle layer of the structure is replaced by corresponding micro-structure voxels. The fully resolved model is visualised in Fig. 1c. Note the periodic repetition of the template 400 3 domain.

Bentheimer Sandstone
Bentheimer sandstone exhibits a more complex geometry and is often used as benchmark rock (Arns et al. 2011;Mitchell et al. 2013;Peksa et al. 2015). While Bentheimer sandstone does not contain a large amount of clay (around 1.9 wt% , the clay type is mainly kaolinite and the reason we chose it for this study. An X-ray micro-CT image was acquired at resolution of 2.2 μm/voxel. More details about the specific sample and image acquisition can be found in our previous work . The simulation presented here is based on an 800 3 central subset of the original image data. The sample is segmented into five phases, namely pore space, clay region (mixture of pore space and kaolinite),

Kaolinite Model
Kaolinite is one of the most common clay mineral constituents of sedimentary rocks, especially in sandstones. It is difficult to capture the booklets' morphology of kaolinite ( Fig. 3a, b) in the micro-CT image due to resolution limitations. We employ a 3D reconstruction technique (Cui et al. 2021) based on SEM images to overcome these limitations. Kaolinite booklets are distributed following a Boolean process with random unit rotation angles (Cui et al. 2021). Individual kaolinite platelets as constituents of the booklets are approximated as hexagonal plates with thickness of about 60 nm and diameter of 2-3 μm. To mimic the morphology of booklet units, we generate stacked hexagonal platelets, separated by a small separation distance (about 30 nm). In this study, we choose a domain size of n 3 f voxel, n f = 400 , with resolution 22 nm, which is 1/100 of micro-CT image resolution, to discretise the kaolinite micro-structure (Fig. 3c, d). The surface area ( S v ) of the model is 18.7 m 2 /g in agreement with the literature (East 1950).

Internal Magnetic Field
The magnetic field at position consists of the uniform external magnetic field 0 and the internal magnetic field i induced by the magnetic susceptibility contrast: The dipole approximation is employed to calculate the nonuniform internal magnetic field (e.g. Song 2003); SI unit will be utilised. The dipole moment of an isolated sphere of radius r D as function of distance r to the centre of the sphere is given by Fig. 2 a A slice of 800 3 tomogram of Bentheimer, 2.2 μm/voxel. b A slice of 800 3 5-phase segmented tomogram of Bentheimer, 2.2 μm/voxel. The phases colour scheme: pore space-blue; clay-light blue; quartz-green; feldspar-brown; high-density minerals-red where 0 = 4 × 10 −7 NA −2 is the magnetic permeability of the vacuum. The magnetic dipole moment of a single cubic voxel is calculated as ≈ 0 v 0 for v ≪ 1 , which is a good approximation for typical rock constituent minerals (Hürlimann 1998). The susceptibility field of a larger body is generated by assigning each voxel an isotropic magnetic volume susceptibility v ( ) according to phase occupancy at location . The induced change in magnetic field is linear in susceptibility, and the internal magnetic field can be calculated as a sum of dipoles. Then the internal magnetic field is given by the convolution of dipole and susceptibility field: The approach is validated by deriving the dipole profile of a sphere discretised on a Cartesian regular grid with spacing = 1 15 r D and comparison with the analytical solution. Figure 4 depicts cross sections of the resultant fields and profiles along the z-axis through the centre of the dipole generated by a solid sphere inclusion ( s = 5.0 × 10 −6 ) in water ( w = −9 × 10 −6 ) with applied field 0 = (0, 0, B 0 ).
Consider now the case of the three-layer model, which consists-at the coarse scale (Fig. 1a)-of one unresolved and two resolved phases. While the resolved phases directly take the magnetic susceptibilities w and s , each unresolved coarse-scale voxel maps to n 3 s fine-scale voxels, a sub-domain of the fine-scale model as discussed previously. Consequently, the coarse-scale susceptibility of the unresolved phase u must vary locally (unless n 3 s is the same size as the fine-scale model), and those voxels are assigned the local volumeweighted susceptibility where c stands for the volume susceptibility of pure clay. Convolution of the coarse-resolution susceptibility field v ( ) with the dipole field, where v ( ) takes the values s , w , or u ( ) , gives the coarse-resolution internal magnetic field i,c . The convolution establishes a macroscopic internal field trend in the unresolved phase as well as near-field perturbations in the resolved phase in regions close to the fine-scale structure.
The fine-scale field i,f is calculated analogous to i,c , using the constituent material susceptibilities w and c . From i,f we derive the coarse-scale averaged trend model i,t covering the fine-scale domain. For our periodic fine-scale structure with n f = 400 (Fig. 1b) and n s = 50 the periodic trend model has a size of n 3 t voxel with n t = n f ∕n s = 8 . In practice, we can map i,t to high resolution by tri-linear interpolation. At places where no fine-scale structure exists because the coarse-scale voxel belongs to a fully resolved phase, the associated trend value and fine-scale field are set to zero. A high-resolution internal magnetic field i covering the whole simulation domain of size (n c n s ) 3 voxel can now be calculated as follows: We note that the field i is never stored in memory but constructed locally when required for either a random trajectory during NMR simulations, to visualise parts of the high-resolution internal magnetic field, or to extract statistics. Given i ≡ (0, 0, B z ) , the corresponding internal magnetic field gradients at image resolution are given by and are calculated using central differences. In Sect. 3, we report histograms and visualisations of the internal gradient fields for different phases with the further distinction of including ( G i ≡ G in ) or excluding ( G i ≡ G ex ) voxels which have different phase neighbours. Not including gradients for voxels with different phase neighbours is reasonable, since near-surface interactions are included in surface relaxivity.
The accuracy of the internal field representation can be further improved by utilising a higher-resolution mesh at the coarse scale. In Sect. 3, we make a distinction between results for a fully resolved calculation where a single fine-scale grid is utilised (possible only for the artificial model), and for coarse-grid refinement factors between 1 and 5. While a factor of 1 indicates that the initial resolution of the coarser grid is utilised, it represents nevertheless a result for the dual-scale approach.

NMR Response Simulation
The NMR T 2 relaxation response is simulated by a lattice random walk algorithm on the segmented image (Arns et al. 2011). In this study, we extend the method to dual-scale samples. Initial positions of the random walkers are chosen uniform randomly across the pore space. Since the dual-scale approach assigns either solid (quartz or kaolinite) or pore space to each fine-scale voxel, we use a simple rejection method. The random walk is carried out on a regular cubic lattice with resolution . Consequently, the time step i of a step i is given by i = 6D( )∕ 2 , where D( ) is the local diffusion coefficient. Since here the structure is fully resolved and we consider a single fluid, D( ) = D 0 , where D 0 is the free-diffusion coefficient of water. Bulk and surface relaxation are implemented as signal weighting factors for each step, with S i = S b S s . Here S b = exp(− i ∕T 2,b ) is the contribution by bulk relaxation, and S s = 1 − 6 2 i ∕( A) the contribution by surface relaxation. For steps within the same fluid S s = 1 . A is a correction factor of order 1, which accounts for the details of the random walk implementation (Bergman et al. 1995). For imaged structures, this value is typically close to A = 3∕2 . For Bentheimer sandstone, the surface relaxivity 2 of the solid takes the value of quartz ( 2q ), feldspar ( 2f ), high-density minerals ( 2h ), or kaolinite ( 2k ).
While 2 ( ) describes the relaxation processes at the surface, the contribution of the internal field to dephasing is calculated explicitly. The evolution of the phase for a spin package in relation to the Larmor frequency at the starting position 0 = B z (0) can be written as follows (Carr and Purcell 1954): where is the gyromagnetic ratio. 0 and j note the location of the start of the random walk and at time t j . The Carr-Purcell-Meiboom-Gill (CPMG) sequence (Carr and Purcell 1954;Meiboom and Gill 1958) is implemented by switching the sign of Δ after t E ∕2 , where t E is the echo spacing, and every nt E , n = 1, 2, … N − 1 thereafter, where N is the total number of echoes recorded (at echo positions t E ∕2 after the phase switch). The phase variation caused by this diffusive relaxation is then incorporated into the calculation of each walker's magnetisation, which becomes Averaging over random walks and recording at times t = nt E ( n = 1, 2, … N ) results in the magnetisation decay M(t). The relaxation time distribution is obtained by fitting a multiexponential decay to M(t). The effective gradient of a particular random path of N s steps can be calculated as follows: with N s = ∑ N s i=1 i and where the sum extends over all terms for which j ≠ 0 and may be considered most meaningful for 2 N s ≤ t E . Thus, G ef f refers to a diffusion-averaged (effective) gradient over part of a magnetisation decay. We later report results for G ef f ( 1 2 t E ) , avoiding refocusing pulses.

NMR Experiments
Bentheimer sandstone core plugs of 25.4 mm diameter and 50 mm length were saturated with 3% NaCl brine and transverse relaxation responses measured with a 2 MHz Magritek Rock Core Analyzer. The standard CPMG pulse sequence was utilised with echo-train length of 10 seconds. Measurements were recorded for two echo times ( t E = 200 μ s and t E = 640 μ s) to study the significance of internal magnetic field gradients. The signal-to-noise ratio (SNR) was set to at least 200. The acquired magnetisation decays were treated as multi-exponential signals and inverted into T 2 relaxation distributions. Like for the NMR simulations, the inversion of the NMR relaxation decay into the relaxation time distribution was carried out utilising a non-negative least square approach with amplitude regularisation term (Hansen 1992).

Results and Discussion
The introduced internal field calculation method and associated NMR response simulation approach are demonstrated for the three-layer model first. We illustrate the field distribution change in the transition region with different degree of discretisation and the accuracy of the dual-scale model. Applications to Bentheimer sandstone including a discussion of internal field distributions follow:

Three-Layer Model
To verify the dual-scale B i calculation method, we calculate the internal magnetic field B i at full resolution directly for the 1000 × 1000 × 2000 resolved micro-structure (Fig. 1c) and compare it to the dual-scale approach with coarse-grid refinement factors of 1, 2, and 5 in Figs. 5 and 6. For concreteness, we present the fields for the common specific field strength of B 0 = 470 Gauss corresponding to a frequency of f = 2 MHz for 1 H NMR. The respective susceptibilities are given in Table 1; we assign a relatively high but not implausible paramagnetic susceptibility for pure clay to highlight the effects from internal fields. The B i cross sections (Fig. 5) illustrate the overall closeness of the dual-scale approach with the full-scale calculation as well as the further improvement of the solution for the dual-scale approach when using grid refinement for the coarse-scale model, particularly in the transition region as visible in the rectangular cut-outs of Fig. 5. A grid refinement factor of 5, resulting in a coarse grid of 100 × 100 × 200 voxel with n s = 10 leads to excellent agreement. This is further illustrated in Fig. 6a for a line crossing the coarse-scale phase boundaries perpendicularly.
We immediately note the large fluctuations in field associated with the fine-scale microstructure, which is fully recovered by all grid refinement factors, as is the far field. From Fig. 6b, magnifying one of the transition regions, we see that the spatial extend of the difference between full and dual-scale approach close to the transition region drops drastically with increasing grid refinement. We note that the resolution of the micro-structure is exactly the same for all grid refinement factors here due to the alignment of the coarsescale phase distribution with the grid. The main effect of the grid refinement is that the coarse-scale trend model of the fine-scale micro-structure has a higher resolution, which impacts on the accuracy of the convolution calculating the coarse-scale field contribution i,c to i (see Eq. 5).
Consider now the simulation of NMR T 2 relaxation responses utilising the dual-scale internal field calculation approach for the three-layer model. NMR simulation parameters are shown in Table 1. Since internal field effects are stronger for higher field strength or larger echo spacing (Kleinberg and Horsfield 1990), we calculate T 2 distributions for the following four conditions: (1) f = 2 MHz, t E = 200 μ s; (2) f = 2 MHz, t E = 640 μ s; (3)  Fig. 7. For all chosen field strength and echo time combinations, the position of the main relaxation time peak is preserved, while the sharpness of the main peak is affected by the internal field discretisation. However, there are significant differences in the width and position of the short relaxation time peak. We note that a coarse-grid refinement of a factor five leads in all cases to excellent agreement between the full-resolution solution and the dual-scale simulation approach. We further remark that the dual-scale approach mimics well the apparent "coupling" of the two main peaks for the case of f = 20 MHz and t E = 640 μ s. The discretisation level of the internal field (recall that the fine-grid resolution remains constant) becomes increasingly important with increasing field strength and longer echo times.

Bentheimer Sandstone
The segmented micro-CT image of Bentheimer (coarse resolution, = 2.2 μ m) is combined with the kaolinite model (fine resolution, = 22nm) to construct the high-resolution internal field B i as well as the high-resolution micro-structure. However, jagged edges along phase boundaries defined by cubic voxels are magnified when refining the coarseresolution sample. To reduce these artefacts, we use a surface interpolation technique to obtain smoother solid surfaces and consider refinement factors of 1, 2, and 4 (Fig. 8). This results in coarse-grid sizes of 800 3 , 1600 3 and 3200 3 , respectively. The fine-scale grid is constructed at a resolution of = 22 nm with a fine-scale (logical) grid size of 80, 000 3 and n s = 100 , n s = 50 , and n s = 25 for the refinement factors of 1, 2, and 4. NMR T 2 response simulations are performed for two different echo spacings, t E = 200 μ s and t E = 640 μ s in order to investigate the significance of internal magnetic field gradients. Here, we only consider a field strength of 470 Gauss, for which we have experimental data. In an earlier study (Cui et al. 2021), we determined the effective parameters for the unresolved clay phase of Bentheimer sandstone by matching simulated NMR responses to experimental data. In this study, the surface relaxivities utilised in simulations were slightly adjusted compared to the previous study (Cui et al. 2021) to account for the integration of the internal field and micro-structure model into a single dual-scale simulation. Table 2 lists the utilised material properties and 2 for each phase. Susceptibilities were set to typical values and agree with the measured volume susceptibility of the crushed sample.
Simulations are compared to experiment in Figs. 9 and 10 depicting magnetisation decay and T 2 distribution, respectively. We immediately note the good match between experiment and simulation for all dual-scale simulations. Grid refinement by a factor of 2 or 4 improves the match further, affecting mainly the peak position of the shorter T 2 peaks. This implies that the internal magnetic field gradient resolution of the transition region between coarse-and fine-scale structure again matters in agreement with the findings for the three-layer model. This is despite the relatively small clay fraction Consider now the spatial internal magnetic field gradient distributions across the pore space. Gradients in the solid are not considered and set zero for visualisation. We further differentiate between including ( G in ( )) and excluding ( G ex ( )) voxels at the solid-fluid boundary.

3
Cross sections of the spatial fields G in ( ) and G ex ( ) for a 15,000 2 sub-slice of the 80, 000 3 volume are given in Fig. 11 and the corresponding histograms for the full sample are depicted in Fig. 12.
The spatial visualisation of internal magnetic field gradients (Fig. 11) highlights the difference between including and excluding fluid-solid gradients in regions of high surface area, e.g. clay regions. For longer-range gradients in the macro-pore space caused by larger homogeneous regions, the fields G in ( ) and G ex ( ) are visually almost indistinguishable. The difference in the respective internal gradient distributions can be seen as sharper offset peaks in Fig. 12 (left). The position and magnitude of this offset changes with increasing grid refinement factor since the discontinuity between solid and fluid is better resolved for higher refinement, while the larger resultant gradients are distributed across fewer voxels. In the right column of Fig. 12, we compare the distribution of G ex to a normal distribution, leading to a good fit with fit parameters given in Table 3 and highlighting the limited shift in peak position.
Most of the values fall into the range of 10 −1 G/cm to 10 2 G/cm, the mean values of G ex are 1.09 G/cm (factor 1), 1.12 G/cm (factor 2), and 1.17 G/cm (factor 4), and maximum internal magnetic field gradients can reach 10 5 G/cm. We note that very high gradients in the fluid can be generated by the high-density paramagnetic inclusions (see red region in Fig. 11).
The histograms of magnetic field gradients explored by NMR simulations ( G ef f ) are shown in Fig. 13 for grid refinement factors of 1 and 4, where Eq. (13) is evaluated for N s = 1 2 t E = 100 μs. The shapes of the distributions are well described by normal distributions (see Table 4). The means of the diffusion-averaged gradients are around 1.07 G/cm. With increasing coarse voxel refinement, the G ef f distributions broaden. There is a relative rise of both tails in Fig. 13 compared to Fig. 12 and associated drop in peak value, which can be explained by diffusional averaging. The distributions are limited to the range of 10 −3 G/cm to 10 3 G/cm since actual histogram contributions outside this range  Fig. 12 Histograms of the internal magnetic field gradients of the water phase including ( G in ) and excluding ( G ex ) the water-fluid boundary for different coarse-grid refinement factors. a, b Factor = 1; c, d factor = 2; e, f factor = 4. Left Comparison of G in and G ex . Right Normal distribution fits to G ex (see also Table 3) Table 3 Normal distribution parameters corresponding to the fits of G ex in Fig. 12  are difficult to discern from zero. Random walks do experience no gradients when interacting with the surface, since near-surface effects are incorporated into the surface relaxivity. Similarly, the very high gradients possible in the fluid when close to a paramagnetic mineral affect a larger region when accounting for diffusional averaging. It is worth noting that the common approach to account for the effect of susceptibility contrast-induced internal magnetic field gradients is based on an analysis of relaxation regimes proposed by Hürlimann et al. (1995). It has been demonstrated that for the two asymptotic regimes ("free-diffusion" or "short time" typically applicable for larger pores; and "motional averaging" for smaller pores), the corresponding analytical expressions (Mitchell et al. 2010;Mitchell and Chandrasekera 2014) could be developed to account for the additional relaxation due to spin dephasing in magnetic field gradients. This enables standard petrophysical interpretations of the transverse relaxation responses in the fast diffusion limit, e.g. estimating pore size distributions, despite the presence of significant internal magnetic field inhomogeneities. For instance, it was applied to great benefit in the analysis of capillary trapping in sandstones subjected to partial desaturation by nitrogen and carbon dioxide (Connolly et al. 2017(Connolly et al. , 2019. However, the analytical approach requires the identification of a "mean" relaxation regime. Furthermore, the expressions so far were developed for two out of three relaxation regimes, being invalid for frequently reported intermediate regimes and coexisting regime mixtures. On the contrary, the dual-scale approach presented here does not require any assumption about governing relaxation regimes, offering potentially a tool to study associated problems. In particular, the coexistence of different regimes and sig-   Fig. 13 Histograms of the diffusion-averaged effective gradients G ef f according to Eq. (13) with N s = 1 2 t E = 100 μ s for different coarse-grid refinement factors. a Factor = 1, b factor = 4. The solid lines are Gaussian distributions over log(G ef f ) with respective parameters reported in Table 4   Table 4 Normal distribution parameters corresponding to the fits of G ef f in Fig. 13

Conclusions
This study introduced an effective NMR relaxation response modelling approach incorporating resolved internal magnetic fields enabling numerical simulations spanning five orders of magnitude in length scales on grids logically containing of the order of 10 15 cells. The internal magnetic fields were calculated via the dipole approximation in a dual-scale approach with trend model. The approach was validated against a model micro-structure and demonstrated for the Bentheimer sandstone as common benchmark rock. For the latter, NMR T 2 responses of a (1.76 mm) 3 domain were derived with grid resolution of 22 nm based on micro-CT images at 2.2 μm resolution and a clay microstructure resolved at 22 nm resolution. Comparison of T 2 distributions between simulation and experiment show excellent agreement. A significant result is the influence of tri-linear interpolation on the position of shorter relaxation time modes. This is associated with reduced errors in the transition regions between macro-and micro-scale regions, where near-field effects become important. We note that the change in peak position cannot be explained by a change in surface relaxivity of the kaolinite model in practice, since the underlying cause is the transition region. This has significant implications for modelling more complex rock samples where clay fractions are larger and paramagnetic impurities more frequent.
Finally, we comment on the observed fine-scale internal field gradient distributions. Spatially we observe that larger gradients are associated with surfaces in an irregular way-with higher gradients at the solid-fluid boundary, which depending on the pore shape may extend more or less into the fluid. The diffusion-averaged effective gradients are well described by log-normal distributions.