Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow

Flow in porous media and the resultant hydrodynamics are important in fields including but not limited to the hydrology, chemical, medical and petroleum industries. The observation and understanding of the hydrodynamics in porous media are critical to the design and optimal utilization of porous media, such as those seen in trickle-bed reactors, medical filters, subsurface flows and carbon sequestration. Magnetic resonance (MR) provides for a non-invasive technique that can probe the hydrodynamics on pore and bulk scale lengths; many previous works have characterized fully saturated porous media, while rapid MR imaging (MRI) methods in particular have previously been applied to partially saturated flows. We present time- and ensemble-averaged MR measurements to observe the effects on a bead pack partially saturated with air under flowing water conditions. The 10 mm internal diameter bead pack was filled with 100 μm borosilicate glass beads. Air was injected into the bead pack as water flowed simultaneously through the sample at 25 ml h−1. The initial partially saturated state was characterized with MRI density maps, free induction decay (FID) experiments, propagators and velocity maps before the water flow rate was increased incrementally from 25 to 500 ml h−1. After the maximum flow rate of 500 ml h−1, the MRI density maps, FID experiments, propagators and velocity maps were repeated and compared to the data taken before the maximum flow rate. This work shows that a partially saturated single-phase flow has global flow dynamics that return to characteristic flow statistics once a steady-state high flow rate has been reached. This high flow rate pushed out a significant amount of the air in the bead pack and caused the return of a preferential flow pattern. Velocity maps indicated that local flow statistics were not the same for the before and after blow out conditions. It has been suggested and shown previously that a flow pattern can return to similar statistics if the preceding flow history is similar.


Introduction
A porous medium is characterized as a structural material composed of void spaces. In many situations, the pores exist in a connected network. In nature, porous materials can harbor water and organisms or offer a means to a lightweight yet rigid structure. Porous media can be utilized in industry for chemical reactions, filtering processes, chromatography and many other uses. The physics of liquid and gas infiltration in porous media is of relevance to biological, geophysical and industrial systems. In many systems of interest, cyclic saturation and partial saturation events control the transport processes in the porous structure. A particular example is the seasonal hydrologic cycling of contaminant radionuclide transport in the earth's subsurface [1]- [3]. An intricate interplay of surface tension, advection, diffusion and buoyancy forces renders the quantitative modeling of partially saturated flow difficult. Cyclic infiltration of liquid into partially saturated porous media results in a redistribution of flow paths due to a redistribution of clusters of saturated pores that form the backbone of advective transport [4,5].
Nuclear magnetic resonance (NMR) has proven effective in determining transport phenomena of porous media [6]- [9]. It has been used to characterize hydrodynamic dispersion, fluid phase interactions and other transport phenomena in porous media for both saturated [10]- [15] and two-phase flows [16,17]. Saturated single-phase flows of water in a bead pack have been well studied, but partially saturated flows have been less studied and are not thoroughly understood. NMR methods have been applied to determine the dynamics of combined liquid and gas flows in porous media [18,19]. Rapid NMR velocimetry methods [18,20,21] provide details of the complex temporal variation in dynamics within the porous structure on a timescale of 100 ms. Two-phase flow has been shown to exhibit hysteresis between imbibition and drainage 3 by using MR imaging (MRI) [18], by recording the pressure drop versus flow rate [22] and by x-ray radiography [23]. MRI experiments [10,18,24] have shown that single phase flow does not have hysteresis if the bead pack is fully saturated. However, if air is present in the pack, then hysteresis effects in regard to pressure drops across the porous media can occur if the beads are randomly packed and there is not a great deal of order [24]. This work shows that a partially saturated single-phase flow has global flow dynamics that return to characteristic flow statistics once a steady-state high flow rate has been reached. This high flow rate pushed out a significant amount of the air in the bead pack and caused the return of a preferential flow pattern. Velocity maps indicated that local flow statistics were not the same for the before and after blow out conditions. It has been suggested and shown previously that a flow pattern can return to similar statistics if the preceding flow history is similar [25].
In this paper, we analyze the stationary dynamics of the transport during spatially and temporally variable processes. Of particular interest is the observation that the ensemble averaged system dynamics are structurally (i.e. pore size) dependent and globally reproducible over a fixed range of displacement scales despite local spatial and temporal variations after cycling of gas and liquid flow rates. A particular model of partially saturated flow in porous media that is consistent with this result is a percolation network. In percolation models, the transport or flux is independent of the exact spatial location of connected pores; rather it depends on the probabilistic percentage of occupied pores [26]. While not a direct quantitative confirmation of this idea, our data presented here indicate that this conceptual model is correct. The immiscible displacements of gas and liquid phases in porous media have provided a rich example for physics theories of percolation [27,28], and fractal interfacial structure and dynamics [29]- [32] scaling [26]. Fluid invasion flows are most studied in the context of these theories and are transient imbibition and drainage-type flows [26]. Steady-state flows have received far less experimental and theoretical attention [33]- [35]. These steady-state flows are of great interest since they represent the stationary stochastic states of a non-equilibrium system with external energy input balancing irreversible energy dissipation processes [33]- [35]. This makes this system of relevance to understanding non-equilibrium thermodynamics [4,36]. In this paper, we present NMR measurements of the liquid-and gas-phase fractions and stationary dynamics for the steady-state flow of water over a range of volumetric rates through a homogenous porous medium of packed spheres in a partially saturated state.

Theory
The modeling of imbibition and drainage of fluids in porous media has led to the development of the theory of invasion percolation [27,37]. The invading fluid generates clusters and structures that exhibit power-law size scaling [28]. Much of the theoretical development has focused on transient fluid invasion processes and is reviewed in chapter 12 of Sahimi's book [26]. The balance of capillary and viscous forces, described by the capillary number [38] Ca = µv w /σ , provides a characterization of the flow in terms of the physical properties of viscosity µ, surface tension σ and velocity v w of the wetting fluid. Transient imbibitions and drainage studies have been conducted for liquid-liquid systems, such as oil and water [26,28,31], and gas-liquid systems, such as air and water [29,30,32]. A primary result of invasion percolation studies is the scaling of the invading fluid saturation S ∝ L (D f −2) , where L is the lattice size and D f is the fractal dimension. The subsequent scaling of cluster size is N ∝ L D f , where D f ∼ 1.89 for two-dimensional (2D) models and D f ∼ 2.52 for 3D models [27,28]. Other transient models deal with the scaling of interface thickness with capillary number [29] and power-law noise scaling of interface fluctuations [30].
The less studied steady-state flow regime represents a continuous interaction between imbibitions and drainage [39]. In this case, cluster size distributions exhibit a power-law probability distribution with an exponential cutoff p(l) ∼ l −2 exp(−l/l * ), where l is the cluster size and l * is the cutoff cluster size [33,34,39]. The cutoff cluster length has been shown to scale as a power law with Ca as l * ∼ Ca −ζ with ζ ∼ 0.98 in 2D [34]. The non-wetting fluid saturation has been observed to decrease with increasing Ca and a theoretical explanation is lacking [34]. The dynamics in these experimental studies have been characterized in terms of the dissipation of energy by the fluctuating interfacial area between the wetting and non-wetting phases at the steady state [33]. This analysis leads to the scaling of pressure drop in the bead pack with capillary number as |∇ P| ∝ Ca β with β ∼ 0.5 in 2D experiments [33]. This result is consistent with a physical model in which the wetting fluid flows subject to an effective permeability which varies as κ eff ∼ √ Ca assuming D'Arcy's law v ω ∼ κ eff /µ ω |∇ P| holds [33]. The Ca scaling results are related to the idea that the cutoff cluster size l * represents an extension length for cluster mobility [33]. This is important since in the work of Tallakstad et al [34] it is assumed that the energy dissipation is localized in a dissipative wetting fluid volume V d = d V total l * , where d is the grain size of the porous media and l * , the cutoff length, gives the spacing between the flowing channels of wetting fluid in channels of width d [33]. In the discussion on the dynamics measured by PGSE NMR, we will examine the role of the pore size in this context.
While the scaling theory of cluster size distribution and mean-field pressure drop has recently been studied experimentally in steady-state partially saturated flows, hydrodynamic dispersion has not been similarly analyzed. de Gennes [4] has presented the detailed theory of hydrodynamic dispersion, the interaction of random diffusive dynamics and coherent advection [40], in partially saturated porous media. The wetting fluid is treated as a percolation cluster with a backbone and dead ends [4]. The theory provides a detailed discussion of the scaling of the dispersion during external steady pressure-driven D'Arcy flow in terms of a critical pressure P c to form an infinite cluster and hence generate fluid flow through the porous medium. The scaling is found to depend on the average fluid velocity of all fluid particles v = v s /s w , where v s is the superficial velocity Q w /A and s w is the wetting fluid saturation level. The dependence onv rather than the backbone average velocityv B = v s /s B , where s B is the probability of a fluid particle being on the backbone and is an important feature of this model of partially saturated flow [4]. As a result, de Gennes [4] finds the longitudinal dispersion scales is the ant in a labyrinth diffusivity, dependent on P = P − P C , the molecular diffusivity D 0 with scaling exponentst− ∈∼ 1.31 and ξ is the correlation length or 'mesh size of a cluster'. The correlation length scales as ξ ∝ d P with d being the pore size diameter and v ∼ 0.9 [4]. The fact that the scaling of the dispersion D is with the average fluid velocityv is important since, in NMR measurements of the propagator of motion, all fluid particles contribute to the measurement. In this model, the correlation length ξ is associated with the maximum size of the dead ends. As a side note, it is often missed in the literature that this paper by de Gennes [4] is one of the first to apply non-equilibrium statistical mechanics in the context of the power-law waiting time distributions of Scher and Montroll [41] and memory functions to hydrodynamic dispersion. This presages the large body of work on the application of non-equilibrium statistical mechanics [42] to porous media transport in the context of continuous time random walks [43] over the last 20 years.

Nuclear magnetic resonance (NMR) experimental techniques
NMR relies on the relationship between the Larmor precession frequency ω o of active nuclei magnetic moments, such as for the proton 1 H, and the large external magnetic field B 0 , ω 0 = γ B 0 , where γ is the gyromagnetic ratio. The additional use of magnetic field gradients G results in a resonant frequency that depends on the position r, ω eff = γ B 0 + G · r [44], and it is possible to encode the NMR signal for spatial location and displacements.

Free induction decay (FID) experiments
The simplest possible NMR experiment involves excitation of the entire sample with a radiofrequency (rf) pulse applied at the Larmor frequency. Acquisition of the resulting signal provides a measure of the total water signal and provides a method for determining relative water saturation levels or non-wetting saturation levels S nw .

Magnetic resonance imaging (MRI) density maps
MRI density maps provide water density distributions in two dimensions within a slice that is excited non-invasively in the third spatial dimension through the timed application of rf excitation pulses and gradient pulses. The exact timing of the gradients in three orthogonal directions (G read , G phase , G slice ) for the MR imaging density maps presented in this paper is shown in the schematic in figure 1(a). During acquisition, the NMR signal can be written as S(k) = ∫∫ ρ(r) exp(i2π k · r)dr, where the spatial wave vector is k = (2π) −1 γ (t read G read +t phase G phase ). Fourier transform of the 2D k space data yields the spin density, ρ(r), averaged over the slice thickness and pixel dimension, which, in the absence of significant NMR relaxation weighting, is the desired MR density map.

Propagators (pulsed gradient spin echo (PGSE) experiments)
The propagator experiment excites a slice in the sample and can provide the distribution of displacements of water molecules within this slice over an observation time of interest within the range of 1-1000 ms. The measurement of translational motion, whether it is due to self-diffusion, coherent velocity or hydrodynamic dispersion, is achieved by combining the resonant frequency dependence on spatial magnetic field variation with the strategic application of a pair of rectangular narrow gradient pulses of amplitude g and duration δ separated by an observation time , termed the pulsed gradient spin-echo (PGSE) sequence (see figure 1(b)) [44]. This sequence encodes the acquired spin-echo signal for displacements that occurred in the observation time and hence accesses averaged and timescale-dependent transport information [45]. In the PGSE experiment, assuming that the gradient g is applied for a sufficiently short time relative to the characteristic time of the process of interest, the data can be analyzed using the propagator formalism. The propagator, or Van Hove selfcorrelation function, is the conditional probability that a molecule starting at position r will move to the position r + R over a time period , while the average propagatorP(R, ) that we present in this paper is defined as the average over the initial molecular distribution ρ(r). The acquired NMR signal can be written as S(q) = ∫∫ ρ(r)P(r, R, ) exp(i2π q · R)drdR = ∫∫p(R, ) exp(i2π q · R)dR, where the displacement wave vector is q=(2π) −1 γ δ g. Fourier transform of the q space data yields the averaged propagator of molecular displacement R over observation time ,P(R, ).

Velocity maps
The MRI density sequence and the PGSE sequence can be combined to allow noninvasive measurement of the velocity map and hence spatial resolution of transport processes (see figure 1(c)) [44].

Methods-general
The NMR experiments were conducted on a Bruker super-wide bore 300 MHz magnet with a Micro2.5 magnetic field gradient probe containing a 10 mm birdcage rf coil and networked to a Bruker Avance III spectrometer with gradient controls. The gradient set is capable of producing orthogonal magnetic field gradients up to 1.5 T m −1 in three orthogonal directions. All experiments were performed on the same porous media bead pack sample. NMR and MRI techniques are inherently non-invasive and therefore the sample did not have to be disturbed in any way to obtain the experimental results presented here.
The 50 mm long porous medium was constructed from 100 µm borosilicate glass beads in a 10 mm ID diameter glass tube. The beads were wet filled with deionized water and compacted in the tube to achieve a porosity of approximately 44%. The partially saturated state was created by setting the water flow rate to 25 ml h −1 and varying the simultaneous flow rate of the air. Two Pharmacia P-500 HPLC pumps were used to deliver the air and deionized water phases against gravity, which combined at a tee junction just below the entrance of the bead pack. This flow resulted in an alternating flow of water and then air bubbles up through the bead pack. The dual air/water flow was continued for 8 h. The air flow was then shut off and the water flow continued at 25 ml h −1 until a steady-state NMR signal indicated that a steady state of partial saturation had been achieved (approximately 4 h). The volumetric rate of the initial air flow was varied between 25 and 100 ml h −1 , but there was no difference observed in the final partially saturated state as observed by the NMR results.
After the partially saturated state had been established, MR experiments were performed as the water flow was increased from 0 to 500 ml h −1 and then decreased to 25 ml h −1 to repeat the velocity map before increasing to 50 ml h −1 and eventually shutting off flow. In table 1, flow rates are listed in the order they were used and the table indicates the experiments that were run at each flow rate. The intention was to monitor the degree of partial saturation as the flow rate was increased from 0 to 500 ml h −1 and to characterize the bead pack before and after the fast flow rate challenge with a density map, propagator and velocity map. The hold time for each flow rate was the time required for completion of all NMR experiments at each respective flow rate. At flow rates 25 and 50 ml h −1 , this was about 20 min, and at all other flow rates, it was 5 min.
For the MRI density maps, which indicate the presence of water (white) and the presence of air or bead (black), at 0 ml h −1 a resolution of 78.1 × 43.0 µm 2 and a slice thickness of 1.5 mm were used. The slice was excited at the center of the bead pack. The echo and repetition times were T E = 16.42 and T R = 1000 ms, respectively, and four averages resulted in an experiment time of about 17 min for a single data set. For the propagator measurements at 25 ml h −1 , a slice thickness of 10 mm was used. The slice was excited at the center of the bead pack. The echo and repetition times were T E = 5 ms and T R = 1000 ms, respectively, and eight averages resulted in an experiment time of about 17 min for a single data set. The displacement observation time was = 400 ms, the pulse gradient length was δ = 1 ms and 128 gradient steps were used with a gradient strength increment of 18 mT mm −1 , resulting in a min/max gradient strength of ±1.134 T m −1 . For the velocity maps measured at 50 ml h −1 , a resolution of 156.2 × 156.2 µm and a slice thickness of 10 mm were used. The slice was excited at the center of the bead pack. The echo and repetition times were T E = 13 ms and T R = 2000 ms, respectively, and four averages resulted in an experiment time of 35 min for a single data set. The displacement observation time was = 9 ms, the pulse gradient length was δ = 1025 ms, and the gradient pulses were g = ±1.06 T m −1 . The free induction decays (FIDs) measured the total water signal from excitation of the entire sample at each flow rate from 0 to 500 ml h −1 , thus providing a method for determining relative water saturation levels or non-wetting saturation levels. For these saturation measurements, a repetition time of 5 s and 30 averages resulted in an experiment time of 35 min for a single data set.

Results and discussion
The MRI density maps (figure 2) and FID signals (figure 3) measured under 0 ml h −1 no flow conditions before and after the maximum flow rate indicated that larger pockets of air were present before but not after the maximum flow rate. For all repetitions, initial partially saturated states show local regions in the MRI density maps that have low or no water signal,  indicating larger multiple pore size air pockets. The MRI density maps of the same region after the maximum flow rate show a return of the water signal intensity throughout the sample and indicate that the larger air pockets are eliminated by the fast flowing water. Although the MRI density maps show no evidence of air after the maximum water flow rate of 500 ml h −1 (figure 2), the slightly lower FID normalized area (figure 3) indicates that smaller air pockets below the resolution of the MRI density maps are present. The FID measurements were repeated at regular intervals as the flow rate was increased from 0 to 500 ml h −1 (figure 3) in time steps of 3-5 min. The FID signal begins to increase at 250 ml h −1 , indicating the flow rate at which air begins to leave the bead pack. This corresponds to a capillary number of 7.6 × 10 −4 . The data in figure 3(a) have been plotted against the capillary number, Ca ( figure 3(b)), to be presented in a manner consistent with data previously published for a 2D system [34] to allow for direct comparison. For the conversion of the flow rate to capillary number, the definition given in the theory section was used [38]. The non-wetting saturation level S nw is calculated using the ratio of the FID normalized area FID NA of the partially saturated condition, FID NA(ps) , to the fully saturated condition, FID NA(fs) , as S nw = 1 − FID NA (ps) FID NA (fs) . The inverse proportion of S nw to Ca has been discussed in the context of compressibility of the non-wetting gas phase (air) in analysis 10

FS
I II III IV Figure 4. Velocity maps from a 10 mm slice excited at the center of the 10 mm ID bead pack. Top row: the velocity maps for the initial partially saturated state. The signal intensity is proportional to the average water velocity in the volume (voxel) centered around the pixel of the image; bottom row: the velocity maps after a 500 ml h −1 water flow has pushed out much of the air. Column 1 shows the initial fully saturated sample and the remaining columns show data from four repetitions I-IV of the initial partially saturated state and then the saturation after the high flow rate. The white regions have been falsely colored to indicate when the NMR signal was below the noise level due to regions of higher dispersion.
of the 2D data of Tallakstad et al [34], and they conclude that the compressibility of the air is not the cause of the decrease in S nw with Ca. Tallakstad et al also note the dearth of data on saturation in steady state. Computer simulations comparing saturation levels and pressure drop in steady partially saturated flows have been undertaken [46] but no theory for S nw scaling with Ca exists, and, to our knowledge, figure 3(b) is the first to report of this in 3D [39]. The velocity maps measured at 50 ml h −1 before and after the maximum flow rate of 500 ml h −1 are shown in figure 4. These velocity maps provide further evidence that some pores below the resolution of the MR image are filled with air as there are significant regions with no flow and other regions with high backbone flows. Such velocity structures are not visible in the initial homogeneous bead pack (figure 4 (column 1)). However, they are present in all of the subsequent partially saturated states and most of the states after the maximum flow rate has blown out the large air pockets. So no flow and high back bone flow structures are expected in a percolation model structure [4]. The variation in the number of spatial regions of no signal, i.e. regions of higher dispersion caused by small air pockets, in the velocity maps for each trial after a maximum flow challenge are consistent with the slight saturation variations in figure 3.
Pressure drops across the sample were measured as the flow rate and hence capillary number were increased, as shown in figure 5. A power-law scaling of approximately 0.33 was observed, which is consistent with the percolation cluster model that predicts | P| ∝ Ca γ constrained by γ 0.5 [34]. In the limiting case of identical viscosities for the two fluids, a scaling of 0.5 is expected and this represents an upper limit to the expected scaling. The propagators measured at 25 ml h −1 before and after the maximum flow rate of 500 ml h −1 are shown in figure 6. The propagator for the initial fully saturated flow is shown in both figures 6(a) and (b). This fully saturated homogeneous bead pack propagator has two peaks, one centered around zero displacement representing water that is initially in the no flow regions of the bead pack and the other centered around the average water velocity, with the dip in the propagator corresponding to the pore size d = d p φ/(1 − φ) and hence nominal bead diameter d p , as has previously been observed and discussed in many detailed MR studies of fully saturated flow in homogeneous bead packs [12]- [14], [16]. The propagators for the initial partially saturated states before maximum flow ( figure 6(a)  resulting from the many large air blockages to the water flow path with regions where the water is trapped and pathways where the water can flow faster [47]. All of the propagators taken after the maximum flow rate of 500 ml h −1 (figure 6(b)) look qualitatively similar to each other and to that obtained for the initial fully saturated state. Although FID measurements ( figure 3) show that there is still some air present, and velocity maps ( figure 4) show that the velocity distribution is no longer spatially homogeneous, these propagator results clearly show that the global flow statistics have returned, or nearly returned, to the fully saturated flowing condition. This also suggests that the flow returned to some characteristic flow pattern that closely resembled the dynamics of the saturated bead pack. This pattern was then stable and could not be altered unless air was introduced again.
For clarity and ease of comparison, the propagators from each repetition are compared individually to the propagators for the saturated condition (figure 7). As expected, the propagators from before and after maximum flow rate for the initial fully saturated condition show identical inflections at the pore size d and similar probabilities of displacements. This corroborated the results of the initial MR image density maps and FID signals that also showed limited difference after the high flow rate. For the initial partially saturated states, it is clear that the flow dynamics converge to the saturated dynamics after the maximum flow rate of water has gone through and removed the large air pockets, even though the FID signal and velocity maps show clearly that more air than in the saturated system is still present. The higher probability of displacement around d ∼ 79 µm indicates the impact that the underlying pore structure has on the dynamics in the bead pack, consistent with the cutoff length l * representing a cluster mobility extension length [33]. That is, there is an increased probability of displacements at the pore scale d, which then decays away until the cutoff length scale is reached and backbone, long displacement, dynamics dominate. The theory of de Gennes [4] indicates that the mean of the entire distribution of dynamics, including the no flow peak around zero displacement, is the appropriate mean displacement for use in analyzing dispersion dynamics. While the propagators for the fully saturated and partially saturated after maximum flows have many similarities, one aspect that requires further study is that on a log scale the partially saturated after maximum flow propagators consistently have a long displacement tail, i.e. higher probability of large displacements.

Conclusions
In this paper, data on the stationary dynamics and saturation of steady-state partially saturated flow in porous media before and after flushing with a high water flow rate to mimic cyclic invasion phenomena have been presented. Data on the saturation of porous media during partially saturated flow are limited and the unique ability of NMR to provide these data is clearly demonstrated. The ability of NMR to measure the propagator of motion yields the stationary dynamics, and connections with the structure of the porous media and percolation models of the flow can be made. In demonstrating these measurements and interpretations, the basis for further studies of steady-state partially saturated flows for a range of Ca numbers, fluids, saturation levels and NMR parameters, such as displacement observation time, has been established.