Accurate phase-shift velocimetry in rock

Spatially resolved Pulsed Field Gradient (PFG) velocimetry techniques can provide precious information concerning ﬂow through opaque systems, including rocks. This velocimetry data is used to enhance ﬂow models in a wide range of systems, from oil behaviour in reservoir rocks to contaminant transport in aquifers. Phase-shift velocimetry is the fastest way to produce velocity maps but critical issues have been reported when studying ﬂow through rocks and porous media, leading to inaccurate results. Combining PFG measurements for ﬂow through Bentheimer sandstone with simulations, we demonstrate that asymmetries in the molecular displacement distributions within each voxel are the main source of phase-shift velocimetry errors. We show that when ﬂow-related average molecular displacements are negligible compared to self-diffusion ones, symmetric displacement distributions can be obtained while phase measurement noise is minimised. We elaborate a complete method for the production of accurate phase-shift velocimetry maps in rocks and low porosity media and demonstrate its validity for a range of ﬂow rates. This development of accurate phase-shift velocimetry now enables more rapid and accurate velocity analysis, potentially helping to inform both industrial applications and theoretical models.


Introduction
Fluid flow through porous media, such as rock or sand packs, is found in a wide range of industrial and natural processes ranging from chemical reactors to petroleum recovery. Knowledge of the flow properties in these media can be crucial in understanding transport processes and developing accurate transport models. Nuclear magnetic resonance based approaches enable the complexity of local flow processes within the system to be characterized, moving our understanding of flow beyond bulk average macroscopic descriptions. NMR based approaches have been used to, for example, explore simultaneous flow of oil and water in sandstone [1], unpick complexities in nanoparticle transport behaviour in rock [2], map organic pollutant transport in fractures [3] and image heavy metal removal in bio-film mediated ion exchangers [4]. The Pulsed Field Gradient Nuclear Magnetic Resonance (PFG NMR) experiment originally proposed by Stejkal and Tanner [5], has long been used to non-invasively study flow and diffusion properties [6]. Furthermore, localised measurement of flow properties can be achieved by combining PFG with an imaging module to give PFG velocimetry also known as Magnetic Resonance Velocimetry (MRV). The resulting spatial maps of velocity provide a rich insight into the transport and structural properties of optically opaque systems.
There are two main methods of PFG velocimetry, namely propagator velocimetry and phase-shift velocimetry. Propagator velocimetry consists of resolving the probability distribution of displacements for each voxel. These are slow to acquire, requiring at least 8 [7] and up to 128 [8] gradient encoding steps (or q values). Phase-shift velocimetry is faster, requiring only two gradient encoding steps to measure the average velocity in each voxel. Indeed, phase-shift velocimetry is at least 4 times faster than propagator velocimetry, and is thus a highly desirable alternative when experiments can have time durations of days.
In the application of PFG velocimetry to porous media, it is useful to distinguish two regimes. In the first regime, where the imaging voxel size is smaller than the typical pore size, e.g. bead packs, both of the PFG methods are used and found to be reliable [9]. In the second regime, where the voxel size is greater than the typical pore size, e.g. sandstone rock [10], though there have been reports of quantitative phase-shift velocimetry [11], it has been generally advised to use the more time consuming propagator method [12], as numerous issues have been reported with the use of phase-shift velocimetry. These issues can be broadly categorised as: 1. Measured average velocity values not agreeing with values calculated from the known flow rate and porosity. Lower values than expected are reported at higher flow rates [13] making the relationship of measured velocity to the imposed flow rate non-linear [14,15].

Standard deviation of voxel velocities exceeds the expected val-
ues. This effect becomes stronger at lower flow rates [16] with a large proportion of the voxels unexpectedly indicating negative velocities [15]. 3. Measured velocity can vary with experimental PFG parameters.
Several authors have shown that at fixed flow rate, different velocity values are measured when different gradient strengths (G) [12] or observation times (D) [15] are used.
These issues have effectively made phase-shift velocimetry unreliable for use with porous media like rocks, where voxel sizes can be greater than the typical pore size. In this work, we clearly characterise the above mentioned problems, identify their underlying causes and propose concrete solutions for producing accurate phase-shift velocimetry measurements in rocks and porous media.

PFG NMR velocimetry
PFG NMR velocimetry consists of making the phase of the NMR signal sensitive to translational motion. This is achieved by applying a pulsed field gradient of amplitude G during a time d, imposing spatially dependent phase shifts to the spins. For a spin moving along the path rðtÞ, the induced phase is given by: After an observation time D, rephasing gradients are applied to the system. By choosing parameters such as d ( D (narrow pulse approximation), one can neglect displacements that occurred during d. Then, for a spin starting at r 0 and ending at r 0 þ R, the resulting phase-shift is given by cdG Á R. At this stage, to describe phase modulation due to molecular motion, it is often convenient to introduce the wave vector q ¼ cdG. The wave vector q is the conjugate of spin displacement in the same way that the wave vector k ¼ R t 0 gðtÞdt is the conjugate of spin position in an imaging experiment [17]. The combination of velocity encoding and imaging allows to measure phase-shift for each voxel in the sample.
The NMR signal resulting from a spatially resolved PFG NMR experiment can be expressed by: Moran [18] showed that, for a spin at position r with a displacement R during the time D; spin density qðrÞ could be generalised to a joint density function, q D ðr; RÞ; defined as: where P D ðR; rÞ is the normalised probability distribution function of spin displacements over the period D, also called a propagator. By applying the velocity encoding gradients along a single direction (for example z) and considering a displacement Z of each spin during the time D, the combination of Eqs. (2) and (3) gives the NMR signal for a voxel situated at position r as: Defining the average velocity of each spin during D as v ¼ Z D , it is possible to rewrite Eq. (4) as: If the time integral of the velocity encoding gradient is zero, this integral is independent of spin position and Sðr; qÞ is the Fourier transform of the velocity-density function P D ð v; rÞ.

Propagator velocimetry
One approach to measure velocity, called propagator velocimetry, consists in acquiring Sðr; qÞ for several q values, or q-steps, and then apply an inverse Fourier transform in order to obtain the propagator P D ð v; rÞ. The number of q-steps and their size has to be selected appropriately so as to cover the velocity range found in each voxel and get the desired propagator resolution. Typically a minimum of q steps has to be used, which leads to significant experimental times, even when using fast acquisition sequences [19].

Phase-shift velocimetry
In another approach, velocity is related to the phase of the signal resulting from a PFG measurement. First, by inserting the expression of the ensemble averaged velocity for a voxel, If the velocity density function is symmetric around the mean velocity v then the integral in Eq. (6) is real and the phase of the resulting signal is found to be proportional to the average velocity, VðrÞ, and the resulting phase is given by: In theory, by subtracting two phase images taken at equal D times, and with equal but opposite G values one can obtain a map with intensities proportional to velocity [20]. The second G value phase image cancels eddy current related phase contributions that are independent of q. This UðrÞ ¼ u 2 ðrÞ À u 1 ðrÞ map is easily transformed into a velocity map using Eq. (8): In practice, the measured phase-shift UðrÞ is affected by additional experimental parameters. The phase-shift effectively measured at any voxel r of a phased image can be expressed as [21]: where VðrÞ is the average velocity of spins, c is the gyromagnetic ratio, aðrÞ corresponds to phase contributions that depend of q and hðrÞ is phase shift caused by noise.
By acquiring a phase-shift map at zero flow, U 0 ðrÞ, it is then possible to remove phase contributions that are not flow related. The resulting phase-shift can then be rewritten as: The noise related phase-shift error is related the uncertainties in the measurement of the x and y components of the nuclear magnetisation in the rotating frame [22]. For an uncertainty DS in each direction, phase error can be estimated by [12]: h is therefore reciprocal of the signal-to-noise ratio (SNR). This makes SNR an important parameter to consider since hðrÞ ( cdGDVðrÞ is a condition for producing accurate velocity maps. Phase-shift velocimetry relies on the linear relation between the phase of the NMR signal and the imposed velocity-encoding gradient, which enables to use Eq. (10) for the production of a velocity map. But it has been shown in rocks that this phasegradient linearity gets compromised as gradient increases [12], with the linear range becoming smaller at higher flow rates. Working at lower flow rates is not a solution, since the imparted phase shift is reduced, making the measurements prone to phase errors that can result in noise-dominated spatial velocity distributions.
Finally, it is important to stress the fact that the application of PFG velocity mapping techniques to rocks is also limited by magnetic susceptibility effects and short spin-spin relaxation times. These are caused by magnetic susceptibility differences between the fluid and the solid phases of porous materials resulting in internal field gradients. The effect of such internal gradients can be minimised by using alternating pulsed field gradients [23] and short echo time sequences [11].

Experimental
The MRI experiments were performed on a horizontal 7 T Bruker Avance Biospec system (300 MHz). A Bruker BGA12SL micro imaging gradient insert (400 mT m À1 ) and 200-A gradient amplifiers were used to provide linear magnetic field gradient pulses. The birdcage Radio-Frequency (RF) volume resonator used for all experiments had an inner diameter of 72 mm.

Experimental setup
All experiments presented in this work were performed on a Bentheimer sandstone sample with a diameter of 3.8 cm and a length of 7.6 cm. The rock core was fitted with inlet and outlet end caps and encapsulated with epoxy resin to make a watertight coating (Fig. 1). It was then vacuum-saturated with deionised (DI) water and flushed with 20 pore volume of DI water to remove particulates and readily soluble salts. Prior to the MRI experiments the rock core was saturated by flowing water for 3 days (saturation was confirmed by T 1 relaxation measurements). The sandstone rock was finally placed in a leak preventing plastic cylinder and positioned inside the RF coil at the centre of the MRI bore. Air bleed outlets located at the end-caps were used to pump out any accumulated air bubbles. Flow was controlled using HPLC isocratic pump (Agilent 1100 series). The flow rate, Q, was varied from 0.5 to 4 ml min À1 .
The porosity of the core was measured by weighing it in both dry and saturated states. The difference of these masses corresponds to the mass of fluid contained in the core. The absolute porosity was then calculated from the ratio of the water volume to the total core volume. This measurement gave a porosity of 17 ± 0.5%.

MR techniques
Relaxation measurements were performed prior to velocimetry ones in order to measure water relaxation times in the rock. The spin -lattice relaxation, T 1 , was measured using an Inversion Recovery sequence and the spin-spin relaxation, T 2 , using a Multi Slice Multi Echo sequence.
Velocity measurements were performed using a combination of the Alternating Pulsed Gradient Stimulated Echo (APGSTE) pulse sequence with a Rapid Acquisition Relaxation Enhancement (RARE) imaging module. The pulse sequence was implemented in-house and calibrated by measuring the velocity of water flowing through an unobstructed tube [15]. The APGSTE pulse sequence has been shown to cancel the cross term between applied gradient and background gradients, related to magnetic susceptibility effects. The stimulated echo (STE) approach ensures that the displacement of spins occurs during a z-storage interval, reducing the signal loss related to fast T 2 decays induced by the internal gradients. Although recent results [24] suggest its limited accuracy when very high flow velocities are encountered, the APGSTE sequence is well suited for measurements of slow translational motion processes.
The voxel size for all experiments was 1 mm 3 . For 2D velocity maps the field of view was of 60 Â 44 mm 2 while for 3D velocity maps it was of 60 Â 44 Â 44 mm 3 . For all velocity measurements the duration of the 90°pulse was 1.2 ms, the duration of the 180°pulse was 2.4 ms, the echo time (TE) was 5.7 ms and the repetition time (TR) was 5000 ms. A RARE factor was of 2 was used for all the experiments.
2D velocimetry experiments were performed on a 1 mm slice along the length of the rock. The duration of the flow encoding alternating gradients was of 1 ms (d = 2 ms), the observation time, D, varied from 25 ms to 200 ms and the gradient strength from À25 mT m À1 to 25 mT m À1 . Phase-shift velocity maps were produced using phase data from 2 q-space points while for the propagator measurements 32 q-space points were acquired. For 2D phase-shift velocimetry maps, 2 averages were used and the total experimental time for a q-space point was of 3 min and 40 s. For the 3D phase-shift velocimetry maps, 6 averages were used and the total experimental time for a q-space point was of 8 h and 4 min. The experimental results presented in this work were obtained from the same region of interest (ROI) of 10 mm in the centre of the rock for 2D maps ($400 voxels) (cf. Fig. 1) and of 40 mm in the centre of the rock for the 3D maps ($50,000 voxels).

Results and discussion
Relaxation measurements performed prior to velocity mapping gave average values within the rock of T 1 = 1500 ms and T 2 = 34 ms. This long water T 2 and the fact that the decay curves were all strongly single exponential, indicates low magnetic susceptibility effects in the system and low levels of paramagnetic impurities in the sandstone. It is worth noting that the Bentheimer sandstone used in this work is an extremely clean outcrop; far from typical of most rocks types, where extremely short relaxation times, particularly at high magnetic field, can prevent the use of conventional MRI pulse sequences. Though, the use of short echo time sequences has been demonstrated to overcome this restriction [11]. Average signal to noise ratio (SNR) was calculated as the ratio of the average signal in the rock region on the average signal in the background. For all experiments presented in this work SNR was superior to 200.

SNR related phase measurement noise
Measurements made at zero flow, with the same parameters as flow experiments, are commonly used for correcting phase measurements by eliminating phase contributions aðrÞ that are depending on the q-value (Eq. (10)). They can also inform about the noise level in the phase measurements.   2a shows the average phase in the ROI against velocity encoding gradients. Using these results and Eq. (10), it is possible to calculate the corresponding velocities at Q = 0 ml min À1 (Fig. 2b). Motion measured at stationary flow is diffusive only, hence average velocity values in the ROI are very low (between À0.3 and À5 lm s À1 ). Fig. 2c shows the standard deviation of voxel velocity values for the same ROI. Note that despite the fact that fairly accurate velocity measurements are produced for all the range of gradient values used here, at low gradient measurements the noise in the phase measurement, h, dominates the resulting velocity standard deviation within the ROI. As a consequence, the resulting maps might produce accurate average velocities if a sufficient number of voxels is considered ( R hðrÞdr ¼ 0), but the individual voxel values can be inaccurate, with measurements errors up to ten times the average value. For example, the velocities obtained for a gradient step of 0.015 T m À1 and D = 50 ms have an average of À5 lm s À1 with a standard deviation of 70 lm s À1 .
As the encoding gradient increases velocity standard deviation is reduced. This is due to the fact that the total imparted phase increases and phase noise h becomes a negligible component. For velocity encoding gradients above 0.05 T m À1 this standard deviation becomes constant at a value that represents the physical standard deviation of velocities within the sample. As the observation time D increases the physical standard deviation is also shown to decrease. This can be related to the fact that the moving protons experience larger trajectories allowing their average velocity to shift towards a common average depending on the local average diffusion. In fact the effect of the observation time is not so straight forward as increasing D increases the imparted phase which reduces noise influence but also reduces SNR which increases the noise effect.
These phase noise effects have to be considered carefully when performing spatially resolved velocimetry. Subtracting the phase from the zero flow experiments will correct from parameter depending phase contributions aðrÞ but might introduce phase noise h 0 ðrÞ that can compromise the accuracy of the resulting velocity maps (cf. Eq. (10)).

Non-linear phase-gradient relation
Once parasitic phase contributions are removed, one is left with Eq. (10) relating measured phase-shift to velocity. The linear relationship between phase and gradient is a condition for this equation to produce accurate velocity calculations. For low gradient values (below 0.03 T m À1 ) the phase varies linearly with the gradient. As gradient increases this linearity gets compromised. This effect is shown to be stronger at higher observation times. At D = 50 ms a linear fit (R 2 > 0.998) reproduces well the behaviour of the curve, but at higher observation times the  data becomes increasingly non-linear. Comparison between experiments performed at Q = 1 ml min À1 (Fig. 3a) and Q = 2 ml min À1 (Fig. 3b) suggest that flow rate increase also compromises the phase-gradient linearity. This same effect was also identified by Chang et al. [12]. The similarity between curves obtained at the same Q Â D (e.g. Q = 1 ml min À1 /D = 200 ms and Q = 2 ml min À1 / D = 100 ms) is also noticeable. The product Q Â D is proportional to the average molecular displacement Z avg ¼ V Â D in the course of the experiment. Fig. 4 shows the same phase data as Fig. 3b plotted against q Â D (with q ¼ cdG). The superposition of these curves suggests that in this system, for a given gradient value the observed phase-gradient behaviour is strongly depending on Q Â D.
Quantitative velocity measurements are possible with the range of parameters that produce a linear phase-gradient relationship. The symmetry of the graph with respect to the ordinates axis underlines the importance of using opposite sign gradients for performing phase velocimetry experiments. The linearity of phase against gradient for low gradient values was exploited by Romanenko et al. [11,25] who performed fast phase-shift velocimetry in rocks by considering only a limited number of q-space values near the origin, relating velocity to the slope of the curve.

Simulation of phase-gradient relation from propagator data
Before phase-shift velocimetry can be used reliably in porous media such as rocks, it is crucial to identify the cause of the nonlinearity of the phase against gradient. The literature contains several proposed explanations, ranging from flow related eddy currents [14], velocity distribution asymmetry within the voxels [12], relaxation effects [15] and other non-identified effects [14].
To answer this question, we started by examining the effect of propagator asymmetry for flow through the rock core. From 2-D images acquired using the APGSTE sequence with 32 equally spaced gradient values (ie positive and negative q values), complex signal was taken from a single voxel in the centre of the rock. From this data, we were able to produce the displacement probability distributions (i.e. displacement propagators) by inverse Fourier transformation of the signal from the 32 q values. Fig. 5a shows normalised propagators in a single voxel in the centre of the rock, measured for flow rates of 1 and 2 ml min À1 with an observation time of 100 ms. As expected, the propagator obtained for Q = 2 ml min À1 is less symmetric than the propagator for Q = 1 ml min À1 .
Using these propagator data it is possible to interrogate phase behaviour by calculating phase in three different ways:

Measured phase
From the complex signal real ðS X Þ and imaginary ðS Y Þcomponents, the experimental phase-shift for each gradient value can be calculated as,

Average phase
This corresponds to the average phase that is actually imparted on the individual spins and it can be simulated for each gradient by using the displacement probability distributions. Each point, i, in the propagator relates a displacement Z i to its probability P i (Fig. 5a). For a given gradient value G the phase imparted by a displacement Z i is u i = cdGZ i (Fig. 5b). The total phase imparted in a voxel presenting the distribution of displacements given by an npoint propagator is therefore:

Simulated measured phase
It is very important to stress that the PFG sequence do not measure the above-mentioned ''average phase" actually imparted on the individual spins but the phase generated from the sum of the real components and imaginary components from all spins. From the imparted phase u i obtained from average phase simulations using the propagator data, it is then possible for each gradient value to calculate the real S Xsim and imaginary components S Ysim using: The resulting phase that will be measured by this simulated experiment is then calculated by adapting Eq. (12),  diagram: when a gradient G is applied, the phase u i(G) will be imparted by a displacement Z i and the corresponding real (S Xi ) and imaginary (S Yi ) components will be measured by the coils.
A MATLAB code was developed to allow the above-mentioned phase calculations to be carried out using experimental PFG data. Fig. 6 shows plots of the experimentally measured phase, the simulated average phase and the simulated measured phase produced using the propagator data shown in Fig. 5a. As expected, the average phase (black line) varies linearly with the gradient strength. In the measured phase-shift obtained both by experiment and simulation the phase-gradient linearity gets compromised above a certain gradient strength that depends on the flow rate. The agreement between experimental and simulated results is excellent for both flow rates presented here. This strongly suggests that propagator asymmetry alone can fully account for the nonlinearity seen in phase-gradient relation.
Furthermore, the equations used for the simulations of the measured phase-shift can inform us on the nature of the non-linearity.
For example, it can be noted that when the angle u is small (i.e. cos (u) ? 1, sin(u) ? u and arctan (u) ? u) Eq. (16) becomes: Hence one expect accurate u average measurements to be performed at small phase angles u. This knowledge allows conclusions to be drawn regarding parameter regions that will provide a linear phase-gradient relation allowing accurate velocity measurements to be performed. To this respect, relative propagator symmetry appears as a necessary condition.

Propagator asymmetry effect on phase-gradient relation
We have demonstrated that velocity distribution asymmetries are the source of non-linearities in the phase-gradient relation and hence of measurement errors encountered in the literature. It is therefore important to discuss the source of these propagator asymmetries. The evolution of the displacement propagator with observation time can be separated into three distinct stages, summarised in Fig. 7.

Diffusion dominated region (symmetric propagator)
In this region the average displacement due to flow (Z avg ¼ V Â D) is negligible compared to those from self-diffusion ). As the diffusive process is symmetric, the measured propagators tend to symmetric.

Intermediate displacement region (asymmetric propagator)
In porous media, especially natural media like rock, there are stagnant zones or pores which do not flow. For water molecules in these stagnant zones to experience flow, the molecules need to diffuse into flowing zones. This process takes time and leads to a prolonged peak in the propagator at zero displacement. When combined with molecules that are flowing, this results in an asymmetric propagator.

Long displacement region (symmetric propagator)
At longer observation times most of the stagnant water molecules can all diffuse into the flowing zones. This combined with mixing due to dispersive processes (e.g. mechanical, Taylor) results in a symmetric displacement propagator at sufficiently high flow rates and long observation times.
As described above, for accurate phase-shift velocimetry it is crucial to avoid intra-voxel propagator asymmetry. Some researchers have suggested working in the long displacement region, by using long observation times [12]. However, this comes at a serious loss of signal from relaxation, which in turn can lead to the introduction of phase noise as described in Section 3.1. Also, in some cases, the displacements required to access this region may exceed those that can be measured by NMR. Alternatively, we propose to work in the diffusion dominated region, which has the twin benefits for accurate velocimetry of both high SNR and symmetric propagators. The main limitation of this approach would be the requirement for short d and D.

Accurate phase-shift velocimetry
Velocity maps for Q = 1 ml min À1 and Q = 2 ml min À1 were produced using the phase data presented in Fig. 3 and Eq. (10). Zero flow phase was used for correcting the baseline. Velocity measurements perpendicular to the flow direction were also performed Fig. 6. Average and measured phase-shift obtained by simulations using the propagator data and experimental phase-shift extracted from the same data. The flow rate was (a) 1 ml min À1 and (b) 2 ml min À1 . Fig. 7. Schematic diagram showing the propagator shape as a function of the average molecular displacements Z avg during D. (Fig. 8a), confirming that any transverse components were orders of magnitude smaller. Fig. 8 shows the average velocity in the ROI against the velocity encoding gradient step at different observation times. The average velocity is almost constant for the lower gradient steps. Note that the regions of constant velocity measurements correspond to the linearity region of the graph of phase against gradient. This explains the fact that these regions are reduced as flow rate and observation time increase. The measurements produced outside the linearity region of the graphs of phase against gradient show a decay of velocity with increasing gradient. This decay is stronger for higher observation times and flow rates, in total agreement with issues previously reported in the literature [12,26].
In the regions where a linear fit gives R 2 > 0.99, the average velocities are of 0.106 ± 0.006 mm s À1 for Q = 1 ml min À1 and of 0.209 ± 0.008 mm s À1 for Q = 2 ml min À1 . These results are in excellent agreement with propagator measurements of velocity obtained by dividing the average displacement of the propagators presented in Fig. 5 by the observation time (less than 4% difference).
At this point, one might want to consider the average displacement for which velocity starts decaying. For Q = 2 ml min À1 and D = 100 ms (Fig. 8b), this seems to happen for a gradient step of 6 mT m À1 . With an average velocity of 0.209 mm s À1 , the average displacement during D is approximately 21 lm. With the pore size distribution in the Bentheimer sandstone ranging from 10 lm to 100 lm [27] one expects to find asymmetric displacement distributions if diffusion is not dominating the propagators. The average displacement caused by diffusion during D is given by . For a water self-diffusion of 2 Â 10 À9 m 2 s À1 this corresponds to a Z diff of 20 lm. This is in agreement with our hypothesis that the accuracy of velocity measurements decreases when Z diff ) Z avg is not valid and the rock structure starts producing non symmetric displacement distributions. We have shown that accurate average velocity measurements can be performed in rock but for accurate velocity maps care needs to be taken with phase noise.
3.6. Velocity maps with negligible phase noise Fig. 9 shows the voxel velocity standard deviation in the ROI against the velocity encoding gradient step DG at different observation times for velocity maps at Q = 1 ml min À1 and Q = 2 ml min À1 . The data used here are the same as for Fig. 8. As for the zero flow measurements, at low gradient measurements (below 2 mT m À1 ) the noise in the phase measurement, h, dominates the resulting velocity standard deviation within the image. Hence, for measurements at low gradient values, phase noise has to be considered carefully since when enough voxels are considered, the average velocity of the ROI might be accurate, but individual voxel velocities maybe highly inaccurate due to phase noise  errors. Such noise errors might constitute a limitation to phaseshift velocimetry techniques based on measurements close to the q-space origin, as the one implemented by Romanenko et al. [25].
It becomes clear that in order to produce accurate velocity maps, one must select low enough gradient values so as to produce accurate velocity measurements (Fig. 8), and high enough gradient value so as to impart enough phase, making phase-noise negligible. For example, for a flow rate of 1 ml min À1 and an observation time of 100 ms, a gradient strength of 3 mT m À1 satisfies both conditions. Fig. 10(a and b) shows axial and radial direction slices from an accurate 3D velocity map produced using this set of parameters. A 3D velocity map was also produced for Q = 2 ml min À1 with an observation time of 50 ms. The velocity distribution within voxels of the 3D velocity maps obtained for Q = 1 ml min À1 and Q = 2 ml min À1 can be seen in Fig. 10c. As expected, average velocity Q = 2 ml min À1 (0.188 mm s À1 ) is the double of the average velocity for Q = 1 ml min À1 (0.0943 mm s À1 ). Fig. 10d shows the average velocity in each one of the planes along the length of the rock. For both flow rates the average velocity measured is constant along the length of the rock indicating that conservation of mass is met in the velocity map. The standard deviation of velocities for Q = 2 ml min À1 (0.0747 mm s À1 ) is the double of the one for Q = 1 ml min À1 (0.0345 mm s À1 ). This indicates that the phase noise component, introducing an additional and constant velocity spread around the mean value, is negligible. Contrary to what has been seen in previous works [15] for comparable flow rates, the proportion of negative velocities is extremely small and concerns voxels situated at the edges of the rock, where partial voxel filling might cause phase measurement errors.
To further validate these measurements, the velocity results were compared with theoretical calculations of the average velocity based on porosity (£) measurements. An estimation of the average velocity in the whole rock can be obtained using the equation V ¼ Q £A where A is the surface of the rock cross section. Average velocities obtained were 0.087 mm s À1 for Q = 1 ml min À1 and 0.173 mm s À1 for Q = 2 ml min À1 . These values are in good agreement with the measured values (less than 10% difference).

Conditions for accurate spatially-resolved phase-shift velocimetry
Sections 3.5 and 3.6 that led to the production of accurate spatially resolved velocity maps allow to draw some general conclusions so as to the conditions for accurate spatially-resolved phase-shift velocimetry:

Symmetric displacement distribution within each voxel
We showed that this could be achieved by ensuring that displacements are either diffusion dominated (small displacements) either dispersion dominated (long displacements) over the observation time of the experiment. For molecular displacements to be diffusion dominated so that the displacement propagator can approach symmetry, the diffusive component has to be greater than the advective ( ffiffiffiffiffiffiffiffiffiffi 2DD p > VD). The value of the observation time needs to respect the following condition: In general, and as seen in Section 3.5, by keeping the mean phase-shift small (small D and GÞ one obtains more symmetric displacement distributions, increasing the accuracy of the measurement. In the case of relatively homogenous media, and for a particular interstitial velocity V, one can relate the observation time directly to experimental parameters as the imposed flow rate, Q . Using Eq. (17) and the relation between the flow rate and average velocity, Q ¼ V£A, where A, is the rock cross section and £ the absolute porosity, one obtains:

Negligible phase noise
A sufficient phase-shift must be imparted so as to neglect phase-noise. This can be expressed as: As seen in Section 3.6, by keeping the mean phase-shift big (big D and GÞ one reduces the phase-noise of the measurement.
It can be clearly seen that conflicts might arise when one tries to fulfil both conditions, and a good compromise between accuracy and dynamic range must be found so as to avoid compromising one of the two. In our case, and for a SNR of 200 achieved with 4 averages, phase-shifts below 0.35 rad were ensuring negligible deviations from linearity for the phase-gradient relation in the range of studied observation times. Also, for a phase-shift of 0.35 phase-noise was below 3% of the imparted phase, ensuring negligible phase-noise effect in the velocity maps.
To further validate the above-mentioned conditions we decided to consider a range of flow rates. In order to ensure that the propagators will stay symmetric (diffusion dominated) and phase noise will not affect the measurements one can use the average displacement analysis presented earlier. In fact, by keeping the average displacement Z avg small compared to the average diffusion displacement Z diff one should obtain the same phase-gradient relation and high enough imparted phase so as to neglect phase noise. The average displacement being proportional to the product Q Â D, it can be controlled for a given flow rate Q by adjusting the observation time D. Fig. 11 shows the resulting phase for different flow rates, measured at constant Z avg equal to the one that allowed producing the accurate velocity maps presented in Fig. 10. As expected the linear regions of the phase to gradient plot are perfectly superposed. Fig. 12 shows the average velocity and velocity standard deviation in the ROI against the velocity encoding gradient step, DG, obtained using the phase data presented in Fig. 11. Small velocity differences for the first 4 points of the plots in Fig. 12a are negligible (less than 4% of the measured velocity) and could either be caused by the phase noise shown in Fig. 12b or by the zero flow phase noise that was shown in Fig. 2a.
The region of constant velocity at low gradient and the region of negligible phase noise are not affected by the flow rate change; it is therefore possible to produce accurate velocity maps at the same gradient value used to produce the velocity maps presented in Fig. 10 (gradient step of 4 mT m À1 ). Fig. 13a shows a plot of the average velocity in the considered ROI against flow rate. It shows the expected linear relationship (R 2 > 0.998) for all the range of studied flow rates and passes from the origin. The accuracy of the average velocity measurement being verified at Q = 1 ml min À1 using propagator measurements (cf. Section 3.5) and the linearity of the resulting plot indicate that accurate velocimetry was achieved for all the range of flow rates. Fig. 13b shows the standard deviation of voxel velocities in the produced maps. It also shows a linear relationship (R 2 > 0.999) for all the range of studied flow rates and passes from the origin. The constant component that would be introduced by phase noise is clearly negligible in these maps. This is also easy to observe in velocity profiles taken from the centre of the velocity map (Fig. 13c). The velocity variations in these profiles are maintained as flow rate increases. These variations are therefore caused by the rock local porosity variations and do not exhibit any noise. To the best of our knowledge these are the first phase-shift velocimetry experiments in rock that are proven to satisfy both velocimetry accuracy and map noise elimination for a range of flow rates.
Depending on the instrument and the PFG sequence used, optimal gradient values might not be easy to obtain at higher flow rates due to the fact that there is a limit in the reduction of the observation time. Considering the set up presented in this work, by reducing the gradient duration, d, to 0.5 ms and D to 5 ms it is possible to achieve Z avg ( Z diff , and produce accurate velocity maps, for flow rates up to 32 ml min À1 . For even higher flow rates the approach could be changed, and much longer Z avg considered. In that case, dispersion will become dominant allowing symmetric velocity distributions to be measured, thus preserving the linearity of phase versus gradient.

Conclusion
Studying water flow through Bentheimer sandstone, we have proposed a complete method for the production of accurate phase-shift velocimetry maps in rocks and other porous media.
Phase-gradient linearity is an essential condition for accurate phase-shift velocimetry. Simulations of measured phase revealed that the phase-gradient relation is wholly dependent on the distribution of intra-voxel displacements. Measured phase was shown to be sensitive to the propagator symmetry, with asymmetric propagators increasingly compromising the phase-gradient linearity as flow rate, gradient strength or observation time are increased. This has allowed us to explain the discrepancies in phase-shift velocimetry measurements that have been identified by previous researchers.
We showed, that accurate velocity measurements can be achieved by ensuring that the average displacement of molecules Z avg during the observation time, is smaller than the average diffusive displacement Z diff , thus producing a symmetric diffusiondominated displacement distribution. In addition, we have stressed the importance of considering phase noise for producing accurate velocity maps, showing how to avoid parameter regions were phase noise becomes dominant, introducing a noise component in the standard deviation of voxel velocities in the map. Finally, we proposed a general approach for identifying experimental parameters that allow accurate spatially-resolved velocimetry with negligible phase noise, and demonstrated its accuracy for a range of flow rates.
It is important to stress that the validity of our results regarding accurate PFG velocimetry is general and would apply to any other porous medium and experimental situation (e.g. experiments at lower magnetic fields or experiments using ultra-short imaging models). As phase-shift velocimetry is at least four times faster than propagator velocimetry, we are confident that this work will contribute to increasing use of phase-shift velocimetry in porous media research, helping to inform both industrial applications and theoretical models.