Article pubs.acs.org/est
Critical Role of the Immobile Zone in Non-Fickian Two-Phase Transport: A New Paradigm Nikolaos K. Karadimitriou, Vahid Joekar-Niasar,* Masoud Babaei, and Craig A. Shore School of Chemical Engineering and Analytical Science, Faculty of Engineering and Physical Science, The University of Manchester, Manchester M13 9PL, United Kingdom S Supporting Information *
ABSTRACT: Using a visualization setup, we characterized the solute transport in a micromodel filled with two fluid phases using direct, real-time imaging. By processing the time series of images of solute transport (dispersion) in a two fluid-phase filled micromodel, we directly delineated the change of transport hydrodynamics as a result of fluid-phase occupancy. We found that, in the water saturation range of 0.6−0.8, the macroscopic dispersion coefficient reaches its maximum value and the coefficient was 1 order of magnitude larger than that in single fluid-phase flow in the same micromodel. The experimental results indicate that this non-monotonic, non-Fickian transport is saturation- and flow-rate-dependent. Using real-time visualization of the resident concentration (averaged concentration over a representative elementary volume of the pore network), we directly estimated the hydrodynamically stagnant (immobile) zones and the mass transfer between mobile and immobile zones. We identified (a) the nonlinear contribution of the immobile zones to the non-Fickian transport under transient transport conditions and (b) the non-monotonic fate of immobile zones with respect to saturation under single and two fluid-phase conditions in a micromodel. These two findings highlight the serious flaws in the assumptions of the conventional mobile− immobile model (MIM), which is commonly used to characterize the transport under two fluid-phase conditions.
■
INTRODUCTION Pore-level characterization of dispersion in porous media has recently received significant attention, because high-resolution fast imaging allows for profound enhancement of our understanding of pore-level processes.1,2 Despite such significant advancements in solute transport in single fluid-phase porous media, our understanding of solute transport at the pore level under two fluid-phase conditions is still in its infancy. Pore-scale analysis and characterization of this challenge will contribute to our understanding of natural and environmental phenomena2−4 and improve our predictive modeling capabilities for engineering applications, such as colloid transport in unsaturated zones and soil remediation,5 enhanced oil recovery,6 and unsaturated soil physics.7 Dispersion in porous media, controlled by advection and diffusion, has been intensively studied in hydrogeology, where usually only one fluid phase exists.8−10 However, the characterization of transport under single fluid-phase conditions is not valid when an additional fluid phase exists.11−16 Although the majority of the literature deals with transport under saturated conditions, in several experimental studies, dispersivity was investigated under unsaturated (two fluid-phase) conditions.17−32 Hydrodynamic dispersion employs Fick’s law to relate dispersive flux to the concentration gradient, with C being © 2016 American Chemical Society
the concentration of the non-reactive solute and x being the longitudinal spatial coordinate in a one-dimensional system. ∂C ∂ 2C ∂C =D 2 −u ∂t ∂x ∂x
(1)
Equation 1 is usually referred to as the advection−dispersion equation, where u is the pore velocity and D is the longitudinal hydrodynamic dispersion coefficient in the x direction. For the continuous injection of a solution with an initial solute concentration of C0 [C(0, t) = C0, t ≥ 0] in an initially solute-free semi-infinite porous medium [C(x, 0) = 0, x > 0] and zero concentration gradient at the outlet boundary [∂C(∞, 0)/∂x = 0, t ≥ 0], the analytical solution of eq 1 reads C(x , t ) =
⎛ x + ut ⎞⎤ 1 ⎡ ⎛ x − ut ⎞ ⎟ + e xu / D erfc⎜ ⎟⎥ C0⎢erfc⎜ ⎝ 2 Dt ⎠⎦ 2 ⎣ ⎝ 2 Dt ⎠ (2)
Dispersion under two fluid-phase conditions is postulated to be a function of flow velocity and saturation. The proposed functional form is analogous to the empirical relation between Received: Revised: Accepted: Published: 4384
December March 24, March 24, March 24,
3, 2015 2016 2016 2016 DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology dispersion and dispersivity under single fluid phase,33 and the only difference is that the parameters are proposed to be functions of saturation as follows:17,34,35 D(S , u) = De(S) + λ(S)un
portion of the water saturation is denoted by Sim, and the mobile portion of the water saturation is denoted by Sm, so that Sim + Sm = S. Furthermore, Dm is hydrodynamic dispersion over the mobile water fraction, so that DmSm = DhS. Similarly, um is the mobile pore velocity, so that umSm + uS. In eq 5, α (T−1) is a first-order mass transfer coefficient, accounting for diffusion of solutes into the immobile, diffusion-dominated region. Equation 5 describes a non-equilibrium mass transfer between the two regions, mobile and immobile. Equation 6 relates Cm (solute concentration of the mobile region) and Cim (solute concentration of the immobile region) to C (solute concentration of the water-filled zone). The correction into the advection−dispersion equation by MIM is also reflected in the dispersion coefficient36
(3)
where S is the saturation, De is the effective diffusion coefficient of the medium, and λ and n are empirically determined coefficients. Note that this relation has no theoretical background. In various studies, the use of this relation has led to the overestimation of the dispersion coefficient for a given saturation.23,27,34,36,37 This mismatch was mostly attributed to the stagnant fraction of saturation. One can reasonably criticize eq 3, because for a given saturation, configurations of many fluids with different tortuosity conditions are possible. Thus, saturation alone may not be sufficient to capture the physics of the transport under two fluid-phase conditions. A universal, non-empirical functional form of dependency of dispersion upon two fluid-phase conditions, if any, has not yet been established. Against this backdrop, new models have been developed to define and establish a relationship between dispersive flux and concentration gradient for a not fully developed Fickian regime. Of such models are the dead-end pore model38 and the mobile− immobile model (MIM).36,39−41 In the MIM, a non-equilibrium mass transfer between the two regions, mobile and immobile, is defined. This is the major difference between the MIM and the dead-end pore model, where the system is allowed to reach equilibrium. This physically consistent assumption makes the MIM model more suitable to describe non-Fickian dispersion, and this is why we chose to investigate this model instead of the dead-end pore model. Until now, the MIM model has been used as a tool to fit the experimental data; however, the validity of the model and consistency of its predefined parameters and their functionalities have never been investigated. To gain insights into the physical processes of transport in two fluid-phase conditions in porous media, we will employ a transparent micromodel experimental setup that allows for realtime imaging of the transport of a solute in water in the presence of another immiscible fluid phase. By image processing, we will be able to (a) identify the mobile and immobile zones, (b) characterize dispersion, and (c) estimate the mass transfer rate between mobile and immobile zones under different flow rates at different water saturations. As such, our analysis is in contrast to the previous studies, which have used fitting the MIM to the experimental data.
DMIM = βDm + (1 − β)2
■
MIM In MIM, a clear distinction between two regions, i.e., flow contributing (mobile) region and stagnant (immobile) region, is made. The solute transfer between these two regions is modeled by a diffusive exchange rate.42 The full set of equations describing MIM for non-reactive transport is given here ∂Cm ∂C ∂ 2C ∂C + Sim im = SmDm 2m − Smum m ∂x ∂t ∂t ∂x
Simφ
∂C im = α(Cm − C im) ∂t
SC = SmCm + SimC im
(7)
where the ratio of mobile to the total saturation is denoted by β, (β = Sm/S). The first term represents the true dispersion in the mobile region (Dm), and the second term is the apparent dispersion as a result of solute exchange between the mobile and immobile regions. The contribution in the second term is only significant when transversal diffusion is dominating the mechanical dispersion. This correction has improved the predictive capability of the macroscale modeling of solute transport under two-phase conditions.23,25,31,43−46 Parameter β has been reported to increase with the increase of water saturation;23,35,36,39,42 however, Toride et al. showed that there is a parabolic relation between β and saturation, with a minimum at intermediate saturations.31 Among the four fitting parameters defined in the MIM, namely, β, α, Dm, and um, the first two parameters capture the interaction between mobile and immobile zones and dictate the non-Fickian behavior of the transport mechanism. These parameters have been only obtained by curve fitting, and no direct physical observation of these parameters has ever been reported.36,47,48 Here, on the basis of 15 sets of measurements for 5 different flow rates for water and 3 different injection rates for ink in a long micromodel, we have been able to determine a qualitative relation between the dispersion coefficient and saturation and pore velocity as well as the values for β and a functional relation between α and time. We have also been able to speculate the effect of the topology formed under two fluid-phase conditions and different flow rates on the formation of mobile and immobile regions.
■
Sm
um 2 αϕS
MATERIALS AND METHODS Experimental Setup. To evaluate the immobile and mobile zones as well as the mass transfer rate between them under different hydrodynamic conditions, we have visualized the real time spreading of ink in water in a polydimethylsiloxane (PDMS) micromodel, in the presence of Fluorinert (FC-43). Fluorinert is a nontoxic, nonflammable, silicon-based fluid, which served as the wetting phase in our experiments, and it is immiscible with water. By design, the topology of the flow network is fully percolating. More information on each component of the experimental setup and the experimental procedure is given inFigure S1 of the Supporting Information. In the experiments, we drained the initially Fluorinert-filled flow network with water, at different flow rates, to establish a range of water saturation from 0.2 to 1 with different
(4)
(5) (6)
where φ refers to porosity and the subscripts m and im refer to the mobile and immobile regions, respectively. The immobile 4385
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology
Figure 1. Processed image of a snapshot of the whole micromodel shown in Figure S2 of the Supporting Information. The micromodel is filled with Fluorinert and water, and ink is injected from the left boundary. Gray- and black-colored areas correspond to the Fluorinert-filled and solid areas, respectively. The areas where ink is present in water are shown with a spectrum of colors between deep blue (zero ink concentration) and red (highest ink concentration).
permeability, and entry capillary pressure do not change with the increase of the averaging window, which turns out to be 5 × 5 mm2. Knowing the REV size, the average resident concentration versus time was calculated for every REV per experiment. For a better illustration, we assigned a different color coding to each phase. More specifically, the gray-, blue-, and blackcolored areas correspond to the Fluorinert-filled, water-filled, and solid areas, respectively. The areas where ink is present in water are shown with a spectrum of colors between blue (zero ink concentration) and red (highest ink concentration). All of the intermediate concentrations are shown as the corresponding colors within the visible spectrum. Figure 1 shows an example of a processed image of the original greyscale image shown in Figure S2 of the Supporting Information. Differentiating Mobile from Immobile Regions. The immobile areas in the area occupied by water in the flow network are practically either rate-limited or totally inaccessible to ink via advection and accessible only via diffusion (they are hydrodynamically stagnant and referred to as dead ends or immobile zones). Identification of such areas by image processing is challenging because the analysis should be performed by applying a global metric considering the variation of the local concentrations over time. To discriminate mobile and immobile regions, fully water-saturated experiments play a key role. By design, there are no dead-end pores, and therefore, the global metric should result in no immobile zone in the fully water-saturated cases. In the first step of our analysis, temporal variation of the concentration was calculated by summation of the variation of pixel-based concentrations with respect to time over the mask nθ region, defined as Ċ (t) = 1/nθ∑(i,j) ∈ θ(dci,j(t)/dt). In our fully water-saturated experiments, where only water exists in the flow network, it is known that resident and outflow concentrations would reach the inlet boundary value because there are no immobile zones. Given this fact and by plotting Ċ (t) versus time, the threshold Ċ (t) was determined such that no immobile saturation in this case could result. Note that Ċ (t) could never reach absolute zero, because there was noise in the images as a result of the invisible to naked eye oscillations of the illumination. The threshold Ċ (t) value was obtained from the fully saturated micromodel experiments, and it was set to be around 10−3 depending upon each case. Using this threshold, the two-phase experiments have been analyzed. In the second step, the temporal variation of the pixel-based concentration was summed up over time as ċi,j = ∑nN= 1dci,j(tn)/dtn. If ċi,j < Ċ (t), that pixel belongs to the immobile zone, and if ċi,j ≥ Ċ (t), the pixel belongs to the mobile zone.
topological configurations. As a result, some parts of the waterfilled area would become hydraulically dead-ends. As soon as steady-state saturation was established, water-based ink was introduced into the flow network as the dispersed phase. Because ink was water-based, it was only soluble in water. Using image processing, we have explored the role of immobile zones in the two-phase non-Fickian transport. Image Processing and Data Analysis. Saturation. As the first step, a “mask” of the steady-state distribution of fluid phases in the flow network was created. This allowed us to calculate saturation and to isolate image processing in the water-filled areas where dispersion of ink occurred. By assignment of different grayscale intensity values to the pixels belonging to the water-filled area and Fluorinert-filled area, these two phases were distinguished and saturation was estimated. Resident Concentration. All of the images were acquired by a set of four monochrome cameras, as described in the Supporting Information. The level of light intensity in such cameras is expressed via grayscale intensity. This scale expands from 1, which corresponds to a completely black pixel, to 255, which corresponds to a white pixel. All images were corrected for the background illumination and cross-correlated to the mask. The average grayscale pixel intensity of the water-filled area of each image was then calculated. As ink was gradually transported into the water-filled area and had a lower grayscale intensity, the average grayscale intensity decreased. The transport continued until ink reached its maximum concentration, and the lowest gray scale intensity was attained. To establish the relation between grayscale intensity and concentration, we measured the grayscale intensity at different known concentrations during a separate experiment. The resulting nonlinear relation is discussed in Figure S3 of the Supporting Information. Denoting the highest grayscale intensity throughout the sequence of images as Imax, the lowest as Imin, and the grayscale intensity of the tth image for a pixel in the mask region with indices (i, j) as It(i,j), the normalized concentration of ink in the tth image was calculated b
θ −(aIi,j′ ) t t by normalization as C(t) = 1/nθ∑n(i,j) ∈ θc(i,j), where ci,j = e and Ii,j′ = (Imax − Ii,j)/(Imax/Imin), where nθ is the number of pixels located in the mask region, a and b are coefficients, and (i, j) is the pixel indices in the mask region. In this way, we were able to calculate the normalized concentration of solutes in every representative elementary volume (REV). An REV is defined as the volume of a homogeneous porous medium above which the system properties, such as porosity, permeability, relative permeability, and capillary pressure curves, are insensitive to the averaging domain size. In this study, the REV size has been selected such that porosity,
4386
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology
Figure 2. (Left) Dispersion coefficient estimated by fitting eq 2 to the average resident concentration versus water saturation for ink rates of 0.05 mL/h (■), 0.1 mL/h (▲), and 1 mL/h (●). The open symbols show dispersion under fully water saturated. (Right) Statistical variation of dispersion coefficients for each ink flow rate. The top and bottom of the box show 25th and 75th percentiles, respectively. Solid symbols show the mean values, and whiskers show the 1 and 99% percentiles.
■
For the highest (1 mL/h) and lowest (0.05 mL/h) ink flow rates, the non-monotonic trend of the dispersion coefficient versus water saturation is more obvious compared to the intermediate ink rate of 0.1 mL/h. Non-monotonicity in the dispersion−saturation relation implies that some hydrodynamic properties, such as immobile zones and flow network tortuosity, have a nonlinear relation with saturation. In a fully watersaturated system, all pores contribute to the transport, leading to a maximum resident concentration equal to the injection boundary value. However, with a decrease of saturation, only the water-filled part of the network can contribute to the transport. Under two-phase conditions, there will be some dead-end clusters that are hydrodynamically stagnant. In these zones, transport can only occur through diffusion in contrast to the flowing zones, where advective transport is dominant. This dual hydrodynamic characteristic of transport imposed by the degree of saturation and the two fluid morphology leads to a larger dispersion coefficient compared to single fluid-phase conditions. However, the variation of the dispersion coefficient versus saturation is not necessarily monotonic because the fraction of stagnant and advective zones is not constant for the range of saturations. There is significant data scattering, especially for the case of the ink rate of 0.1 mL/h. This data scattering is not due to the experimental or image processing error because (a) it is variable with the ink rate and (b) is absent in the case of fully watersaturated experiments and (c) the curve-fitting practice shows very small errors (quantitative analysis of the data is shown in Table S7 of the Supporting Information). That highlights the complex fundamentals of transport in two fluid-phase porous media that implies that dispersion may not be only a function of fluid saturation and other hydrodynamic properties, such as tortuosity, in a two fluid-phase system should be included. It is known that, for a given saturation, many configurations of fluids are possible. Each morphology will dictate a specific topology that will lead to a different dispersion coefficient for a given injection rate.
RESULTS AND DISCUSSION Two-Phase Dispersion Coefficient. By processing the images when the steady-state saturation was attained, water saturation in the flow network for each transport experiment at each REV was calculated. Results are shown in Table S1 of the Supporting Information. To estimate the dispersion coefficient, we estimated the resident concentration spatially averaged over immobile and mobile zones under two-phase conditions, for every snapshot. To calculate the dispersion coefficient, eq 2 is fitted to the experimental data (C versus time) assuming two fitting parameters (D and u). C0 is the maximum average resident concentration calculated from image processing and is variable depending upon the hydrodynamic conditions and saturation profile. The maximum value of the total resident concentration may not necessarily reach the injection concentration. In a few cases, mostly for the ink rate of 1 mL/h, a change in the saturation occurred as a result of the viscous forces. Thus, as a result of this change of the hydrodynamic conditions, experimental results after this point were not considered in the analysis. Therefore, C0 was chosen as a fitting parameter only for those cases highlighted in parts a−c in Table S1 of the Supporting Information. The fitting results turned out to be reliable because the normalized mean square error (NMSE) in all cases was larger than 0.97. NMSE was defined as NMSE = 1 − ∥(Xi − X̂ i)/(Xi − X̅ i)∥2 where, ∥ ∥2 indicates the 2-norm of a vector and Xi and X̂ i are the experimental and fitting results, respectively. Fitted results are also shown in Table S1 of the Supporting Information. The dispersion coefficient versus saturation and the statistical variation of the dispersion coefficient for different ink rates are shown in Figure 2. In all three different ink rates, the dispersion coefficient is variable with respect to saturation. The results are interesting from two aspects: (a) a non-monotonic trend between dispersion and saturation and (b) the scattering of the data, given that the uncertainty in the analysis turned out to be negligible. 4387
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology
Figure 3. Processed images from the experiments visualized at the third REV. The first row shows transport in a fully water-saturated micromodel. The second and third rows show transport in two fluid-phase micromodels, where ink is injected in the water-filled area from the left boundary. All of the images correspond to average concentrations of approximately 50%. The gradient in the concentration is depicted as the gradient in color between deep blue (zero ink concentration) and red (ink concentration of 1). Fluorinert-filled and solid areas are shown in gray and black, respectively. It is clear that some pockets of water-filled areas (in blue) in the second and third rows appear not to contribute significantly to the flow and transport. These areas form dead-end clusters (immobile zones) with almost close to zero advective transport flux.
saturation, and the bottom row shows the flow network at the low water saturation. Snapshots have been selected such that the average resident concentration is around 50%. The first row for fully water-saturated experiments shows that, with increase of the flow rate, the concentration front becomes parabolic as a result of the network side boundary effect, which has smaller pore connectivity. With an increase of the rate, i.e., higher Péclet number, the advective transport becomes stronger and transport along the network is faster than the lateral direction. The transition from red to blue indicates the thickness of the transition zone of the concentration. When the ink flow rate is increased, this transition zone becomes smaller, as expected. In all of these cases, there is no stagnant zone (blue pockets of water-filled areas by the end of the experiments) because everywhere will be filled by ink by flow. The second row shows snapshots with saturation values in the intermediate-to-high region (0.71, 0.61, and 0.79 from left to right for ink flow rates of 0.05, 0.1, and 1 mL/h, respectively). It is clear that some pockets of water-filled areas (in blue) appear not to contribute significantly to the flow and transport, even if the snapshots were taken when the maximum concentrations were achieved by the end of experiments. These areas form dead-end clusters (immobile zones) with almost close to zero advective transport flux. These areas are not fully accessible to ink by advection, which means that an immobile zone is created. Eventually, ink will be transported to these zones by diffusion but after an excessively long time. These dead-end clusters are responsible for the deviation of the transport behavior of the system from a typical Fickian transport, and consequently, the average resident concentration would not reach the injection boundary value. In the third row, where water saturation is 0.22, 0.39, and 0.39 from left to right for ink rates of 0.05, 0.1, and 1 mL/h, respectively, the immobile areas are small and limited. The
It is known that, for small pore velocities (small Péclet numbers), where advective transport is in the order of diffusive transport, the dispersion coefficient will be in the same order as the diffusion coefficient. The Péclet number is defined as Pe = Lu/D, with L being a characteristic length, in our case an REV size, u being the average pore velocity, and D being the diffusion coefficient. Because estimation of the dispersion coefficient is obtained by fitting the advection−dispersion equation, the macroscopic concentration versus time is needed to be estimated. Averaging the concentration over the REV will provide the macroscopic concentration, and thus, length of an REV has been chosen for estimation of the Péclet number. From Figure 2, one can conclude that for the lowest ink rate (0.05 mL/h), the variation of the dispersion coefficient is smaller than the variation of the dispersion coefficient for the other two ink rates (0.1 and 1 mL/h). This is due to the fact that advective transport for the smallest flow rate is weak and the contrast in transport in the two different zones is not as significant as in high flow rate experiments.49 In the case of the ink rate of 1 mL/h, there is a significant increase of the dispersion coefficient in the intermediate saturation compared to the fully water-saturated condition. Note that, for a given ink injection rate, because fluid saturation values are different, the pore-scale Péclet number would be variable because pore velocity is a function of saturation. To visualize the non-monotonicity in immobile zones as an underlying reason for the non-monotonic dispersion−saturation relation, the concentration profiles for different saturation values and ink rates from a number of processed images are shown in Figure 3. Columns 1, 2, and 3 from left to right correspond to the ink flow rate equal to 0.05, 0.1, and 1 mL/h, respectively. The top row shows the fully water-saturated micromodel. The middle row shows the flow network at the intermediate water 4388
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology
Figure 4. Dispersion coefficient plotted against pore velocity for fully water-saturated (open symbols) and two fluid-phase filled (solid symbols) micromodels for three different injection rates. The data points are attained by fitting eq 2 to the average resident concentration versus time data.
immobile zone for small water saturations is small compared to the intermediate-to-high saturation (third row versus second row). Because by design there is no immobile zone in a fully watersaturated system, the observation of these three rows indicates a non-monotonic relation between immobile zone saturation and total water saturation. At a very small percolating water saturation, the immobile zone is small, while under fully saturated conditions, it is zero. The immobile zone saturation reaches a maximum value at intermediate water saturations. As discussed before, with a constant flux at the inlet boundary, saturation influences the pore velocity and dispersion coefficient. To establish a more concrete relation, the data points showing the dispersion coefficient versus pore velocity (both obtained by fitting) have been plotted for fully watersaturated (open symbols) and two fluid-phase (solid symbols) conditions in Figure 4. Again, for a given pore velocity, there is a significant increase in the dispersion coefficient for a twophase situation compared to a single-phase situation. A linear scaling could be a reasonable fit of the log−log plot of dispersion versus pore velocity shown in Figure 4. However, the data scattering increases for high and low ink flow rates. This is translated to a power law relation in a linear presentation. Interaction between Mobile and Immobile Zones. The interaction between mobile and immobile zones can be quantified by their amounts and the mass transfer rate between these two zones. Parameters β and α capture these two aspects, underlying the non-Fickian transport under two-phase conditions. As discussed earlier, all of the functionalities proposed for β and α in the literature have been proposed on the basis of fitting MIM to the experimental data. For the first time in the literature, we use direct data extraction from the images to calculate β and α. The first parameter, β, is calculated by dividing the number of pixels in the mobile region over the number of pixels in the water-occupied region. For a better illustration, the relation between immobile and mobile saturation is shown in Figure 5. Our results for three different
Figure 5. Immobile saturation versus the mobile saturation for three different ink rates in this experiment and results reported in the literature20,28,32 shown in the inset.
ink rates have been compared to the β values resulting from fitting MIM to the experimental data in the literature.23,31,35 The data from the literature were obtained in column-scale experiments. The range of variation of immobile saturation in the literature is smaller than 0.2, whereas we have calculated immobile zones of up to 0.3 saturation. The higher immobile saturation in our micromodel experiments can be due to the two-dimensionality of the micromodel. From our experiments and the literature data, it is evident that the immobile saturation is more likely to be high when the mobile saturation is around the intermediate values. When β was plotted against Sw, no clear trend could be deduced. This can be attributed to the fact that data points that have similar values for saturation do not necessarily have similar spatial configurations in the flow network. These variations in the spatial configuration of each saturation implicitly cause corresponding variations in the tortuosity of the available 4389
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology porous medium. It is known that tortuosity is related to the effective diffusivity, creating different configurations for mobile and immobile regions. The effect of tortuosity becomes more evident at intermediate saturations, where the statistical distribution of available flow paths is wider, creating higher dispersion in the values for β. Tortuosity potentially has an effect also on the average pore-scale velocities and dispersion coefficients. However, we are not in a position to quantify tortuosity at this stage because it would require intensive pore network modeling of each individual case. The results also indicate a general increase in immobile saturation with the increase of the ink rate. Thus, it implies that the immobile saturation is a dynamic property of the system. The results also show that, at intermediate-to-high saturations (Sw ∼ 0.6−0.8), the maximum immobile saturation resulted. This is in agreement with the results of the dispersion coefficient versus saturation in Figure 2, with the maximum dispersion coefficient at a similar range of saturation. That supports the crucial role of immobile water saturation in the deviation of the transport behavior from Fickian to non-Fickian regime under two-phase conditions. Another important aspect of immobile zones is their mass exchange rate with the mobile zones under different dynamic conditions. The mass transfer rate between these two zones has been calculated by the volumetric averaging of concentrations over these regions, denoted by Cm(t) and Cim(t), respectively (refer to the Supporting Information). Therefore, the mass exchange rate between the two regions, α, is calculated as a function of time using eq 5. In addition to the physically based direct estimation of the parameters, we can even evaluate the evolution of the mass transfer rate under transient conditions, which has never been reported before. For a better illustration, results have been plotted against a dimensionless time. To define the dimensionless time, actual time has been divided by the characteristic time. The characteristic time is defined as the time required at a given ink rate to fill that part of the pore space that is occupied by water. Thus, the characteristic time is tchar = VSϕ/qink, and the dimensionless time is defined as τ = qinkt/VSϕ, where V denotes the volume of an REV. Panels a and b of Figure 6 show the mass transfer rate versus time for water injection rates of 0.1 and 1 mL/h, respectively. These panels clearly show that assuming a constant transfer rate is not valid for transient transport conditions. The values of α are higher for the ink injection rate of 1 mL/h by almost an order of magnitude from the two other rates of 0.05 and 0.1 mL/h. Similar proportionality between the mass transfer rate and pore velocity has been reported in a few studies.36,48 With an increasing ink injection rate, concentration fingering and inhomogeneous concentration front take place. That leads to an increase of the interfacial boundary between the mobile and immobile zones. An increase of the interface will eventually lead to the increase of the mass transfer rate. All of the experimental measurements were also fitted using STANMOD50 (results shown in Figures S4 and S5 of the Supporting Information). Even though there is a good correlation between experimental data and fitted results (example shown in Figure S4 of the Supporting Information), the underlying assumptions for the parametrization of MIM are incorrect. The fitted values for β show an unrealistic range, as low as 0.09 (Figure S5 of the Supporting Information), which is not supported by the experimental data (Figure 5). Furthermore, the MIM employs a constant α, which is not
Figure 6. Time evolution of the mass transfer rate between mobile and immobile zones (α) for experiments with water flow rates of (a) 0.1 mL/h and (b) 1 mL/h. In each panel, black, blue, and green represent the ink rates of 0.05, 0.1, and 1 mL/h for REVs 1, 2, and 3. The inset shows the evolution of total, mobile, and immobile concentrations.
the case for transient conditions (Figure 6). The drawbacks of the parametrization assumptions of MIM can be summarized as follows: (a) The immobile saturation fraction is not constant and is changing non-monotonically with the change of saturation. (b) The mass transfer rate is not constant and is definitely changing under transient transport conditions. (c) With the increase of the injection rate (Péclet number), the mass transfer rate increases. (d) With the increase of the injection rate (Péclet number), a larger immobile saturation would be potentially induced.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.5b05947. Visualization setup, fluids and flow rate control, experimental procedure, example of a grayscale image taken from the whole micromodel (using four cameras) 4390
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology
■
(Figure S1), transport of ink in water in a two fluid-phase filled micromodel (Figure S2), relation between the grayscale value and concentration value (Figure S3), comparison between the breakthrough curves resulting from image processing and curve fitting using STANMOD (Figure S4), comparison between the β values estimated from image processing and β values calculated using STANDMOD (Figure S5), estimation of the dispersion coefficient and pore velocity using fitting eq 2 (Table S1), and statistical analysis of the results (Table S2) (PDF)
(8) Levy, M.; Berkowitz, B. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. J. Contam. Hydrol. 2003, 64, 203−226. (9) Bijeljic, B.; Blunt, M. J. Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. Water Resour. Res. 2006, 42, W01202. (10) Seymour, J. D.; Callaghan, P. T. Generalized approach to NMR analysis of flow and dispersion in porous media. AIChE J. 1997, 43, 2096−2111. (11) Karimi, A.; et al. Wettability alteration in carbonates using zirconium oxide nanofluids: EOR implications. Energy Fuels 2012, 26 (2), 1028−1036. (12) Mahani, H.; Berg, S.; Ilic, D.; Bartels, W.-B.; Joekar-Niasar, V. Kinetics of low-salinity-flooding effect. SPE J. 2015, 20 (01), 8−20. (13) Aoudia, M.; Al-Maamari, R. S.; Nabipour, M.; Al-Bemani, A. S.; Ayatollahi, S. Laboratory study of alkyl ether sulfonates for improved oil recovery in high-salinity carbonate reservoirs: A case study. Energy Fuels 2010, 24 (6), 3655−3660. (14) Šimůnek, J.; Jarvis, N. J.; van Genuchten, M. T.; Gärdenäs, A. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 2003, 272 (1), 14−35. (15) Butters, G. L.; Jury, W. A. Field scale transport of bromide in an unsaturated soil: 2. Dispersion modeling. Water Resour. Res. 1989, 25 (7), 1583−1589. (16) Joekar-Niasar, V.; Ataie-Ashtiani, B. Assessment of nitrate contamination in unsaturated zone of urban areas: The case study of Tehran, Iran. Environ. Geol. 2009, 57 (8), 1785−1798. (17) Yule, D. F.; Gardner, W. R. Longitudinal and transverse dispersion coefficients in unsaturated plainfield sand. Water Resour. Res. 1978, 14 (4), 582−588. (18) Khan, A. U.-H.; Jury, W. A. A laboratory study of the dispersion scale effect in column outflow experiments. J. Contam. Hydrol. 1990, 5 (2), 119−131. (19) Maciejewski, S. Numerical and experimental study of solute transport in unsaturated soils. J. Contam. Hydrol. 1993, 14 (3), 193− 206. (20) Forrer, I.; Kasteel, R.; Flury, M.; Flühler, H. Longitudinal and lateral dispersion in an unsaturated field soil. Water Resour. Res. 1999, 35 (10), 3049−3060. (21) Russo, D. Stochastic modeling of macrodispersion for solute transport in a heterogeneous unsaturated porous formation. Water Resour. Res. 1993, 29 (2), 383−397. (22) Russo, D. Stochastic analysis of flow and transport in unsaturated heterogeneous porous formation: Effects of variability in water saturation. Water Resour. Res. 1998, 34 (4), 569−581. (23) Padilla, I. Y.; Yeh, T-CJ; Conklin, M. H. The effect of water content on solute transport in unsaturated porous media. Water Resour. Res. 1999, 35 (11), 3303−3313. (24) Wildenschild, D.; Jensen, K. H. Laboratory investigations of effective flow behavior in unsaturated heterogeneous sands. Water Resour. Res. 1999, 35 (1), 17−27. (25) Bromly, M.; Hinz, C. Non-Fickian transport in homogeneous unsaturated repacked sand. Water Resour. Res. 2004, 40 (7), W07402. (26) Guillon, V.; Fleury, M.; Bauer, D.; Neel, M. C. Superdispersion in homogeneous unsaturated porous media using NMR propagators. Phys. Rev. E 2013, 87 (4), 043007. (27) Maraqa, M. A.; Wallace, R. B.; Voice, T. C. Effects of degree of water saturation on dispersivity and immobile water in sandy soil columns. J. Contam. Hydrol. 1997, 25 (3−4), 199−218. (28) Nützmann, G.; Maciejewski, S.; Joswig, K. Estimation of water saturation dependence of dispersion in unsaturated porous media: Experiments and modelling analysis. Adv. Water Resour. 2002, 25 (5), 565−576. (29) Sato, T.; Tanahashi, H.; Loáiciga, H. A. Solute dispersion in a variably saturated sand. Water Resour. Res. 2003, 39 (6), WR001649. (30) Haga, D.; Niibori, Y.; Chida, T. Hydrodynamic dispersion and mass transfer in unsaturated flow. Water Resour. Res. 1999, 35 (4), 1065−1077.
AUTHOR INFORMATION
Corresponding Author
*Telephone: +44-16130-64867. E-mail: vahid.niasar@ manchester.ac.uk. Author Contributions
Nikolaos K. Karadimitriou designed the setup and performed the experiments and part of image processing. Vahid JoekarNiasar defined the research objectives and produced the curvefitting code to extract dispersion coefficients. Masoud Babaei developed MATLAB codes for image processing, estimated mobile and immobile model parameters from image analyses, and produced the colored snapshots. Nikolaos K. Karadimitriou estimated the mobile and immobile model parameters by STANMOD. Craig A. Shore fabricated the experimental parts. The manuscript was written through contributions of Vahid Joekar-Niasar, Nikolaos K. Karadimitriou, and Masoud Babaei. All authors have given approval to the final version of the manuscript. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors are thankful to the School of Chemical Engineering and Analytical Science of the University of Manchester for funding the micromodel laboratory and this research project.
■
REFERENCES
(1) de Anna, P.; Jimenez-Martinez, J.; Tabuteau, H.; Turuban, R.; Le Borgne, T.; Derrien, M.; Méheust, Y. Mixing and reaction kinetics in porous media: An experimental pore scale quantification. Environ. Sci. Technol. 2014, 48, 508−516. (2) Yang, X.; Zhang, Y.; Chen, F.; Yang, Y. Interplay of natural organic matter with flow rate and particle size on colloid transport: Experimentation, visualization, and modeling. Environ. Sci. Technol. 2015, 49, 13385−13393. (3) Rodrigues, S. N.; Dickson, S. E. The effect of matrix properties and preferential pathways on the transport of Escherichia coli RS2-GFP in single, saturated, variable-aperture fractures. Environ. Sci. Technol. 2015, 49, 8425−8431. (4) Gvirtzman, H.; Gorelick, S. M. Dispersion and advection in unsaturated porous media enhanced by anion exclusion. Nature 1991, 352 (6338), 793−795. (5) Sang, W.; Morales, V. L.; Zhang, W.; Stoof, C. R.; Gao, B.; Schatz, A. L.; Zhang, Y.; Steenhuis, T. S. Quantification of colloid retention and release by straining and energy minima in variably saturated porous media. Environ. Sci. Technol. 2013, 47, 8256−8264. (6) Jenneman, G. E.; Knapp, R. M.; McInerney, M. J.; Menzie, D. E.; Revus, D. E. Experimental studies of in-situ microbial enhanced oil recovery. SPEJ, Soc. Pet. Eng. J. 1984, 24 (01), 33−37. (7) Parker, J. C. Multiphase flow and transport in porous media. Rev. Geophys. 1989, 27 (3), 311−328. 4391
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
Article
Environmental Science & Technology (31) Toride, N.; Inoue, M.; Leij, F. J. Hydrodynamic dispersion in an unsaturated dune sand. Soil Sci. Soc. Am. J. 2003, 67 (3), 703−712. (32) Roth, K.; Hammel, K. Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow. Water Resour. Res. 1996, 32 (6), 1653−1663. (33) Scheidegger, A. E. General theory of dispersion in porous media. Journal of Geophysical Research 1961, 66 (10), 3273−3278. (34) Kirda, C.; Nielsen, D. R.; Biggar, J. W. Simultaneous transport of chloride and water during infiltration. Soil Sci. Soc. Am. J. 1973, 37 (3), 339−345. (35) De Smedt, F.; Wauters, F.; Sevilla, J. Study of tracer movement through unsaturated sand. J. Hydrol. 1986, 85 (1), 169−181. (36) De Smedt, F d; Wierenga, P. J. Solute transfer through columns of glass beads. Water Resour. Res. 1984, 20 (2), 225−232. (37) Matsubayashi, U.; Devkota, L. P.; Takagi, F. Characteristics of the dispersion coefficient in miscible displacement through a glass beads medium. J. Hydrol. 1997, 192 (1), 51−64. (38) Coats, K. H.; Smith, B. D. Dead-end pore volume and dispersion in porous media. SPEJ, Soc. Pet. Eng. J. 1964, 4 (01), 73−84. (39) Gaudet, J. P.; Jegat, H.; Vachaud, G.; Wierenga, P. J. Solute transfer, with exchange between mobile and stagnant water, through unsaturated sand. Soil Sci. Soc. Am. J. 1977, 41 (4), 665−671. (40) Van Genuchten, M. T.; Cleary, R. W. Movement of solutes in soil: Computer-simulated and laboratory results. Dev. Psychiatry 1979, 5, 349−386. (41) Bond, W. J.; Wierenga, P. Immobile water during solute transport in unsaturated sand columns. Water Resour. Res. 1990, 26 (10), 2475−2481. (42) Van Genuchten, M. T.; Wierenga, P. J. Mass transfer studies in sorbing porous media I. Analytical solutions. Soil Sci. Soc. Am. J. 1976, 40 (4), 473−480. (43) Mayes, M. A.; et al. Transport of multiple tracers in variably saturated humid region structured soils and semi-arid region laminated sediments. J. Hydrol. 2003, 275 (3), 141−161. (44) Zinn, B.; Harvey, C. F. When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian hydraulic conductivity fields. Water Resour. Res. 2003, 39 (3), 1051. (45) Gao, G.; Zhan, H.; Feng, S.; Fu, B.; Ma, Y.; Huang, G. A new mobile−immobile model for reactive solute transport with scaledependent dispersion. Water Resour. Res. 2010, 46 (8), W08533. (46) Tran Ngoc, T. D.; Lewandowska, J.; Vauclin, M.; Bertin, H. Two-scale modeling of solute dispersion in unsaturated doubleporosity media: Homogenization and experimental validation. Int. J. Numer. Anal. Methods Geomech. 2011, 35 (14), 1536−1559. (47) Haggerty, R.; Harvey, C. F.; von Schwerin, C. F.; Meigs, L. C. What controls the apparent timescale of solute mass transfer in aquifers and soils? A comparison of experimental results. Water Resour. Res. 2004, 40 (1), W01510. (48) Bajracharya, K.; Barry, D. A. Nonequilibrium solute transport parameters and their physical significance: Numerical and experimental results. J. Contam. Hydrol. 1997, 24 (3), 185−204. (49) Rolle, M.; Hochstetler, D.; Chiogna, G.; Kitanidis, P. K.; Grathwohl. Experimental Investigation and Pore-Scale Modeling Interpretation of Compound-Specific Transverse Dispersion in Porous Media. Transp. Porous Media 2012, 93, 347−362. (50) Feinstein, D. T.; Guo, W. STANMOD: A Suite of WindowsBased Programs for Evaluating Solute Transport. Groundwater 2004, 42 (4), 482−487.
4392
DOI: 10.1021/acs.est.5b05947 Environ. Sci. Technol. 2016, 50, 4384−4392
