FLUIDS 3 , 044102 ( 2018 ) Fluid flow in a porous medium with transverse permeability discontinuity

Magnetic resonance imaging (MRI) velocimetry methods are used to study fully developed axially symmetric fluid flow in a model porous medium of cylindrical symmetry with a transverse permeability discontinuity. Spatial mapping of fluid flow results in radial velocity profiles. High spatial resolution of these profiles allows estimating the slip in velocities at the boundary with a permeability discontinuity zone in a sample. The profiles are compared to theoretical velocity fields for a fully developed axially symmetric flow in a cylinder derived from the Beavers-Joseph [G. S. Beavers and D. D. Joseph, J. Fluid Mech. 30, 197 (1967)] and Brinkman [H. C. Brinkman, Appl. Sci. Res. A 1, 27 (1947)] models. Velocity fields are also computed using pore-scale lattice Boltzmann modeling (LBM) where the assumption about the boundary could be omitted. Both approaches give good agreement between theory and experiment, though LBM velocity fields follow the experiment more closely. This work shows great promise for MRI velocimetry methods in addressing the boundary behavior of fluids in opaque heterogeneous porous media.


I. INTRODUCTION
Flow in porous media is important in many areas of science and technology including biology, soil science, reaction engineering, waste treatment and separation science [1][2][3].Flow in a composite region, where flow in a porous medium is coupled to the flow in a region free of a porous medium, is of particular interest since it occurs in a wide variety of practical applications.Some examples of such practical systems in industry include thermal insulation and crude oil extraction, while blood flow in large vessels represents a biological foundation.
Usually, permeability in a porous medium varies only microscopically, or on a length scale of a few pore sizes.Hence, within a macroscopic volume-averaging element, which is much greater than a few pore diameters, the permeability is considered to be uniform.In a composite region, such that a macroscopic volume-averaging element with uniform properties cannot be defined, the permeability varies, thus complicating the definition of volume-averaged flow equations and relevant boundary conditions and, consequently, the analytical flow description in such systems.Beavers and Joseph used the Darcy law to describe flow in the porous medium [4].They circumvented the formal mathematical treatment of the fluid momentum transport across the interface by constructing an empirical so-called hydrodynamic boundary condition [4].Further steps were made when non-Darcian effects in the composite medium were accounted for by using the Brinkman-Forchheimerextended Darcy equation in the porous medium [5].The boundary condition at the interface with clear fluid was studied by Ochoa-Tapia and Whitaker [6,7], who were able to match the Brinkman-extended Darcy law to the Stokes equations by requiring discontinuity in the stress but retaining continuity in velocity using a sophisticated averaging procedure [6,7].The boundary condition was termed the stress jump boundary condition and was further utilized by Kuznetzov [8] for derivation of the governing flow equations for fluid flow in the channel bound to a cylindrical porous medium.
The refinement of these theoretical formulations requires experimental flow data.Major velocimetric methods, such as particle image velocimetry and laser Doppler anemometry, cannot satisfactorily provide measurements of flow fields in composite media because of the mismatch in the refractive indices between the medium and the fluid.Therefore, velocimetry data, in the interfacial regions between a bulk fluid flow and a flow in a usually opaque saturated porous medium, are extremely limited.
Pulsed field gradient nuclear magnetic resonance (NMR) methods have been shown to be extremely useful in studying fluid transport and flow in opaque porous media [9][10][11][12][13][14]. Spatial encoding provided by magnetic resonance imaging (MRI) methods allowed for localization of flow parameters including at the composite, or interface region.Current hardware developments allow spatial mapping of velocities with resolution up to 10 μm [15].Thus the technique is very well suited for measurements of transport phenomena between bulk fluid flow coupled to a saturated porous medium.
In this work, MRI velocimetry methods are applied to measure the radial velocity profiles in a clear channel flow bounded by a saturated porous medium.To validate experimental results obtained by the MRI methodology, velocimetry data are compared to the theoretical flow fields derived in [8] and using pore-scale lattice Boltzmann modeling (LBM) of flow in a clear channel bound by a saturated porous medium.The LBM approach is also used to see if it can provide a numerical foundation for the description of coupled flows while avoiding the need for hypothetical assumptions about the interface used in previous theories [4,6,7].

II. BACKGROUND AND THEORY
In homogeneous porous media, or media where permeability of the medium is uniform within a chosen averaging volume, the flow is usually governed by Darcy's law where μ is the local fluid viscosity, K is the local permeability of the porous medium, u D is the Darcy velocity vector, defined as an averaged velocity over a representative volume of the porous medium, and p f V is the pore averaged pressure.Analytical description of flow for a composite region has been attempted in a number of simple geometries.For one-dimensional flow through a channel, of width h, bounded by a saturated porous medium on one side, and by an impermeable wall on the other side [Fig.1(a)], the resulting velocity fields are usually estimated from two different theoretical approaches.In the first approach, often termed Brinkman's equation [16], the interface between the saturated porous medium and the adjoining channel is treated as a part of a continuum.A macroscopic shear term is added to Darcy's law to account for the velocity gradients present at the interface where u V is a volume-averaged velocity for a local representative elementary volume located away from the interface, μ is the local viscosity, μ is the effective viscosity, and μ = μ in the original article [16].The rest of the notation is kept the same as in Eq. (1).In the second approach, developed in [4], velocity profiles in the channel are found by solving the Navier-Stokes equation with a no-slip boundary condition imposed on the velocity at the impermeable wall and with a slip, or hydrodynamic, boundary condition imposed on the velocity at the interface with porous medium.The hydrodynamic boundary condition relates the interfacial shear to the slip velocity at the interface according to where α, the empirical slip coefficient, is a dimensionless quantity that supposedly depends on the structure of the porous medium, y = 0 + indicates that the velocity gradient is evaluated at the interface and towards the adjoining channel in the direction transverse to flow, u A is the area-averaged component of velocity in the x direction (i.e., along the flow), u A,i is the area-averaged interfacial or slip velocity evaluated at y = 0, u D is the Darcy velocity, and K is the permeability of the porous medium.The resulting velocity distribution for this model is given elsewhere [4,17,18].The validity of both approaches in the interpretation of experimental flow data has been extensively studied [19].For a complete analysis of the application of either approach to the description of flow in spatially heterogeneous model porous media the reader should refer to [19][20][21].The major problem in these two approaches was the overdetermination of the physical problem in an attempt at the continuous matching of shear across the interface [21].This issue was overcome in the approach proposed in [6,7], where the use of sophisticated averaging technique matched the Brinkman-extended Darcy law to the Stokes equations by requiring a discontinuity in the stress but retaining continuity in velocity.The boundary condition was termed the stress jump boundary condition and the use of discontinuity in the stress at the interface circumvented the overdetermination of the physical problem as was pointed out by Neld [22].

A. Fully developed axially symmetric flow in a clear channel bound to a saturated porous medium
The above approach was utilized by Kuznetzov to derive governing flow equations for fluid flow in a clear channel bound to a cylindrical porous medium [8].The model geometry used to obtain governing flow equations for this type of flow consists of a cylindrical porous medium in which permeability is discontinuous at a channel coaxial with a central axis of the medium.The radius of this open channel, R, is much larger than the average pore size of the medium.The flow is in the direction parallel to the central axis of the medium.The pictorial representation of this model flow with relevant definitions is given in Fig. 1

(b).
The problem was treated before [8] and basic flow equations are stated here rather than derived.Closely following [8], we see that in the open fluid region, or η region, the equation of fluid motion expressed in dimensionless variables is given by while in the porous medium, or ω region, the dimensionless expression holds, where where r is the radial coordinate (in m), r is the dimensionless radial coordinate, R is the radius of the open channel (in m), R * is the outer radius of the cylindrical porous medium (in m), Da R is the Darcy number for a cylindrical channel, R is the dimensionless pressure gradient in the direction transverse to flow for a cylindrical channel, μ f is the fluid viscosity (in kg m −1 s −1 ), μ eff is the effective viscosity in the Brinkman term (in kg m −1 s −1 ), U is the reference velocity (in m s −1 ), p f f is the intrinsic average pressure (in Pa), and x is the streamwise coordinate (in m).The following boundary conditions are utilized: The solution of Eqs. ( 4) and ( 5) with the boundary conditions ( 7) is given by where the coefficients A 3 and B 3 are given by and respectively, with To evaluate the dimensionless interfacial velocity u i , the stress boundary condition given by is used [6,7].Then the expression for u i becomes complicated and is given by [8]

B. Pore-scale LBM of axially symmetric flow in a clear channel within a saturated porous medium
To avoid difficulties involved with a specific treatment of the interface and an effective viscosity parameter (as described above), one can apply pore-scale LBM to treat this problem.One usually starts from the microscopic approach that is based on the Navier-Stokes equations governing the dynamics of fluid flowing both in the pure fluid region and in the heterogeneous porous region on the pore scale.This is equivalent to the single domain formation on the microscopic scale across the whole domain.One of the advantages of this approach compared to the approaches described in Sec.II is the elimination of any specific treatment of the interface.However, the challenge is in the accurate representation of the complicated pore structure on the pore-scale level.This requires experimental and demanding computational resources.Experimentally, one needs to use x-ray microcomputed tomography (XμCT) to obtain and to further reconstruct the complicated 3D structure of a porous medium under study.Computationally, the LBM method has been established as an efficient algorithm for simulating fluid flows, especially for flow of fluids through complex porous media [23].The reconstructed 3D digital structure, where each voxel represents a portion of the original porous material, is usually mapped into computational lattice nodes in LBM [24].
In this study, LBM employs the single-relaxation-time lattice Bhatnagar-Gross-Krook collision model [25,26] that can be written as where x denotes the position of lattice node, τ denotes the relaxation time, and f i denotes the particle distribution function which streams along the discrete lattice velocity vector e i ; the right-hand side of Eq. ( 15) represents the rate of change of f i resulting from collisions between the incoming f i and the equilibrium particle distribution function f eq i .The relaxation time τ that controls the rate of approach to equilibrium is chosen to be unity to reduce the viscosity dependence [27] in the whole fluid simulation.Here δt is assumed to be unity.The function f eq i depends on only the fluid local density ρ and velocity u and has the form (the three-dimensional d3q19 model) where The fluid density ρ and velocity u on each lattice node are evaluated through the moments of the particle distribution functions with respect to e i , A Chapman-Enskog procedure can be applied to Eq. ( 15) to derive the macroscopic governing equations of fluid flows with second-order accuracy.They are given by the continuity and Navier-Stokes equations Therefore, through pore-scale lattice Boltzmann modeling we can directly simulate the behavior of fluids flowing inside both the clear fluid region and the adjacent porous medium.

A. NMR methods
NMR and MRI velocimetric methods and their applications to study fluid flow have been extensively reviewed [9,12,15,28].In the present work, the motion encoding was achieved by a phase method [29].Briefly, a pair of field-gradient pulses of different effective polarity is incorporated into an NMR pulse sequence in addition to rf excitation which produces transverse magnetization.These field-gradient pulses of duration δ and magnitude g are usually separated by a time delay .The application of the first gradient pulse changes the local magnetic field according to B = B 0 + rg, where B is the local magnetic field and r is the position vector.If a collection of magnetically active spins was at position r 0 flowing with velocity v, then transverse nuclear magnetization undergoes a phase shift of γ gr 0 δ.During a time interval , the nuclei move to a position r = r 0 + v .As a result, at the end of the second gradient pulse the transverse magnetization has the phase shift of γ g(r 0 − r)δ, or −vγ gδ .The transverse magnetization produced by this collection of magnetically active spins averaged through a volume element d 3 (r 0 − r) can be expressed as where P(r 0 − r; ,r) is the probability for an ensemble of nuclei to be displaced by r 0 − r from r 0 during , or a propagator, and q = γ δg is a wave vector conjugate to the space of spin displacement r 0 − r.For the case of unrestricted isotropic diffusion superimposed on uniform flow Eq. ( 22) can be solved and a solution is given by [30] E(q,r, where E 0 is the initial intensity or the initial spin density, D is the apparent dispersion coefficient, and v(r) is the local velocity.Analyzing experimental intensities for a few values of γ gδ by Eq. ( 23), one can estimate values of E 0 , v, and D.
The combination of this method with standard MRI methods [28] allows for spatial mapping of velocities under flow conditions.The mapping is usually done by three-dimensional data acquisition where two dimensions are used for spatial encoding and the third dimension encodes the flow.Specifically, a set of planar images in the direction transverse to flow is acquired for a few values of field gradient γ gδ.In such a set of images, the intensity of an individual area element, or a pixel, is a function of γ gδ and obeys Eq. ( 23).Therefore, velocity values can be determined for each individual pixel by fitting the field-gradient-dependent intensity of the pixel to Eq. ( 23).The intensity of each pixel is then replaced by a corresponding velocity value to yield two-dimensional velocity maps.The same approach is used to produce spin density E 0 and dispersion D maps.Two-dimensional velocity maps can be further analyzed to produce radial velocity profiles.The analysis is simply based on azimuthal averaging of velocities extracted from the maps, assuming axial symmetry of the flow.A complete description of this method is given elsewhere [31].
Timing diagrams of velocimetric MRI protocols used in the present work are displayed in Fig. 2. The spin-echo protocol shown in Fig. 2(a) was used to characterize fast fluid flow in the open channel and in the vicinity of the interface with the porous medium.To characterize fluid flow in the porous medium where T 2 is greatly shortened by diffusion of fluid in the porous medium, the stimulated-echo protocol, shown in Fig. 2(b), was preferred.Since this protocol allows for z storage of magnetization during the NMR flow experiment, the observation time in this protocol is limited only by spin-lattice relaxation, T 1 .Since the dephasing of spin coherence is not a limiting factor in this protocol, the NMR data can be measured for larger values of , the time between field-gradient pulses.As a result, the measured velocity fields will reflect the motion of spin particles that traveled two or three pores before data were sampled [10].This is beneficial for the characterization of temporally averaged percolating flow which occurs in a porous medium.
For both spin-echo and stimulated-echo experiments, planar sample cross sections were imaged transverse to the direction of flow.Flow encoding was performed for 10-12 and for 6-8 different values of applied field gradient q in spin-echo and stimulated-echo experiments, respectively.Image size was kept at 256 × 256 pixels for spin-echo and 128 × 128 pixels for stimulated-echo experiments, thus the total three-dimensional data acquisition required from 10 to 15 and from 18 to 24 h, respectively.The larger experimental time in the latter experiment is explained by longer (approximately 100 times) values, the time between field-gradient pulses, and larger signal averaging.The field of view (FOV) was 3.0 × 3.0 cm 2 in both experiments.
Bulk displacements were measured using the timing diagram displayed in Fig. 2(c).The spectra were measured with 256 flow encoding gradients at different values of flow encoding times.The resulting two-dimensional data were processed by Fourier transformation in spectroscopic and q directions with no zero fill in either direction of the Fourier transforms.The resulting spectrum in the q direction (second dimension) displays all coherent displacements that occur during the chosen flow evolution time.Details of this approach can be found in [28].
All experiments were performed with a 25-mm "birdcage" 1 H resonator tuned to a proton resonance frequency of 600 MHz.The field drift was negligible during experiments; therefore, no deuterium lock was used.Shimming was performed directly on the sample since a sufficient signal-to-noise ratio was achieved after one scan; a 20 Hz linewidth of a 1 H single water peak was readily achieved after successful shimming.The lengths of π/2 hard pulse 1 H pulse and sinc-shaped selective pulse were 20 μs and 1.5 ms, respectively.

B. Experimental setup
A reticulated polyurethane foam was used as a model porous medium.The material is designed for water filtering and is characterized by a spongelike structure where a nearly rigid polymer network comprises large and open pores.The pore count is approximately four pores per cm (d av p = 0.25 cm) and permeability data for a similar material can be found elsewhere [5].A sample porous medium was manufactured from bulk material in the shape of a cylinder which was 55 mm in length and 24 mm in diameter and was used to characterize displacement in the bulk of this porous medium as it happens well far away from the interface with an open channel.Another porous medium was manufactured from this material with exactly the same shape, though a cylindrical channel of 7 mm in diameter was cut through the sample coaxial with its central axis.To prevent fluid momentum transport across the interface a rigid tube [(7 ± 1) mm outer diameter and (5 ± 1) mm inner diameter] made from impermeable polymer (G10) was inserted into the channel.Bulk flow measurements were performed in the porous medium without the channel, while MRI velocimetry measurements were performed in the porous medium with the channel and with the polymer tube inserted into the channel.
Boiled distilled water with a trace amount of antiseptic was used in this study as a model fluid.Water was boiled to minimize bubble formation in the porous medium during flow experiments and traces of antiseptic (sodium azide) were also added to prevent bacteriological contamination of the water.The porous media were incorporated into a custom-made plastic column with outlets on the top and the bottom of the column to provide circular flow through the porous medium.The column was mounted in the middle of the superconducting magnet and Tygon tubing was used to connect the column with a water reservoir and a pump.The desired flow rate was delivered by the pump placed between the entrance to the column and water reservoir.The flow in the column was directed from the bottom to the top in order to facilitate bubble removal from the porous medium during the experiment.

C. X-ray microcomputed tomography
The sample was imaged using a high-resolution x-ray 3D computed tomography microscope instrument, model VeraXRM-510 (manufactured by Xradia Inc, Pleasanton, CA, USA).The voxel resolution was 5.5 μm.The 3D structure of the porous medium was further constructed by stacking together individual 2D images.

IV. RESULTS AND DISCUSSION
The character of transport in a bulk porous medium is usually reflected in a spectrum of molecular displacements P av (r, ) measured over a time interval using the generalized NMR protocol displayed in Fig. 2(c).The shape of the spectrum of molecular displacements strongly depends upon the quantity v/d p , where is the time allowed for displacements to occur, v is the interstitial fluid velocity, and d p is the averaged pore size [27; 28; 29].For small v/d p , the propagator is very asymmetric and reflects the instantaneous velocity distribution in channels and pores of the medium.For large v/d p , the propagator usually becomes nearly Gaussian.For intermediate v/d p , the displacement spectrum can be very complex with multimodal propagators reflecting different modes of transport in the porous medium [9,10].
The character of transport in the bulk porous medium sample used in this study can be understood from examining the temporal dependence of the longitudinal probability displacement distributions of fluid molecules in a single-phase incompressible flow of water in a reticulated foam cylinder at Re ≈ 13 determined using the average pore velocity over a length scale of d av p . Figure 3 z is the net fluid particle displacement over time in the direction of flow.It was found that at = 1000 ms the displacement probability distribution can be simulated with a bimodal Gaussian propagator, as shown in Fig. 3(b).The first peak is located at around zero value and corresponds to fluid molecules located in low-velocity zones such that displacements occurring there during the measurement time are small compared to the pore size.Using a velocity of 0.0026 cm/s and a time of 1 s, we estimate that the length scales of these small local areas are on the order of 4.3 × 10 −2 cm.The average pore size of the reticulated foam, being 2.5 × 10 −1 cm, is approximately one and half orders of magnitude higher.Since the presence of such small pores can hardly be expected in this material, the appearance of the first peak can be plausibly explained by the displacements of water molecules inside the polymeric material that comprises the foam.The fact that water penetrates into the polymeric core of the foam is also evident from the swelling of the foam upon prolonged exposure to water.Similarly shaped propagators were also observed for flow in a bed of permeable spheres [32].The faster displacing peak can be described by Gaussian distribution and this can be explained by reaching the → ∞ limit where diffusing molecules probe the connectivity of the pore space only.In regions where the displacement probability distribution lies between short-time and long-time behavior, the transport appears to be intermediate in character between dispersion and diffusion and between convection and dispersion, leading to a complicated form of the propagator.Similar effects have been reported before [10].The bimodal Gaussian behavior is apparent for all < 1000 ms and, in principle, reflects similar transport character in the bulk reticulated foam within the probed flow evolution times.
While insights into bulk transport properties within open-cell porous media can be deduced from volume-averaged velocimetry in a fairly straightforward fashion [33], the near interface fluid transport requires more demanding velocimetry experiments to fully assess the influence of the fluid momentum transport across the interface on the boundary behavior and the depth of the boundary layer.A 3D visualization of the solid matrix of the porous medium used in this work, based on reconstruction of the segmented XμCT images is displayed in Fig. 4(a) and confirms the complexity of the interconnected pore network structure.Structural variations contributing to the varying porosity of the medium are also apparent.However, the porosity as a function of depth or radial position can A 3D visualization of the cylindrical porous foam with a coaxial clear-cut circular channel reconstructed using a series of the segmented 2D XμCT images displayed in (a).The solid matrix is shown in yellow, while transparent designates pore spaces.The radius of the clear channel in this porous medium produced specifically for XμCT was approximately 0.5 cm.(c) Azimuthally averaged porosity profile produced from the reconstructed porous medium displayed in (a).be determined using segmented binary XCT images.Furthermore, the porosity in the radial direction can be determined through azimuthal averaging.Figure 4(b) shows an azimuthally averaged porosity profile that was obtained using an ensemble average for ten different cross-sectional images.The porosity inside the clear fluid region attains a value of unity, while the porosity inside the adjacent porous region fluctuates, which indicates the expected variability and complexity of the pore space, especially in the near interface region.
Two-dimensional flow-encoded images of flow in the channel with impermeable and porous interfaces as a function of the applied field-gradient strength q are displayed in Fig. 5. Images are characterized by sharp contrast and a high signal-to-noise ratio.Distortions in image intensity on the outside edges are induced by the rf inhomogeneity of a birdcage coil [34]; however, the effect of these distortions on velocity measures using the phase method are usually minimal.The intensity of each individual pixel in each series of images is modulated according to Eq. ( 23).As a result, velocity gradients cause additional phase shifts that lead to the formation of alternating circular rings in image intensity.As the strength of the applied field gradient q increases, the rings become more compact as the phase gets accumulated, and the number of phase cycles is increased toward the center of the 044102-10 FIG. 5. Phased flow images as a function of the applied field gradient q (a) in the absence of the fluid momentum transport across the interface (the rigid tube is inserted into the channel) and (b) with the momentum transport across the interface (the rigid tube is taken out).Flow encoding was performed by spin-echo protocol in the direction of flow for a few cross sections across the height of the porous medium.The raw matrix size was 256 × 256 and the FOV was 3 × 3 cm 2 with a slice thickness of 0.1 cm.The flow encoding time was kept to 8 ms as no outflow effects influenced the shape of velocity profiles when a volumetric flow rate Q of 1.4 cm 3 /s was used.
image.This behavior is most apparent at the highest field gradient for the flow in the clear channel in the absence of the fluid momentum transport across the interface [Fig.5(a)], where highly symmetric phase rings are the consequence of well-developed Poiseuille flow.When fluid momentum transport across the interface occurs, the symmetry in the phase ring structure gets distorted [Fig.5(b)] as the presence of momentum transfer across the permeable interface creates additional velocity gradients.
Two-dimensional velocity maps are displayed both in the absence [Fig.6(a)] and in the presence [Fig.6(c)] of the fluid momentum transport across the interface.A Gaussian distribution of dynamic displacements over a chosen time interval was confirmed by bulk displacement measurements (Fig. 3) and was used for velocity calculations as described in Sec.III.The radial velocity profiles were produced from these maps by azimuthal averaging with a fixed bin size in the radial direction (0.226 mm), as explained in [31].The error bars, calculated as standard deviations in each bin, are also shown for each averaged velocity value.
As can be seen from Fig. 6(b), the radial velocity profile of flow in the open channel is nearly parabolic in the absence of the momentum transfer across the interface with the porous medium.Deviation from Poiseuille flow is explained by the incomplete flow development.However, the fluid velocity at the boundary is zero, as expected for the flow coupled to the impermeable interface.The volume-averaged velocity distribution in the porous medium is consistent with Darcy's law (1).As can be seen from Fig. 6(b), the velocity increases and then reaches a constant value within the porous medium.The velocity at the tube wall is also equal to zero, as expected from the no-slip boundary condition.The boundary velocity at the column wall is obscured by larger error induced by the reduced signal-to-noise ratio closer to the edges.Nevertheless, the overall flow behavior is consistent with theoretical predictions [35,36].Flow behavior in the presence of the momentum transfer across the interface with the porous medium drastically changes the overall flow behavior as can be seen from the radial velocity profile displayed in Fig. 6(d).In the absence of the rigid wall the transport of fluid momentum across the interface occurs and results in a slip in velocity at the boundary with the porous medium.The velocity reaches a constant value at r greater than R, the radius of the channel.This flow behavior is consistent with theoretical considerations outlined in Sec.II.
To establish a steady flow field within the samples during LBM simulations a constant pressure gradient that ensures a laminar flow was applied along the longitudinal direction.When a rigid tube is inserted into the channel, flow fields inside the middle pipe and outside the pipe in the porous sections are independent of each other as the two systems are separated by an impermeable wall as shown in Fig. 7(a).As a result, the flow inside the middle pipe is a classic Poiseuille flow where an analytical solution for the velocity profile is explicitly defined.A comparison of the simulated result with the analytical solution is shown in Fig. 7(b).It should be noted that the simulated velocity profile inside the middle pipe region is in very close agreement with the analytical solution, thus indicating that the resolution of spatial discretization based on reconstructed x-ray computed tomography (XCT) images is sufficient to resolve a given Poiseuille flow.As expected, the velocity inside the impermeable wall falls to zero.The azimuthally averaged velocity profile inside the porous region during this flow 044102-12 is also displayed in Fig. 7(b).Velocities in the profile vary across the porous region with these variations corresponding to fluctuations in the sample porosity in the radial direction, as shown in Fig. 4(b).Furthermore, these variations in velocities are expected, as in this work, as velocities inside an individual pore are solved directly and should be different from the volume-averaged uniform velocities profile derived from the phenomenological Darcy law (1).Moreover, the velocity attains its highest value in the profile at around r/R ∼ 2.5 and this corresponds to the highest porosity achieved after the azimuthal averaging, as evident from the radial porosity profile displayed in Fig. 4(b).
The LBM simulated and normalized two-dimensional velocity map and a corresponding azimuthally averaged velocity profile for the flow in the presence of the fluid momentum transport across the porous interface are displayed in Figs.7(c) and 7(d), respectively.The flow field was imposed by applying the same pressure gradient as in the case for the flow where momentum transport across the interface was absent.The results of the pore-scale direct modeling of this flow field indicate that the effects of the viscous shear layer appear to penetrate into a boundary layer region of the porous medium within, as shown in Fig. 7(c).In this boundary layer the velocity rapidly varies, resulting in a slip at the interface, as can be clearly seen from Fig. 7(d).Although the boundary condition on the tangential velocity component of the free fluid at the interface is not explicitly specified, the pore-scale modeling clearly predicts slip at the permeable interface with the porous medium.Further, it can be observed from Fig. 7(d) that the slip velocity value of the free fluid at the interface is considerably greater than the averaged velocity value within the Darcy section.This indicates that the flow rate through the free fluid region adjacent to a permeable porous material under the same condition is significantly greater that that of the fluid region adjacent to a nonpermeable wall.

A. Spatial averaging and mass balance assessment
To assess the influence of structural fluctuations in the porous bed on velocity fields measured in the MRI flow experiment, planar velocity maps were collected as a function of slice thickness at a flow rate of Q = 2.25 cm 3 /s.Azimuthally averaged radial velocity profiles determined from this set of experiments are displayed in Fig. 8. Figure 8 shows that radial velocity profiles are independent of the chosen slice thickness used in this series of flow experiments.The independence of velocity distributions on the slice thickness is probably explained by the relative structural homogeneity of the porous medium.Moreover, the result also suggests that any fluctuations in velocities that might be induced by structural inhomogeneities are not measurable during the chosen flow encoding time = 8 ms, in fact the shortest flow encoding time used in this study.
In addition, the velocity profile obtained from the 1 mm slice was used to check the mass balance in an MRI flow experiment.The calculated flow rate obtained by radial integration of the 1 mm profile was Q = 2.5 cm 3 /s compared to Q = 2.25 cm 3 /s determined by timed collection.The difference between the two values can be explained by difficulties in determining correct values of velocities near the column wall induced by the poor image quality at the edge due to rf inhomogeneity of the coil, but in all the difference between the measured and the calculated mass transport was within 10%.

B. Temporal averaging in MRI flow experiments
To explore the influence of percolating motion of a particle in a porous medium (an observed spin particle is allowed to travel over two or three pores before velocities fields are measured) on the fluid velocity in the porous medium and in the vicinity of the porous interface with the open channel, the velocity data were collected flow at = 8 ms and at = 500 ms and then were compared.The measurements were first performed with the impermeable tube inserted into the channel to decouple flow in the channel from flow in the porous medium.No difference within experimental error [Fig.9(a)] between the velocity profiles measured at = 8 ms and = 500 ms was observed.The absence of temporal effects in velocity profiles can be explained as follows.When flow in the clear channel is no longer coupled to the porous medium and is in principle fully developed at Q = 1.4 cm 3 /s no difference in velocity fields in the clear channel associated with the entry effects should be observed by the time the fluid reaches the detection region.This is supported by a simple estimate of the entry length L h of 10.692 cm for Re ch = 324 determined using the diameter of the clear channel at this volumetric flow rate.This means that it takes about 2 s for the flow to get established and this clearly happens as the detection region is located inside the magnet at a distance of 52 cm from the bottom magnetic plate.In the porous medium, a spin particle displaces approximately 0.2 cm at Q = 1.4 cm 3 /s during = 500 ms.This distance is less than the average pore size d av p = 0.35 cm.To observe temporal flow fluctuations in the porous medium a particle should be allowed to percolate distances of the order of at least two or three d p ; therefore, profiles in the porous medium should not be different within the chosen flow encoding times.This assumption is clearly confirmed by the velocity distributions in the porous medium displayed in Fig. 9(a).
Temporal averaging in the presence of transverse permeability discontinuity results in a slightly different picture.In the clear channel measured peak velocities are higher for = 500 ms than those for = 8 ms.Nonetheless, these differences are within 10% of each other and might simply be explained by instabilities in the flow rate that become more observable at = 500 ms.The velocity FIG. 10.Azimuthally averaged radial velocity profiles (a) in the presence of transverse permeability discontinuity and (b) in the absence of the discontinuity in permeability as a function of the adjusted flow rate.The velocity fields were measured using the stimulated-echo protocol, with a raw acquisition matrix size of 128 × 128, a FOV of 3 × 3 cm 2 , and a slice thickness of 0.5 cm at = 500 ms.The measured slip velocities in (a) were 1.0 ± 0.3, 1.6 ± 0.5, and 2.5 ± 0.8 cm/s for flow rates of 1, 2.5, and 5 cm 3 /s, respectively.fields in the porous medium also differ within 10% of each other and are probably associated with similar instabilities in the flow rate.Moreover, temporal averaging for longer time results in a more homogeneous velocity profile and can be explained by averaging individual velocity fluctuations over a longer observational time.Despite these fluctuations, temporally averaged, rather than snapshot, velocity measures are more appropriate for these complex systems, as they allow for a better estimate of velocities in the boundary region where additional fluctuations in velocity values might be induced by, on the average, less than d av p in the porous medium, well away from the boundary region pore size.This decrease in d p in the boundary region is associated with the modification of the surface of the porous medium during the manufacturing of the clear channel.

C. Azimuthally averaged radial velocity profiles under varying mass transport
Velocity distributions were also imaged with and without transverse permeability discontinuity for a few flow rates as shown in Figs.10(a) and 10(b).The flow rate was varied to measure the effect of fluctuations in mass transport on slip velocity values at the interface and on the general character of flow in both cases.In both cases the uncertainty in velocity measurements of faster displacing species was increased with the increase of both the flow rate and the diffusion time .This most This also explains the distorted shape of the velocity profiles in the open channel in both cases with no reliable velocity measures with the increased mass transport.In addition, a greater error in Fig. 10(b) in the open channel is also due to the increased bin size in the radial averaging procedure as the bin size was decreased to 100 for the clarity of presentation.In the absence of the momentum transport across the interface [Fig.10(b)] the boundary velocities at the open channel side of the interface and at the porous medium side of the interface were zero and were independent of the flow rate, as expected when the no-slip boundary condition is imposed.No slip at either side of the interface was observed.
Transverse permeability discontinuity promotes fluid momentum transport across the interface that manifests in the apparent slip in velocities at the interface.Measurable nonzero fluid velocities [Fig.10(a)] at the interface were present at all flow rates studied.Slip velocities were also increased as a function of flow rate, as shown in Fig. 11.
The averaged Darcy velocities were also plotted as a function of the varying mass transport.In the case of an impermeable boundary Darcy velocities were increasing linearly with the increase in mass transport, as expected.In the case of a permeable interface the slip velocities were following a linear dependence while Darcy's velocities were best described by a quadratic dependence (including the zero point).The linear dependence in this case was attempted first and is shown as a red dashed line in Fig. 11.The quality of the fit improved significantly, assuming a quadratic dependence, shown as a black dashed line in Fig. 11.The presumed deviation from linear behavior in this case suggests the increase in the momentum flux between the clear channel and the porous medium with increased velocity.This observation might be in line with what would be expected due to increasing eddy production at the very rough interface created by the porous medium at the interface with the clear channel.8) at r < R, Eq. ( 9) at r > R (red), and LBM (blue).Velocity fields were temporally averaged within = 8 ms.The inset shows zoomed in boundary behavior of the exact solution to the problem and LBM modeling.Note a better match to the experiment for LBM.

D. Comparison of experimental velocity fields in the presence of the transverse permeability discontinuity to theoretical and LBM velocity fields
The velocity field at Q = 1.4 cm 3 /s was compared to the velocity fields constructed using Eqs.( 8) and ( 9) and LBM.These computational velocity fields were compared with experimental data as displayed in Fig. 12.Both approaches provided a flow description that coincided very well with the experiment in both the open channel and the porous medium.The major difference between the two models was in the description of flow behavior in the boundary region.
As can be seen from the inset in Fig. 12, the theoretical flow behavior derived using both Eqs. ( 8) and (9) shows an almost stepwise change in the velocities when the boundary changes to the open channel.This reflects the sensitivity of this combined approach to the choice of the boundary condition or to the averaging procedure as discussed in the Introduction.The velocity profile derived using LBM demonstrates a smooth transition in velocities from the boundary to the open channel, thus confirming the assumption that the boundary velocities can be evaluated using this method without any assumption of the boundary velocities.

V. CONCLUSION
MRI velocimetry methods were shown to be very efficient in quantifying the boundary behavior of fluids in the presence of transverse permeability discontinuity.The encouraging agreement between theory and experiment proves that MRI velocimetry methods can be used to extend the availability of experimental flow data in heterogeneous opaque porous media.The larger availability of such data may greatly contribute to the further theoretical development of these types of systems and help in refining the present theories and simulation modalities.As a prospect for future work in this area, MRI velocimetry methods should be used to probe flow in heterogeneous porous media of smaller pore size.If successful, these methods should be extended to studying flow of more complex rheology fluids in porous media.The establishment of MRI methodology to study generally complex flows in opaque systems could further contribute to benchmark experimental data in order to improve the fundamental understanding of many industrial (catalysis) and nature (geological) relevant heterogeneous flows, especially in the boundary region.

FIG. 1 .
FIG. 1. Diagrams and relevant definitions describing (a) one-dimensional fluid flow and (b) in flow plane cross section of a fully developed axially symmetric fluid flow in a porous medium with a model sample geometry.See the text for further explanations.

044102- 6 FIG. 2 .
FIG. 2. Timing diagrams of flow measurements using (a) spin-echo and (b) stimulated-echo MRI protocols.(c) Timing diagram for measuring bulk displacement spectra by NMR.Motion encoding gradients are applied in the direction of flow.Bars indicate that the time axes are not continuous in (b).
FIG.3.Displacement spectra (propagators) for water flow in a bulk reticulated foam at Re ≈ 13 for different times: (a) measured normalized propagators using NMR protocol displayed in Fig.2(c) (solid lines) and (b) measured (symbols) and simulated (dashed lines) with a bimodal Gaussian function and deconvoluted normalized propagator at = 1000 ms.

FIG. 4 .
FIG. 4. (a) Example of an individual 2DXμCT image of the foam with relevant length scales.(b) A 3D visualization of the cylindrical porous foam with a coaxial clear-cut circular channel reconstructed using a series of the segmented 2D XμCT images displayed in (a).The solid matrix is shown in yellow, while transparent designates pore spaces.The radius of the clear channel in this porous medium produced specifically for XμCT was approximately 0.5 cm.(c) Azimuthally averaged porosity profile produced from the reconstructed porous medium displayed in (a).

FIG. 6 .
FIG. 6.(a) and (c) Velocity maps and (b) and (d) corresponding azimuthally averaged radial velocity profiles (a) and (b) in the absence and (c) and (d) in the presence of the fluid momentum transport across the interface.The greater error in velocities in (b) is due to the larger bin size used in azimuthal averaging.The increase of errors in velocities near the wall in both experiments is induced by the distortions in image intensity due to the rf inhomogeneity in the coil as mentioned elsewhere in the paper.Larger errors in velocities in (a) in 0.85 < r/R < 1.14 are due to the absence of the collected NMR signal from the location of the rigid tube placed inside the channel.

FIG. 7 .
FIG. 7. (a) and (c) LBM simulated normalized velocity maps and (b) and (d) corresponding normalized azimuthally averaged radial velocity profiles (a) and (b) in the absence and (c) and (d) in the presence of the fluid momentum transport across the interface.The blue line in (b) represents the analytical solution of the velocity profile for Poiseuille flow.The irregular white shapes in (a) and (c) are the foam.

FIG. 8 .
FIG.8.Effect of spatial averaging on flow in the channel with a permeable wall at the volumetric flow rate Q = 2.25 cm 3 /s.Radial velocity profiles were produced by azimuthal averaging with a bin size of 400 from velocity maps produced by spin-echo flow experiments, acquired in 256 × 256 matrices for 12 different values of field gradients with = 8 ms and a FOV of 3 × 3 cm 2 .

FIG. 9 .
FIG. 9. Effect of temporal averaging on flow in the channel (a) decoupled from and (b) coupled to the flow in the porous medium at a volumetric flow rate Q = 1.4 cm 3 /s.The = 8 ms profiles were produced from the spin-echo experiments acquired in 256 × 256 matrices with 12 flow encoding gradients.The = 500 ms profiles were produced from the stimulated-echo experiments acquired in 128 × 128 matrices with 8 flow gradients.The FOV in both experiments was 3 × 3 cm 2 with a slice thickness of 1 cm.

FIG. 11 .
FIG.11.Slip and averaged Darcy velocities as a function of the increased mass transport in the presence of transverse permeability discontinuity (red) and in the absence of the discontinuity in permeability (blue).Note the quadratic Darcy velocity dependence on the increased mass transport (dashed black line) in the presence of transverse permeability discontinuity.The quadratic analysis was performed assuming that no slip velocities are observed at zero flow rate.Hence the zero point (black open triangle) has no error bars.
FIG.12.Comparison of radial velocity profile in the presence of fluid momentum transport across the permeable interface at Q = 1.4 cm3 /s (open symbols) with Eq. (8) at r < R, Eq. (9) at r > R (red), and LBM (blue).Velocity fields were temporally averaged within = 8 ms.The inset shows zoomed in boundary behavior of the exact solution to the problem and LBM modeling.Note a better match to the experiment for LBM.