Suppressing viscous fingering in structured porous media.

Finger-like protrusions that form along fluid-fluid displacement fronts in porous media are often excited by hydrodynamic instability when low-viscosity fluids displace high-viscosity resident fluids. Such interfacial instabilities are undesirable in many natural and engineered displacement processes. We report a phenomenon whereby gradual and monotonic variation of pore sizes along the front path suppresses viscous fingering during immiscible displacement, that seemingly contradicts conventional expectation of enhanced instability with pore size variability. Experiments and pore-scale numerical simulations were combined with an analytical model for the characteristics of displacement front morphology as a function of the pore size gradient. Our results suggest that the gradual reduction of pore sizes act to restrain viscous fingering for a predictable range of flow conditions (as anticipated by gradient percolation theory). The study provides insights into ways for suppressing unwanted interfacial instabilities in porous media, and provides design principles for new engineered porous media such as exchange columns, fabric, paper, and membranes with respect to their desired immiscible displacement behavior.

The unstable growth of fluid−fluid interfacial perturbations has been the subject of a large literature owing to its many applications: for example, fluid mixing in microfluidics (1), chromatographic separation of solvents (2), infiltration of water into soils (3), oil recovery from underground reservoirs (4, 5), carbon dioxide sequestration (6, 7), and the formation of plumes in midocean ridges (8), to list a few. For viscously dominated flows, Hill (9) and Saffman and Taylor (10) were the first to quantify the highly ramified morphology of an interface resulting from displacement of a viscous fluid by a fluid of lower viscosity and so document the emergence of finger-like invasion patterns [viscous fingering (VF)]; Hill (9) investigated the process using a packed bed, while Saffman and Taylor (10) employed fluid-filled Hele-Shaw cells to study VF. An excellent review on VF is provided by Homsy (11). Although the fundamental principles governing interfacial instability are relatively well understood, their manifestation in porous media with rich morphologies of displacement fronts remains an active field of research.

Fluid VF during immiscible displacement in porous media is relevant to a variety of applications. In oil recovery from geologic reservoirs, VF can result in early breakthrough of the invading fluid (often water or brine), thus diminishing the efficiency of oil recovery and at times rendering it uneconomical (4, 5). In environmental applications, VF has been implicated in the potential for early arrival of pollutants to underlying groundwater resources. The technological challenges presented by VF have prompted numerous theoretical and experimental studies (12–20). Some of the studies have shown that the use of non-Newtonian fluids (13) or nonlinear control of injection rate (17) stabilize the fluid−fluid interface. Other studies (14–16) suggest that alteration of wetting properties of the porous medium offers a potential remedy for eliminating VF. However, for many applications, the alteration of wetting characteristics of the porous medium is not trivial; hence other solutions must be developed to control VF.

In this research, we demonstrate the influence of regular pore size variations in a porous medium as a means for suppressing the growth of viscous fingers during immiscible displacement. Such a statement may appear counterintuitive at first glance, because, in the literature, pore size variations are considered to be a factor that enhances the frequency of fluid front tip splitting and thus intensifies the fingering phenomenon (15). We report a structure in the form of an “ordered porous medium” in which pore size varies monotonically along the direction of flow. From the physical point of view, such an ordered porous medium allows simultaneous control over viscous and capillary forces in the same direction, which is otherwise rare in random porous media and has not been studied before. The prescribed structure of the porous media is inspired by the theory of percolation under a gradient introduced by Wilkinson (21) which has been used to describe displacement patterns in porous media (22–24). The work of Xu et al. (23) combines gradient percolation with conventional invasion percolation to derive the now classic phase diagram of Lenormand et al. (25) for the fluid front stability in random porous media. In addition, Yortsos et al. (24) extended the approach to model flow profiles in spatially correlated pore networks. In the studies above, a percolation gradient was introduced in the form of an externally applied pressure drop using a set of local pore filling rules. In particular, Yortsos et al. (24) conjectured that a gradient of pore network correlations could act to trigger or suppress VF, yet, to date, no experimental or theoretical evidence supports this conjecture. Although recent studies (18, 19) have indicated that gradual variation in the thickness of a Hele-Shaw cell can significantly restrain VF, no studies have shown how gradual pore size variation would affect VF in porous media (which is a more complicated system than a Hele-Shaw cell).

In this study, we combine experiments, numerical simulations, and theoretical analysis to demonstrate how a porous medium with ordered pore sizes controls (i.e., triggers or suppresses) VF during immiscible flow in porous media. We begin with a porous medium (uniform and ordered) saturated with high-viscosity defending fluid where low-viscosity invading fluid displaces the defending fluid at various flow rates. Our results show that fluid fronts traversing a porous medium where the pore size is gradually reduced along the flow direction results in a velocity-dependent morphological transformation of the front from unstable to stable. Moreover, we show that stabilization of the invasion front at high injection rates requires an increase in the pore size gradient along the flow path. These results provide a means for inhibiting or triggering VF and interfacial instability in engineered porous materials. The insights gained from this study pave the way to new designs of chromatographic columns, membranes, sensors, and other porous media such that the displacement front morphology is unconditionally stable (under prescribed operations conditions), and improve fundamental understanding of VF in porous media. This study is limited to drainage conditions only.

Results and Discussion

Experiments.

We conducted fluid displacement experiments in micromodels of ordered porous media [Fig. 1; all four sides were made from polydimethylsiloxane (PDMS)]. We first filled the device with silicone oil (viscosity = 200 mPa⋅s), and then injected water (colored with red dye, viscosity = 1 mPa⋅s) at a constant flow rate to displace the silicone oil, which wetted the surface. The experimental results for two different capillary numbers = 7.5 × 10−6 and 1.4 × 10−5 are shown in Fig. 1. The capillary number is defined as , where is the area-averaged (Darcy) velocity, and σ = 28.2 mN/m is the interfacial tension between the two fluids. The displacement is unstable at both for a uniform porous medium (), as expected according to the traditional VF criterion. Note that represents the pore size gradient where is the pore width at the outlet, = 880 μm is the pore width at the inlet, and 29.5 mm is the length of the ordered region. However, for a nonuniform porous medium ( −6.1 × 10−3), stable displacement is achieved at the lower , and becomes unstable at higher . These results suggest that there is a capillary number range where the reduction in pore size along the flow direction can stabilize the displacement of a fluid with higher viscosity during injection of a fluid with lower viscosity. We next turn to numerical simulations to develop more understanding of the macroscopic and microscopic features of these observations.

Fig. 1.

(A) Design of the microfluidic device. The lengths of the ordered medium l and the pore radius r are indicated. (B) Experimental results for Ca = 7.5 × 10−6 and 1.4 × 10−5 for a uniform porous medium with λ = 0 and a nonuniform medium with λ = −6.1 × 10−3. The invading fluid is water with red dye, and the displaced fluid is transparent silicone oil; the flow is from the bottom to the top.

Numerical Simulations.

Details about the numerical setup and boundary conditions are presented in . Fig. 2 shows the morphology of the displacement patterns obtained from the numerical simulations at different values of and . Note that, in the present study, capillary fingering is absent and therefore does not play a role in controlling the displacement patterns.

Fig. 2.

Simulation results demonstrating displacement fluid front morphologies for different values of the capillary number Ca and the pore size gradients λ at the time where the invading fluid reaches the outlet (the direction of displacement is from bottom to top). The white, orange, and black colors represent invading fluid, defending fluid, and the interface, respectively. The viscosities of the invading fluid µ1 and defending fluid µ2 were kept constant at 10−3 Pa⋅s and 10−1 Pa⋅s, respectively. The viscosity ratio of defending fluid over invading fluid M = 100. The results illustrate that fingering is suppressed as λ becomes more negative. Moreover, for λ = −1.0 × 10−3, −3.0 × 10−3, and −4.7 × 10−3, the critical capillary numbers Ca at which the crossover from VF to CD takes place are 3.2 × 10−7, 3.2 × 10−6, and 9.6 × 10−6, respectively.

The morphology of the invading fluid−fluid interface reflects the combined effects of and , as shown in Fig. 2. Results show that, at high displacement rates with = 3.2 × 10−5, there is no stabilizing effect of negative on VF (for the range of values considered), resulting in almost similar displacement patterns for all cases. However, at lower capillary numbers, there exists a critical capillary number () that indicates a transition in front propagation from VF to compact displacement (CD) for a prescribed value of the pore size gradient . The simulation results in Fig. 2 suggest that is related to approximately as (we will expand upon this in Theoretical Analysis).

The experiments and simulations shown in Figs. 1 and 2, respectively, highlight the interplay of capillary and viscous forces on multiphase flow (immiscible) displacement in porous media and the role of the pore size gradient that may affect both forces simultaneously and the resulting front morphology. To provide additional insights, we performed simulations where all of the properties are kept constant except for changing the sign of (in essence, reversing the direction of fluid injection). Typical examples in Fig. 3 depict displacement front patterns in porous media with = −5.6 × 10−3 and = 5.6 × 10−3, respectively. Based on its definition, positive and negative values of correspond to cases where, respectively, either smaller or larger pores are present at the injection location. Inspection of the patterns in Fig. 3 illustrates the dramatic effect of displacement front flow direction with respect to the pore size gradient. Although the pore size distribution, porosity, wetting, and fluid properties were identical in the two cases, Fig. 3 shows that the gradient in pore size relative to the front flow direction resulted in significantly different displacement patterns. For the scenario where = 5.6 × 10−3 with continuously increasing pore sizes along the direction of flow, VF is accentuated as the lower viscosity fluid preferentially flows through the least resistant pathway. In contrast, when the flow direction is reversed and = −5.6 × 10−3, the fluid−fluid interface becomes more stable (compact), and the front spans the entire width of the domain. For a range of negative , we observe local short fingers on the order of pore sizes (an example is presented in Fig. 3), which is referred to as “microfingering” in this study.

Fig. 3.

(A and B) Effects of the sign of λ on the morphology of displacement patterns (white) in porous media. Negative (A) and positive (B) values of λ correspond to the cases when the large and small pores, respectively, are placed at the injection point. In both cases, the capillary number was kept constant at Ca = 3.2 × 10−6, and the viscosities of the invading and defending fluids were constant at 10−3 Pa⋅s and 10−1 Pa⋅s, respectively. The direction of displacement is from bottom to top. The white, orange, and black colors represent invading fluid, defending fluid, and interface, respectively. The observed contrast in the displacement patterns indicates the substantial impact of the sign of λ on the interface morphology and dynamics and the resulting macroscopic response. (C) Fractal dimension D, (D) interface length L normalized with respect to r, (E) displacement efficiency E (%) defined as the ratio of the recovered volume of the defending fluid to its initial volume, and (F) normalized fingertip velocity V.

To systematically quantify front behaviors observed in Fig. 2, we computed four metrics aimed to characterize front displacement patterns as functions of the prescribed and : (i) front fractal dimension (Fig. 3), computed using the box-counting method following Shokri et al. (26), which measures the interface roughness; (ii) the fluid−fluid interface length (Fig. 3) spanning the length of the interface between invading and defending fluids normalized with respect to ; (iii) displacement efficiency (%), defined as the ratio of the recovered volume of the defending fluid relative to its initial volume (Fig. 3); and (iv) normalized fingertip velocity (Fig. 3). To calculate and , each image shown in Fig. 2 was segmented in ImageJ (27) and used to calculate these parameters. The fingertip velocity was directly measured as the ratio of the distance traveled by the interface over time, and then normalized with respect to the injection velocity.

As shown qualitatively in Fig. 2, all metrics corresponding to = 3.2 × 10−5 remained insensitive to for , while, for other values of , the pore size gradient exerted a significant impact on each of the metrics presented in Fig. 3 . Closer inspection of the results displayed in Figs. 2 and 3 reveals that the maximum value of as a function of corresponds to the case when a transition from a stable to an unstable displacement pattern is observed. This maximum value is used to classify the observed patterns as stable or unstable, with the corresponding critical capillary number denoted .

Theoretical Analysis.

Our experimental and numerical results confirm that a prescribed gradient in pore size () along the fluid displacement flow direction can significantly affect the onset of VF. In this section, we develop a generalized capillary number that incorporates the contribution of the pore size gradient and quantifies the relative importance of viscous and capillary forces. Furthermore, using linear stability analysis, we derive a criterion for the conditions that separate stable and unstable displacement patterns. Such a stability criterion can aid in the design of structured porous media to achieve desired results (fingering or stability, depending on the application) under given boundary conditions. A schematic of the porous medium used to formulate the theoretical framework is presented in Fig. S1. Note that, while inertia was considered in the direct numerical simulations (DNSs) using computational fluid dynamics (CFD), following Saffman and Taylor (10) and Al-Housseiny et al. (18), in the analytical analysis presented below, effects of inertia on two-phase displacement were neglected (justified by the relatively slow flows in porous materials). It should be noted that the gravitational force is included in the theoretical analysis.

The generalized capillary number is derived by applying a force balance (28–30) as the interface moves from position a to position b (see Fig. S1). The stress balance is expressed as (the subscripts correspond to the positions a and b)where is the viscous stress over a unit interfacial length, is the gravitational stress that drives the interface forward, and is the local resistive capillary stress. Viscous stresses over a characteristic length can be evaluated as , while capillary stresses are described by the Young−Laplace equation as and gravitational stresses . Expressing (with being the porosity), , and implementing the boundary conditions = 0, = , , and enabled us to develop an equation for a generalized capillary number that quantifies the relative importance of each stress.

The generalized capillary number is defined asSubstituting and , respectively, with the capillary number and Bond number (ratio of gravitational to capillary forces) in the derived expression for and using gives rise toEq. represents the generalized capillary number in the form of the conventional capillary number and Bond number .

Furthermore, using linear stability analysis, we derived an analytical solution capable of distinguishing between the stable and unstable displacement patterns that takes into account the gradient of pore size along the flow direction (among other parameters). The approach we adopted to derive the stability criterion is similar to that described by Saffman and Taylor (10). However, we have modified the dynamic boundary conditions to include the effect of on VF; see for further details about the derivation of the stability criterion represented by Eq. . This is an analytical tool that enables us to predict the critical value of the generalized capillary number that distinguishes stable and unstable displacement patterns. In particular, we findwhere represents the ratio of gravitational to viscous forces, , indicating viscosity ratio of defending fluid over invading fluid, and is the number of pores perpendicular to the displacement direction (transverse direction). When , the displacement pattern is stable and is unstable. It is important to note that, although contact angle is not included in Eq. , it will strongly influence (Eq. ) and therefore whether or will be dependent upon . The theoretical results illustrating how varies as a function of (both positive and negative) under different capillary numbers along with the stability criterion are presented in Fig. 4.

Fig. 4.

Semilog plot showing the relationship between the generalized capillary number Ca* and λ. The solid line indicates the analytically predicted critical value of the generalized capillary number, which is a solution of the stability criterion derived from linear stability analysis. The symbols are the results obtained by the direct numerical simulations, where filled, half-filled, and open symbols represent the stable front, critical value of the generalized capillary number, and unstable front, respectively. Yellow and red regions mark analytically predicted stable and unstable regions, respectively, separated by the solid line predicted by Eq. . Insets illustrate the displacement front morphology at Ca = 3.2 × 10−6 computed by the direct numerical simulation.

Fig. 4 illustrates that the analytically predicted agrees well with the numerical results in distinguishing between stable and unstable displacement fronts for all combinations of the capillary number and pore size gradient . The validation of Eq. at different viscosity ratio M, contact angle , and length scale is also presented in Figs. S2–S4 respectively. The slight discrepancy in the classification of some points is attributed to the simplifying assumptions made for the derivation of the analytical stability criterion [e.g., ignoring thin wetting films (18) and trapped fluids behind the displacement front]. Examination of Eqs. and suggests that during imbibition (the displacement of a nonwetting phase by a wetting phase), a positive would delay the onset of VF, whereas unstable fronts would always persist for negative ; such a conclusion is experimentally supported by the results of Al-Housseiny et al. (18).

An important result of our simulations is that, for the same capillary number, when , the competition between capillary and viscous forces remains the same (as reflected in the constant value of for in Fig. 4); therefore, VF continues. However, for and the same capillary number, increases, which results in a stable displacement front. The results presented in Fig. 4 further show that there is an analogy between viscosity of the invading fluid and negative According to the phase diagram introduced by Lenormand et al. (25), increasing μ1 transforms the invasion behavior from VF to a stable regime. Similarly, our results indicate that increasing the gradient of pore size (more negative values) stabilizes the displacement front. Therefore, the overall trend observed in Fig. 4 suggests that it is the viscous dissipation that governs stability of the displacement front (due to increase in as gets more negative), in agreement with the previous investigation (20). Negative (meaning decreasing pore size along the flow) enhances the viscous forces. This further promotes lateral growth in disturbances and results in development of a smooth invasion pattern.

Conclusion

Our results demonstrate the impact of λ on the nature of immiscible displacement in porous media. We show that the VF, which is traditionally considered as a function of flow rate, viscosity ratio, and wetting properties of porous media, is controlled by the pore size gradient λ as well. Depending upon the wettability of the porous medium, for a given of the invasion process, both positive and negative λ can inhibit or trigger the growth of viscous fingers. Our numerical and experimental analyses at the pore scale enabled us to identify two pore-scale invasion mechanisms responsible for suppressing VF. More detailed discussions are presented in .

In this research, we have employed a design of a porous medium in the form of an ordered structure to suppress VF. This study has implications in a number of industrial applications, from the design of stable exchange porous columns for analyses and separation science to designing new membranes and porous products for suppression of spurious VF. We envision potential applications related to optimization of reactant transport and phase distribution in fuel cells, sensors and control of fluid flow in spacecraft under microgravity (31), and more. In addition, this research may also contribute toward reconciling pore-scale flow behavior with capillary dispersion phenomena observed during immiscible displacement at the continuum scale (6).

Materials and Methods

Experimental Setup.

The PDMS microfluidic device was made by photolithography. Positive photoresist and plasma etching were used to make the silicone mold for the PDMS to obtain uniform height of the channels. The ratio between the cross-linker and the elastomeric base was chosen as 1.5:10 to enhance the stiffness of the channels. The finished channel was hydrophobic and oleophilic. The triangular area at the inlet (Fig. 1) was designed for stabilizing the interface before it reached the porous medium. The displaced fluid was phenylmethylsiloxane oligomer (PDM-7050) purchased from Gelest Inc. The invading fluid was deionized water mixed with 0.1 wt % food dye for visualization. Considering the small weight ratio of the dye, its effects on the water viscosity and the water−oil interfacial tension were negligible.

The microfluidic device consists of pillar arrays with height H = 160 μm and variable pillar diameter spanning the width of the ordered region w = 30 mm. The pillar diameters and pores were ordered along the direction of the flow, with a pore size gradient . We denote by a reduction in pore size along the flow path, and vice versa for . The capillary number is defined as , where is the area-averaged (Darcy) velocity and σ = 28.2 mN/m is the interfacial tension between the two fluids. We started the experiment at a low capillary number = 6.1 × 10−7 until a stable interface reached the first row of the pillars. Then, the flow rate was increased to a specified value, and the time evolution of the displacement process was recorded by a Nikon camera.

Numerical Setup.

DNS where volume-of-fluid method (interface tracking approach) is coupled with a Navier−Stokes equation has emerged as a powerful tool for diagnosing pore-scale multiphase flow problems with complex boundary conditions (32–35), enabling parameterization of macroscopic quantities (36). In the present study, we utilized DNS within a CFD framework to investigate how the proposed pore size arrangement influences the general dynamics of two-phase flow in porous media and stability of the displacement front. Additional details regarding the numerical algorithm employed in this study are provided in Deshpande et al. (37) and Rabbani et al. (35).

For the 2D simulations performed in the present study, we assumed an invading fluid of viscosity = 10−3 Pa⋅s, displacing an immiscible fluid (defending fluid) of viscosity = 10−1 Pa⋅s. The resulting viscosity ratio of defending fluid with respect to invading fluid was = 100. The contact angle between interface and the solid surface measured along the defending fluid was kept uniform at 30° (i.e., the defending fluid acts as the wetting phase). The values of and were kept constant at 8 mm and 0.17 , respectively, while was allowed to vary based on the simulated value. In the case of simulations, the width of ordered region is equal to . The simulations were performed at several capillary numbers ranging from 3.2 × 10−7 to 3.2 × 10−5 and the pore size gradients ranging from 6.5 × 10−3 to −6.5 × 10−3, respectively. The data, code, and materials used in this analysis will be available freely via sending a request to the corresponding author.