Magnetic Resonance Velocimetry for porous media: sources and reduction of measurement errors for improved accuracy

Nuclear Magnetic Resonance (NMR) is a powerful tool for noninvasively investigating static and dynamic properties of fluids inside complex and opaque structures. Thus, Magnetic Resonance Velocimetry (MRV) allows for the in situ analysis of the local flow velocity of fluids. Such an analysis characterises mass transport properties and helps to validate or improve numerical predictions for porous media. The benefit of such validations is lowered by several problems worsening the measurement accuracy. This publication addresses systematic errors and the influence of noise, which may reduce the accuracy of MRV measurements. Techniques are described to minimise or correct displacement and phase errors. A new multi-echo MRV-sequence is proposed as a compromise between Velocity-to-Noise Ratio (VNR), total measurement time, spatial resolution and displacement errors. For improved VNR, a dual-VENC encoding scheme was used and repetition time TR\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{\text{R}}$$\end{document} was optimised. As an example, the proposed MRV sequence was used to measure the three-directional flow velocity of water in an Open Cell Foam (OCF) with 10 pores per inch with three-dimensional isotropic spatial resolution.


Introduction
Understanding transport phenomena within monolithic catalyst carriers requires numerical and experimental tools. Morphological properties of monolithic catalysts such as Open Cell Foam (OCF) structures control the transport phenomena and have a considerable influence on the performance of the reactive system (Sinn 2021). Generally, OCFs are classified based on the number of pores per inch. The structure of OCFs is determined by their pore size, window size, porosity and specific surface area. For measuring flow velocity different methods are used like Laser Doppler Velocimetry, Particle Image Velocimetry or Planar Laser Induced Fluorescence (Northrup 1991;Saleh 1992;Lehwald 2008;Mihailova 2015), which need optical access to the fluid. A non-optical technique is Positron Emission Particle Tracking, where a tracer particle is added to the fluid (Langford 2016). The benefit of noninvasive and non-optical Phase Contrast (PC) Magnetic Resonance Velocimetry (MRV) is that, without adding any tracer particles, the three-directional velocity of a fluid in opaque structures can be measured with three-dimensional spatial resolution. Therefore, PC-MRV is widely used in medical application, for example to measure cerebral (Dumoulin 1989;Marks 1992) and aortic (O'Donnell 1985;Nayler 1986;Dulce 1992;Kvitting 2004;Frydrychowicz 2011) blood flow and cerebrospinal fluid pulsation (Enzmann 1991(Enzmann , 1993Linninger 2007). PC-MRV is also an important tool to analyse flow velocity in engineering applications like porous media (Sederman 1997;Ren 2005;Onstad 2011;Huang 2017;Penn 2018;Sadeghi 2020;Clarke 2021). Velocity fields are measured by encoding the signal phase of moving spins (Moran 1985;O'Donnell 1985). For measuring one-directional velocity typically two steps are done with velocity encoding gradients applied in opposite direction. This is done to subtract phase errors originating, e.g. from B 0 inhomogeneities, concomitant gradients or eddy currents (Bernstein 2004). Threedirectional velocity encoding is realised by using at least four steps with flow encoding gradients applied in different spatial directions. A simple four-step encoding scheme uses a single reference scan with velocity encoding gradients set to zero and three scans with velocity encoding gradients applied in x-, y-or z-direction (Hausmann 1991). Alternatively, a balanced four-step encoding scheme can be used where the direction of flow encoding gradients is altered in pairs for each step (Dumoulin 1991). By linear combination of the four datasets the velocity components in x-, y-and z-direction are reconstructed. A comparison of the different flow encoding schemes can be found in ; Conturo (1992); Dumoulin (1995). In  it was shown that simple and balanced four-step encoding have the same average noise performance. However, in contrast to the simple scheme, the noise in the measured velocity components is uncorrelated and velocity encoding for the same dynamic range is achieved with lower gradient strength.
One challenge is the optimal choice of the Velocity ENCoding (VENC), which determines the velocity range that can be uniquely assigned. In systems with a wide range of velocities, single-VENC measurements result either in aliasing for fast flow or a low Velocity-to-Noise Ratio (VNR) for slow flow. In the past years several aliasing-correction algorithms were presented (Herráez 2002;Jenkinson 2003;Schofield 2003;Abdul-Rahman 2007. However, such numerical algorithms may fail for many wrappings or for small pores, which is a problem for analysing porous media. Unwrapping techniques can also be realised by measuring additional encoding steps. In Lee (1995) the two-step encoding (one-directional velocity) was expanded by adding an extra step such that a higher VENC-value was achieved. This step was used to unwrap the aliased low-VENC data and therefore improvement of the VNR was achieved. A similar dual-VENC technique was used for three-directional velocity (Nett 2012;Schnell 2017;Carrillo 2019;Callahan 2020).
In clinical application it is of central interest to keep the measurement time short to reduce costs, avoid patient discomfort and minimise blurring of the images caused by movements (e.g. breathing). Due to their shorter measurement time, gradient-echo sequences with small flip angles and short repetition time T R are the common choice for in vivo measurements. Also they enable a low specific absorption rate (heating of tissue) (Foster 2007). However, a disadvantage of this technique is the dependence on the effective transverse relaxation time T * 2 and thus the susceptibility to magnetic field inhomogeneities. The advantage of Spin Echo (SE) sequences is their independence of T * 2 and therefore better image quality is achieved in case of magnetic field inhomogeneities. Thus, SE sequences are a good choice for measuring flow velocity in porous media (Sederman 1997;Ren 2005;Huang 2017;Sadeghi 2020;Williamson 2020;Clarke 2021).
Eddy currents depend on the strength and direction of the applied gradient. Thus, phase errors originating from the velocity encoding gradients do not subtract completely in the phase difference image and cause a spatial dependent offset (Bernstein 2004). To correct such offsets, an additional measurement with same parameters but without flow can be done (Sederman 2004;Sankey 2009;Shiko 2012;John 2020;Clarke 2021). The offset can be corrected by directly subtracting the phase values of the no-flow dataset voxelwise. Doing so, the VNR of the velocity map is reduced by √ 2 . Alternatively the no-flow dataset can be fitted by a polynomial lowering the VNR-loss (Lingamneni 1995;Busch 2017;Wüstenhagen 2021). However, it must be considered that for an inaccurate fitting systematic errors may occur lowering the benefit of this technique. To our knowledge this technique was not shown before for MRV measurements on porous media and also not for dual-VENC data.
Flow through porous media like OCFs, which have different sized pores, usually has a large difference between slow and fast flow regions. For both, slow and fast velocities, there are characteristic problems which lower the accuracy of measured velocity fields. This is especially disadvantageous if MRV measurements are used for cross-validation with Computational Fluid Dynamics (CFD) simulations. Slow flow regions, which are especially located at the pore surface, are often negatively influenced by partial volume effects (Tang 1993) while displacement errors have a negative effect on fast flow regions. Displacement errors may occur because the fluid moves between encoding events, in particular if the delay between spatial encoding and velocity encoding is long. Thus, velocity values would be assigned to wrong positions (Larson 1990;Nishimura 1991;Kouwenhoven 1995;John 2020). To reduce this effect, the voxel size of the measurement and timings of the sequence should be adjusted to occurring velocities. In Bruschewski (2019) a single-point PC sequence was proposed to minimise displacement errors. However, the point-wise acquisition leads to long acquisition times and is therefore not suitable for 3D imaging.
For reduced measurement time in Sederman (2004) an echo-planar imaging sequence with velocity encoding was proposed. In Shiko (2012) and Huang (2017) Rapid Acquisition with Relaxation Enhancement (RARE)-PCMRV-sequences were presented for reduced measurement time. After one excitation pulse several echoes are acquired. Because of displacement errors such sequences are limited to slow flow velocity. To reduce the displacement effect, the FLow Imaging Employing Single-Shot ENcoding (FLIESSEN) sequence was proposed, which is a single-shot multi-echo sequence with velocity encoding and decoding for each echo (Amar 2010).
The objective of the present work was to describe techniques to reduce artefacts and increase the accuracy of MRV measurements for porous media using water flow through an OCF as an example. To reduce measurement time, a multiecho sequence was designed where each echo pair experiences a different velocity encoding. Thus, per excitation pulse for one line in k-space (Paschal 2004) a whole eight-step dual-VENC balanced encoding is done. For minimising artefacts in the fast flow regions, such as displacement errors, timings of the sequence and voxel size were carefully chosen. To improve accuracy in slow flow regions and minimising aliasing problems in fast flow regions, a dual-VENC method was used. Additionally, an algorithm is described to correct phase errors by polynomial fitting of the no-flow dataset to reduce VNR-loss.

MRV
Flow velocities can be measured using PC-MRV. Here velocity encoding gradients are applied which cause a phase evolution of moving spins proportional to the velocity while for static spins the initial phase is retained after encoding. The phase Φ of the Nuclear Magnetic Resonance (NMR) signal is described by where is the gyromagnetic ratio and ⃗ v is the velocity. ⃗ M 1 is the first moment of the gradient applied for the duration T defined by: Φ 0 results from background phase effects due to eddy currents, magnetic susceptibility and concomitant gradients. To filter out such effects at one-directional velocity, two measurements are done with different first moments ⃗ M (1) 1 and ⃗ M (2) 1 (e.g. using inverted gradient polarities) but under otherwise same conditions. Subtraction of the signal phases results in ΔΦ , which is dependent on velocity ⃗ v and 1 . Nonetheless, eddy current errors originating from the velocity encoding gradients are not eliminated. Therefore, an additional measurement is done with same parameters but without flow and the phase images are subtracted.
In this work the four-step balanced encoding was used . Z 1 to Z 4 are the complex images acquired for individual encoding steps. The velocity v is determined as The phase can only be uniquely measured over an interval of 2 . If there is flow in positive and negative direction, it is common that the VENC is thus defined as the velocity that produces a phase shift ΔΦ = , i.e.
For too high velocities, the velocity dependent phase shift can exceed ± and phase aliasing occurs. The phase shift can be calculated as ΔΦ = arg(Z 1 ) − arg(Z 2 ) = Φ 1 − Φ 2 , which is computationally costly because of two arg operations. Also additional aliasings may be introduced. Instead, the phase shift is calculated by ΔΦ = arg(Z 1 ∕Z 2 ) = arg(Z 1 ⋅Z 2 ) (Bernstein 2004). The phases in x-, y-and z-direction are calculated as follows: Analogously to Lee (1995) the standard deviation of the velocity is determined as where N is the number of encoding steps. For each voxel the Signal-to-Noise Ratio (SNR) is calculated from the magnitude I of the complex signal and the standard deviation : In magnitude images the noise distribution in regions without NMR signal is governed by a Rayleigh distribution (Gudbjartsson 1995). If the standard deviation 0 is determined in the signal free part of a magnitude image, the standard deviation , that is used to characterise the VNR, Page 4 of 20 is given by = 1.5264 ⋅ 0 . In Eq. (8), it is visible that the VENC value should be set as low as possible to achieve a low noise value v (i.e. good VNR). However, it must be noted that for VENC < v true aliasing occurs. For correction of the aliasing a reference scan with a higher VENC value is used and difference D = v high − v low is determined. The number of wrappings w is then determined as where N.I. is the nearest integer function. The corrected (unwrapped) velocity is calculated as: In case of three-directional velocity, an eight-step measurement can be done (four steps for low-and four steps for high-VENC measurement). This doubles the number of measurements in comparison with a single-VENC measurement but results in a VNR improvement by a factor It should be noted that for a large R unwrapping may fail due to large noise values in the difference D. Neglecting pulsation effects and assuming stable conditions (temperature, electrical components) during the whole measurement operation, the condition from Lee (1995)  where f(x) is the probabilty density function and erf the Gaussian error function. For a given ratio R and SNR, the fraction P of wrongly unwrapped voxels can be calculated with Eqs. (14) and (15).
An additional VNR-improvement is achieved by weighted averaging of the unwrapped low-VENC and the high-VENC maps:

MRV Sequence
A SE and PC-based MRV pulse sequence was implemented ( Fig. 1). Spatial resolution was achieved by one read and two phase encoding gradients. The spatial encoding by the readout gradients was sufficient to avoid signal wraparound for the measured samples and the chosen Field Of View (FOV). Therefore, non-selective rectangular RF pulses were used for excitation and refocusing to reduce the interecho delay and avoid signal losses. To suppress unwanted signals, spoiler gradients were applied before and after the refocusing pulses. In Shiko (2012) and Huang (2017) the velocity encoding was done before the multiecho loop. However, in this case measurements are limited to slow velocities. Otherwise large displacement errors would arise for the later echoes. To avoid this, velocity encoding was done similar to the FLIESSEN sequence (Amar 2010). Before and after each echo velocity encoding pairs with same strength but opposite direction were applied to set the phase shift to zero before the next Radio Frequency (RF) refocusing pulse. This increases the echo time but reduces displacement effects. In total eight different velocity encoding gradients were applied to achieve a full dual-VENC encoding scheme for acquisition of threedirectional velocity in one measurement. For every echo pair consisting of two consecutive echoes (odd and even) the same velocity encoding was done. Datasets from odd and even echoes are combined to suppress additional phase differences caused by different coherence pathways (for details see Shiko 2012; Huang 2017). Between velocity encoding and decoding, spins may be accelerated or decelerated (e.g. for spins moving near a narrowing). In this case, the signal phase is not completely eliminated after decoding and an acceleration-dependent phase shift remains. Combination of odd and even echoes cancels this effect (Amar 2010) (provided that the acceleration remains constant for both echoes). The velocity with VENC low was reconstructed from the first eight echoes and the velocity with VENC high was reconstructed from the last eight echoes. Associated velocities are called v 1 and v 2 . Hence, in total 16 echoes were acquired. As discussed by John (2020), the maximum displacement in number of voxels is given by where t displ is the time between the centre of velocity encoding and the centre of acquisition. Even though t displ should be chosen as short as possible it shall be noted that the gradient amplitude is limited and therefore for a given VENC the duration and delay Δ cannot be chosen arbitrarily short. Since gradients with large amplitude may produce eddy currents which worsen the image quality, these should only be applied if compensation is provided to correct for eddy currents. If the compensation is not perfect, which is usually the case for experimental data, the remaining artefacts scale with the gradient amplitude. Additionally, doing measurements with large gradient amplitude may shorten the lifespan of the gradient system. Using = 6.0 ms , a duration of 192 ms results from excitation pulse until acquisition of the last of the 16 echoes. Experimentally a transversal relaxation time T 2 = 0.41 s was determined. Thus, due to transversal relaxation the SNR of the last echo is reduced by 36.5% in comparison with the first echo. The SNR of the final (unwrapped) velocity map is mainly determined by the VENC low . Therefore, the VENC high encoding steps are measured in the latter echoes. The acquisition time per echo is 1.5 ms . The rectangular flow encoding gradients have a = 0.78 ms and a delay of Δ = 2.71 ms . The time between the centre of velocity encoding and the centre of acquisition was t displ = 3.58 ms . It should be noted that this value may differ for other sequences, since it depends on several factors (hardware, acquisition time, voxel size, VENC, flow velocity, etc.) (Table 1).

Optimised repetition time
For a given sequence, the SNR can be improved by increasing the number of accumulations or increasing the repetition time T R . For a constant measurement time, it shall be analysed at which ratio between accumulations and T R the optimal SNR improvement is achieved. Therefore, the longitudinal magnetisation is considered. After applying the 90 • excitation pulse at t 0 = 0 the longitudinal magnetisation is M z0 = M z (t 0 ) = 0 . Due to longitudinal relaxation the magnetisation after applying the first 180 • pulse (at t 1 ) is In total 16 inversion pulses are applied in the echo train. The longitudinal magnetisation directly after the k-th inversion pulse is calculated as follows: For T 1 = 2 s , which was experimentally determined, a quite low M z16 = −0.003 ⋅ M 0 results. Therefore, M z16 is neglected In eight loops (dotted rectangle) different velocity encoding gradients are applied for achieving a full dual-VENC encoding scheme. Spatial encoding is done by gradients in read-, phase1-and phase2-direction. Velocity encoding gradients with duration and time delay Δ are applied. Spins are refocused by using 180 • pulses after duration . SE is acquired after an additional . Each echo is acquired for 1.5 ms Table 1 Encoding steps (dual-VENC balanced encoding scheme) for the echo pairs. The factors a and b are set depending on VENC 1 and VENC 2

Echo pair
Page 6 of 20 in the following. We consider the time between the last inversion pulse and the excitation pulse where longitudinal relaxation takes place. With M 0 = 1 and normalisation to unit time, the longitudinal magnetisation is approximated by (Fig. 2) Plots for Eq. (20) at different T 1 values are shown in Fig. 3. For long T 1 a broad maximum is visible. However, an optimal T R value should be chosen to avoid SNR-losses, which is of particular importance for short T 1 .

Samples and flow setup
Flow measurements were done for two different samples: A flow phantom and an OCF. The flow phantom was manufactured in the university workshop by heating and turning a glass tube (inner diameter 5 mm ) in helical shape. Since more heat was added from one side during the manufacturing process, the helix is flattened on the outside of the windings. Also the cross section was narrowed. The cross-sectional area is decreased by 32% in comparison with the straight section (the value was determined on the basis of the 3D NMR data). A picture of the flow phantom is shown in the appendix (Fig. 10).
The OCF with 10 pores per inch ( 25 mm in diameter) consists of polylactic acid and was 3D-printed by a Digital Light Processing printer (Elegoo Mars, Shenzhen, China) with 2 k resolution ( 12.5 m × 12.5 m ) in XY direction and layer height (Z) of 50 m . To prevent unwanted bypass flow, seal tape was first wrapped around the OCF before it was pushed into the glass vessel. The properties of the OCF are given in Table 5 in the supplementary material. To cover similar velocity ranges, flow rates of 9.9 ± 0.1 ml min (flow phantom) and 206 ± 2 ml min (OCF) were realised using a peristaltic pump (Standard Digital Drive with Easy-Load II Pump Head; Masterflex, Vernon To minimise pulsation effects, a pulse dampener was used Hills, USA) with pulse dampener and pump tubing also from Masterflex. PVC-tubing was used for connection to the flow phantom and the glass vessel with the OCF (Fig. 4). Tap water was used as the operating fluid since this is widely used as flow model for MRV measurements. The appropriate entrance length ( L e ) was considered to provide full development of flow before reaching the flow phantom or OCF: Here D in is the inner tube diameter, Re is the Reynolds number, f is the fluid density and f is the dynamic viscosity of the fluid. The average inlet velocity u z was calculated depending on the Volumetric Flow Rate (VFR). Reynolds numbers were Re = 42 (flow phantom) and Re = 175 (glass vessel with OCF). Entrance length was L e = 12 mm (flow phantom) and L e = 245 mm (OCF). For calculating the pore-Re number, the pore diameter 5.8 mm (Table 5) was used instead of the tube diameter. Since it is difficult to estimate the average velocity in a pore, the high-VENC was chosen for u z leading to a pore-Re of 261. Therefore, laminar flow can be assumed.

Measurement parameters
All measurements were performed on a horizontal 7 T scanner (BioSpec 70/20 USR, Bruker BioSpin MRI, Ettlingen, Germany), which is equipped with a 114 mm bore gradient system (BGA 12S2, maximum gradient strength of 440 mT m , maximum slew rate of 3440 mT m⋅ms ). A horizontal 72-mm bore birdcage 1 H quadrature transceiver RF coil (Bruker BioSpin MRI, Ettlingen, Germany) was used in all measurements. ParaVision 5.1 was used as the user interface to perform the measurements. The repetition time was T R = 2.9 s (according to the maximum in Fig. 3). The gradient amplitude was set depending on the desired VENC value. Further measurement parameters are listed in Table 2.
The measurements for the flow phantom were done at VENC 1 = VENC 2 = 40 mm s to analyse the reproducibility of the measured velocity v 1 and v 2 . For improving the VNR, the dual-VENC method was used for the OCF. VENC values were VENC 1 = 20 mm s and VENC 2 = 45 mm s . The final (unwrapped) velocity is mainly determined by the VNR of the low-VENC data. Thus, the low-VENC value was chosen for VENC 1 because the first echoes experience a lower SNRloss (respectively VNR-loss) due to transversal relaxation. Measurements were done with and without flow. Considering the no-flow data, noise parameters were determined for the high-VENC data (last eight echoes) after averaging the two accumulations. For each pair of consecutive odd and even echo the magnitude images were combined (averaging) and the respective SNR was calculated according to Eq. (9). The mean value was SNR = 199.6 which was used to calculate v = 0.144 mm s (Eq. (8), number of encoding steps N = 4 ). Post-processing of the data was performed using a self-developed MATLAB script (R2021a MathWorks, Natick, MA, USA). Binary masking was based on Otsu's thresholding (Otsu 1979) for the corresponding magnitude images.

Offset correction by polynomial fitting
If offset correction is done by subtracting the no-flow velocity map, the VNR is lowered by √ 2 . For reduced VNRloss another correction method was used. Analogously to Lingamneni (1995) a polynomial fit was applied to the noflow dataset for reducing the effect of random noise. The following algorithm was used: 1. A binary mask BM 1 was calculated based on Otsu's thresholding for the corresponding magnitude images. 2. To ignore edge effects, another binary mask BM 2 was calculated. For each voxel in BM 1 it was checked if all 26 neighbouring voxels have the value 1. If yes, the considered voxel was set 1, if not it was set 0. Examples are shown in Fig. 8a and b. 3. Each no-flow velocity map was multiplied by the binary mask BM 2 . 4. For each no-flow velocity map, polynomial fitting was done. The three-dimensional polynomial can be written as follows, where the order O is defined as the highest exponent occurring in the sum:

Flow phantom
In Fig. 5, the velocity maps ( v z 1 ) for axial slice 61 (out of 120) of the flow phantom are shown. Corresponding velocity maps for v z 2 are shown in Fig. 11 in the supplementary material. For the maps without flow an offset is visible (depending on the spatial position up to 10% of the VENC). For a quantitative analysis for each of these maps the VFR was calculated in a Region Of Interest (ROI) (straight tube) by where L x and L y are the voxel size in x-and y-direction. The bar plot with the resulting VFR for v z 1 and v z 2 is shown in Fig. 6. For all cases without offset correction (odd, even, combined) the VFR is between 12 and 14 ml min and therefore does not fit the expected 9.9 ± 0.1 ml min . Better agreement is visible for the data that are corrected by subtraction of a reference scan without flow. Similar results were obtained for all axial slices. Considering only odd or only even echoes, the VFR is still a bit too low or too high. Best agreement is obtained for the combination (averaging) of odd and even echoes and correction with the no-flow measurement. However, these values are always between 10.1 and 10.2 ml min , thus overestimating the expected VFR of 9.9 ± 0.1 ml min by 1 to 4% . An analogue analysis was done for velocities v 1 and v 2 in x-and y-direction (see appendix). A plot with the VFR calculated for each axial slice can be found in the appendix Figs. 13 and 14. Also here the best agreement is visible for combination of odd and even echoes and correction with the no-flow measurement. Considering the whole 3D dataset, maximum velocities were found in the centre of the glass tube and had values of 23.5 mm s (x-/y-direction) and 17.3 mm s (z-direction). Since the cross section is slightly flattened due to the manufacture process, it is not easy to predict the flow behaviour. However, for the straight tube with a round cross section the velocity profile can be predicted. For a laminar tube flow it is v max = 2 ⋅ v mean . Regarding the tube diameter and the VFR given in Sect. 3.3, it was therefore calculated v max = 16.8 mm s . This value is 3% lower than the measured one. For laminar flow, the mean axial velocity v r at radius  For one exemplary axial slice, the calculated axial velocity is compared with the measured one (Fig. 7).

OCF
A dual-VENC measurement was done for the OCF. Complex valued images for each encoding step were combined (averaging odd and even echoes), before velocity values in x-, y-and z-direction were calculated. In a measurement under stable conditions, the velocity maps of two accumulations would only differ by their noise values. In the following, the reproducibility of a measurement is analysed by comparing two accumulations using the example of the high-VENC velocity maps. After subtracting the two accumulations from each other, the standard deviation within the pores was determined and compared with the value calculated according to Eq. (8) (for details see Sect. 3.4). The results with and without flow are compared to determine if flow effects (such as unwanted pulsation) have an additional effect. Corresponding values are listed in Table 3. Without flow (left column), similar values are visible in x-, y-and z-direction which are also roughly twice as high as the value 0.144 mm s calculated according to Eq. (8). This corresponds to the expectations because this 0.144 mm s was determined after averaging both accumulations (noise reduction by √ 2 ) and the values shown in Table 3 were calculated by subtracting both accumulations (noise increased by √ 2 ). Values resulting from the velocity maps with flow (right column) are larger, between 2.2 and 3.7 fold increased.
Offset correction was done by direct subtraction of the no-flow velocity map and by using the algorithm described in Sect. 3.5. As an example, the procedure is shown for one axial slice and one velocity component ( v y 1 ), which was calculated from the first accumulation of the measurement. The corresponding binary masks BM 1 and BM 2 are shown in Fig. 8a and b. For the no-flow velocity map an offset up to 10.5% of the VENC is visible depending on the spatial position (Fig. 8c). After applying the binary mask BM 2 , polynomial fitting for the no-flow velocity map was done by using an opensource MATLAB toolbox (D'Errico 2022). The polynomial fits with first and second order are shown in Fig. 8d and e, respectively. To quantify the goodness of the fits, the Root-Mean-Square Error (RMSE) was calculated. Using the first order fit to correct the no-flow velocity map results in Fig. 8f. Areas with similar nonzero values can be observed locally. In contrast, the values in Fig. 8g (second order fit) are more randomly distributed around zero. Calculating the no-flow velocity map from each of the two accumulations and subtract the two from each other, results in Fig. 8h. The standard deviation of the polynomially fitted maps is normalised to the standard deviation of Fig. 8h. Values can be found in Table 4. From first to second order fit a significant reduction can be observed for both, standard deviation and RMSE. For higher order fits only a slight reduction is visible. For increased VNR, both accumulations were averaged and velocity maps calculated. Polynomial fits (first to sixth order) were done for each velocity component. Corresponding RMSE values are shown in Fig. 18 in the supplementary material. The offset of the velocity maps were corrected by the third order polynomial fit. A higher order fit was not used since it does not deliver a significantly larger accuracy and to avoid the risk of overfitting. The dual-VENC method was used for unwrapping. For an Velocity v z 1 ( ) and v z 2 ( ) calculated for the low and high VENC dataset, respectively, calculated in a central axial slice after combining odd and even echoes and subtracting the no-flow data. In addition, the expected curve according to Eq. (25) is plotted Table 3 For each velocity map reconstructed from the high-VENC data (with and without flow) both accumulations were subtracted and the standard deviation was determined inside of the structure (first accumulation) an offset up to 10.5% of the VENC-value is visible depending on the spatial position (c). This map was polynomially fitted and afterwards BM 1 was applied to the resulting fits. Shown are polyno-mial fits with first (d) and second (e) order. Subtraction of (c) and (d) results in (f); subtraction of (c) and (e) results in (g). To evaluate the reproducibility, a second accumulation was acquired. No-flow velocity maps (no fitting) from first and second accumulation were subtracted from each other (h) Fig. 9 Velocity maps were offset corrected by a third order polynomial fit and unwrapped using the dual-VENC method. Unwrapped low-VENC and high-VENC maps were weighted averaged. Exemplarily, several axial slices are shown for all velocity components (x, y, z). For fixation in the glass tube, seal tape was wrapped around the OCF. No tape was present at the beginning and end of the OCF. The water filled gaps between OCF and glass vessel result in the "rings" in slice 21 and 61. Velocity values inside of these "rings" are about zero additional VNR-gain, unwrapped low-VENC and high-VENC data were weighted averaged based on the ratio R = 2.25 . Velocity maps are shown in Fig. 9.

Discussion
High accuracy of MRV measurements on porous media is essential, especially if such measurements are used for validation of CFD-simulations. For improved VNR, a new multi-echo sequence with optimised repetition time T R was proposed and a dual-VENC technique was used with weighted averaging of high-VENC and unwrapped low-VENC data. Furthermore, systematic errors were analysed, which may occur in MRV measurements, in particular displacement errors and phase errors originating from eddy currents or different coherence pathways. Techniques were described to minimise or correct such errors.

Pulse Sequence
A flow phantom with predictable flow behaviour was measured to test the new multi-echo sequence. An offset was observed and the measured VFR was overestimated if only odd or only even echoes are considered. This confirms that a combination (averaging) of odd and even echoes is necessary to correct for different coherence pathways. An offset was also observed if only the velocity maps with flow are considered. This confirms the assumption that the phase difference image (Eqs. (5) to (7)) does not eliminate all phase errors. These findings are consistent with Huang (2017). Eddy currents depend on the strength and direction of the applied magnetic gradients. Thus, the eddy currents originating from the velocity encoding gradients behave additive in the phase difference image. Therefore, it was necessary to subtract flow and no-flow measurements for correction of remaining offsets. The corresponding VFR was slightly overestimated by 1 to 4% . It can be assumed that the slow velocities, which are located at the surface of the flow phantom, were slightly overestimated due to partial volume effects. Details about such effects are described in Tang (1993); Bernstein (2004).
The measured axial velocity was compared with theoretically predicted values for laminar flow. A good agreement was visible for the flow profile. The maximum velocity determined from the measured axial velocity was slightly larger than the predicted value due to superimposed noise.
After one excitation pulse, several echoes were acquired. For the transversal relaxation time T 2 = 0.41 s it was possible to measure 16 echoes with acceptable signal loss to achieve full dual-VENC encoding. Velocity maps v 1 and v 2 were reconstructed from the first and last eight echoes, respectively. For setting VENC 1 = VENC 2 , reproducibility was observed, even though v 2 has lower SNR due to longer echo times.
Velocity encoding and decoding gradients were applied before and after each echo. Thus, the echo time was increased, respectively the SNR decreased. In Huang (2017) the interecho delay was shorter because the velocity encoding gradients were applied before the echo train. However, this approach increases displacement errors. Thus, if the same echo time is assumed for both sequences, the displacement error for the last echo differs by a factor 50. Single-point imaging sequences, as proposed by Bruschewski (2019), are able to reach a great reduction in displacement errors but in comparison with the sequence shown here the measurement time is strongly increased and not suitable for 3D imaging. The reduction in displacement errors by updating the flow encoding for each echo was already proposed in Amar (2010) for the FLIESSEN sequence. Since this sequence aimed at short total measurement time, each echo of the long echo train was measured with different spatial phase encoding. In contrast, the sequence proposed in the present work uses all echoes for flow encoding. Therefore, the minimum total measurement time is longer than in FLIESSEN or other sequences that have a long echo train length with short T R . However, for the sequence shown here a reduced spatial resolution in phase encoding directions is avoided as each image is measured at a fixed echo time.
The proposed approach assumes constant velocity and does not compensate for nonzero acceleration. To minimise the displacement errors, the time t displ between velocity encoding and acquisition was kept as short as possible. Thus, also the influence of acceleration on the signal phase is limited. Most importantly, only two consecutive echoes (odd, even) are used for the calculation of the same image as different velocity encoding steps are performed along the echo train. Thus, summation of acceleration effects (via the consecutive echoes) would not lead to incorrect velocity encoding, but to reduced signal due to intra-voxel dephasing. However, a comparison of the data measured with and without flow yielded no notable differences in the signal intensity, neither for the flow phantom nor the OCF. Therefore, acceleration effects, although being not suppressed, were negligible in the experiments. This was also proven by the fact that the measured VFR corresponds well with the expected values.
For the proposed sequence a good compromise was found to achieve reduced displacement errors, acceptable total measurement time and sufficient VNR.

Displacement Errors
For the flow phantom maximum velocities of 23.5 mm s (x-/ydirection) and 17.3 mm s (z-direction) in the centre of the glass tube were determined. Higher velocities in x-/y-direction occurred due to the reduced cross-sectional area at the windings. Displacement errors are always a problem if the direction or magnitude of the velocity vectors change. This is generally the case for turbulent flow. Since for the measurement with the flow phantom the Reynolds number was Re = 42 , laminar flow can be assumed. For the special case that only the straight tube of the flow phantom is considered (only axial velocity), in each axial slice a similar velocity should be present and therefore displacement errors would not have a large impact on the velocity maps. However, large displacement errors would negatively influence the measured velocity fields in the windings of the flow phantom.
In this measurement a voxel size of 0.30 mm was chosen. Considering the maximum velocity in x-/y-direction and the time between velocity encoding and readout, with Eq. (17) a displacement = 0.28 voxels can be calculated, which corresponds to 0.08 mm . Even though this value is below one voxel, displacement errors may occur, namely for spins located at the border of a voxel. However, it should be noted that the estimated displacement error occurs for the maximum velocity and therefore in most voxels the displacement error is much lower.
For the OCF, the maximum velocity was 44.7 mm s and the voxel size 0.45 mm . Therefore, displacement errors up to displacement = 0.36 voxels (corresponding to 0.16 mm ) were present. This displacement error is small in comparison with the pore size of 5.8 ± 1.9 mm . Additionally, the displacement error will be lower in most regions of the OCF as the displacement value of 0.36 voxels was calculated for the maximum velocity.

Dual-VENC
When applying the dual-VENC technique, often a ratio R between two and three is chosen (Schnell 2017;Aristova 2019;Callahan 2020;Nakaza 2022;Mahinrad 2022). In the present study, an optimal R with negligible wrong unwrappings ( m = 5 , P = 6 ⋅ 10 −7 ) was calculated. Therefore, the optimal R = 22 could be chosen (Eqs. (14) and (15), N = 8 due to subtraction of flow and no-flow velocity maps). For this calculation the mean SNR was assumed. However, voxels located at the edge areas are only partly filled with fluid and therefore have a lower SNR. For example, if a voxel is only filled to a quarter, the calculated optimal ratio would be R = 5.4 . Due to large specific surface area, especially in porous media such edge voxels should be taken into account. Moreover, the calculation for such optimal R is only true for fully reproducible measurements. For a real measurement changes may occur due to fluctuations of the fluid field or changes of the temperature or electrical components of the NMR system. Especially for 3D measurements, which usually run over longer time (the measurement for the OCF needed 20.6 h ), such variations may be a significant factor. In Table 3, the standard deviation of the velocity maps with flow was up to 3.7 fold increased in comparison with the ones without flow. This indicates the presence of larger variations during the measurements with flow, which are presumably due to pulsation effect caused by the peristaltic pump used. Therefore, a rather low R = 2.25 was chosen. However, a more systematic analysis of these effects should be done to achieve a better understanding and reduce flow induced phase fluctuations to such a level that similar v values are obtained as in no-flow measurements. This should preferably be done by improving the pulsation dampening (Blythe 2015) or using a pulsation-free pump system. Alternatively, triggering techniques could be used that are widely applied in medical applications of MRV (Selby 1992;Brown 2014). Then larger R values could be chosen to use the full potential of dual-VENC encoding. Each slice (and each velocity component) was visually inspected to determine if wrong unwrapping was present. Therefore, inside of the pores it was searched for voxels, which show a large difference (about 2 ⋅ VENC ) in comparison with the surrounding ones. None of such voxels was found.
Flow through porous media, like the used OCF, usually has a large dynamic range. Therefore, the dual-VENC approach is of great use to prevent aliasing in fast flow regions but still realise a sufficient VNR in slow flow regions. For the chosen ratio between high-and low-VENC, the VNR was improved by R = 2.25 in comparison with a single-VENC measurement. Instead of doing a dual-VENC measurement, in the same time two accumulations of the high-VENC measurement could be done and averaged, which increases the VNR by √ 2 . Therefore, the effective VNR-improvement of the dual-VENC technique with R = 2.25 , is 2.25 √ 2 =1.59. An additional VNRimprovement was achieved by weighted averaging unwrapped and high-VENC data, as described by Eq. (16). For R = 2.25 this is 9.4%.

Offset-Correction by polynomial fitting
For both samples (flow phantom and OCF) an offset was observed in the velocity images underlining the necessity of using a no-flow measurement for corrections. However, if velocity maps of flow and no-flow measurements are subtracted, this decreases the VNR by √ 2 if the same v values occur in datasets measured with and without flow. Instead, a polynomial fit can be done for the no-flow map and be subtracted from the flow map to decrease the influence of noise. However, it should be noted that for inaccurate fitting systematic errors may be introduced. Porous media, like the used OCF, usually have a large specific surface area. An algorithm was used (see Sect. 3.5) to ignore edge voxels, which could worsen the quality of the fit (e.g. due to partial volume effects). Thus, only voxels inside of the pores were considered for the polynomial fit. For the correction of v y 1 by a first order polynomial fit, locally significant systematic errors occur (Fig. 8f). In contrast, when using a second order polynomial fit, velocity values are more equally distributed without large areas that could indicate systematic fitting errors ( fig. 8g). Analogously, also the RMSE decreased for higher order fits (see Table 4 and Fig. 18). Both offset correction methods -polynomial fit and direct subtraction -were evaluated in relation to the standard deviation. For correction with the first order polynomial fit the standard deviation was 6% larger in comparison with direct offset correction for each voxel. This can be explained by the systematic errors introduced by the fit. For higher order fits the influence of such systematic offsets was lower. However, the maximum improvement of √ 2 was not achieved. For the second order fit the standard deviation was 9% lower than for the direct offset correction and for the sixth order even 12% . Therefore, the VNR can be improved if the phase offset is corrected by using at least a second order polynomial fit. One should note that for a too high order overfitting could take place and have a negative impact.
The fact that the VNR improvement was considerably lower than the improvement by √ 2 expected from theory is consistent with the findings of larger v values for flow data than for no-flow data (see Table 3). Since the VNR was thus dominated by the data measured with flow, the VNR improvement by polynomial fitting of the no-flow data will increase if the flow induced phase fluctuations are reduced to a minimum, either by better pulsation dampening or avoiding the use of a peristaltic pump.

Conclusion and outlook
In this publication, methods were described to achieve a good accuracy for MRV measurements and reduce systematic errors for slow flow of water through an OCF with 10 pores per inch. A new multi-echo sequence with dedicated timings was proposed to achieve a compromise between VNR, total measurement time, spatial resolution and displacement errors. A dual-VENC encoding scheme was used and the repetition time was optimised for improving the VNR. Phase errors originating from different coherence pathways could be eliminated by averaging the datasets of odd and even echoes. Phase errors resulting from eddy currents could be eliminated by subtracting flow and no-flow velocity maps. The VNR-loss could be reduced if the no-flow map is fitted with a polynomial of second or higher order.
In the future, the presented MRV sequence as well as the correction techniques could be used for measurements to cross validate CFD simulations. In this context, other effects should also be taken into account that have not been dealt with in this publication. These are for example partial volume effects, which occur at the surface and therefore reduce the measurement accuracy especially for porous media with large specific surface area. One technique for correcting partial volume effects was described in Tang (1995). To reduce the measurement time of MRV, different acceleration methods could be used like compressed sensing (Tayler 2012;John 2020) or parallel imaging (Peng 2010). The choice of the optimal ratio R between high-and low-VENC should be analysed systematically to achieve the best VNR-improvement at a negligible number of wrong unwrappings. Here, Eqs. (14) and 15) provide only an upper limit for R. However, a considerably lower ratio R may be required because of unwanted experimental changes during the total measurement time (pulsation of the flow, temperature changes, changes in electrical components) as well as lower signal intensities (lower SNR) in voxels that are only partially filled with the fluid. See Figs. 10,11,12,13,14,15,16,17,18 and Table 5.

Parameter Value
Open porosity 0.77 Pore diameter a 5.8 ± 1.9 Window diameter a 3.3 ± 0.9 Specific surface area b 521.3 Fits were done for different orders of the polynomial. Before fitting, velocity maps of both accumulations were averaged. The velocity maps corrected by a sixth order polynomial fit were used to calculate the maps shown in Fig. 9