# Suppressing viscous fingering in structured porous media

^{a}School of Chemical Engineering and Analytical Science, The University of Manchester, Manchester M13 9PL, United Kingdom;^{b}Soil and Terrestrial Environmental Physics, Department of Environmental Sciences, ETH Zurich, 8092 Zurich, Switzerland;^{c}Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544;^{d}Department of Chemical and Biological Engineering, Princeton University, Princeton, NJ 08544

See allHide authors and affiliations

Edited by Tom C. Lubensky, University of Pennsylvania, Philadelphia, PA, and approved March 30, 2018 (received for review January 13, 2018)

## Significance

Viscous fingering commonly takes place during injection of one fluid that displaces a resident fluid in a porous medium. Fingering normally is promoted where the injected fluid is less viscous than the resident fluid being displaced. We propose a design of a porous medium in the form of an ordered structure to suppress or trigger (depending on the application) viscous fingering in porous media without modifying fluid properties or wettability. We utilize pore-scale direct numerical simulations, state-of-art experiments and analysis to derive predictive tools to evaluate effects of various parameters on controlling viscous fingering in porous media. Moreover, we propose generalized analytical solutions and a phase diagram for the parameter space affecting viscous fingering patterns.

## Abstract

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.

- suppressed viscous fingering
- structured porous media
- microfluidics
- direct numerical simulation
- analytical model

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*A*; all four sides were made from polydimethylsiloxane (PDMS)]. We first filled the device with silicone oil (viscosity ^{−6} and 1.4 × 10^{−5} are shown in Fig. 1*B*. The capillary number is defined as *σ* = 28.2 mN/m is the interfacial tension between the two fluids. The displacement is unstable at both ^{−3}), stable displacement is achieved at the lower

### Numerical Simulations.

Details about the numerical setup and boundary conditions are presented in *Materials and Methods*. Fig. 2 shows the morphology of the displacement patterns obtained from the numerical simulations at different values of

The morphology of the invading fluid−fluid interface reflects the combined effects of ^{−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 (*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 *A* and *B* 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 *A* and *B* 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 *A* and *B* 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*A*), which is referred to as “microfingering” in this study.

To systematically quantify front behaviors observed in Fig. 2, we computed four metrics aimed to characterize front displacement patterns as functions of the prescribed *i*) front fractal dimension *C*), computed using the box-counting method following Shokri et al. (26), which measures the interface roughness; (*ii*) the fluid−fluid interface length *D*) spanning the length of the interface between invading and defending fluids normalized with respect to *iii*) displacement efficiency *E*); and (*iv*) normalized fingertip velocity V (Fig. 3*F*). To calculate _{,} and then normalized with respect to the injection velocity.

As shown qualitatively in Fig. 2, all metrics corresponding to ^{−5} remained insensitive to λ for *C*–*F*. Closer inspection of the results displayed in Figs. 2 and 3 reveals that the maximum value of

### Theoretical Analysis.

Our experimental and numerical results confirm that a prescribed gradient in pore size (

The generalized capillary number *a* to position *b* (see Fig. S1). The stress balance is expressed as (the subscripts correspond to the positions *a* and *b*)

The generalized capillary number **3** represents the generalized capillary 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 *Supporting Information* for further details about the derivation of the stability criterion represented by Eq. **4**. This is an analytical tool that enables us to predict the critical value of the generalized capillary number **4**, it will strongly influence **3**) and therefore whether

Fig. 4 illustrates that the analytically predicted **4** at different viscosity ratio *M*, contact angle θ, and length scale l 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. **3** and **4** 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 *μ*_{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

## 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 *λ* 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 *Supporting Information*.

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 *σ* = 28.2 mN/m is the interfacial tension between the two fluids. We started the experiment at a low capillary number ^{−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 ^{−3} Pa⋅s, displacing an immiscible fluid (defending fluid) of viscosity ^{−1} Pa⋅s. The resulting viscosity ratio of defending fluid with respect to invading fluid was M = 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 l and ^{−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.

## Acknowledgments

We acknowledge the UK Engineering and Physical Sciences Research Council for providing PhD Studentship EP/M506436/1 (to H.S.R.). We also acknowledge the assistance given by IT Services and the use of the Computational Shared Facility at The University of Manchester.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: nima.shokri{at}manchester.ac.uk.

Author contributions: H.S.R. and N.S. designed research; H.S.R., D.O., Y.L., C.-Y.L., N.B.L., S.S.D., H.A.S., and N.S. performed research; H.S.R., D.O., H.A.S., and N.S. contributed new reagents/analytic tools; H.S.R., Y.L., C.-Y.L., N.B.L., and N.S. analyzed data; and H.S.R., D.O., Y.L., C.-Y.L., N.B.L., S.S.D., H.A.S., and N.S. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: Data presented in the figures are available on Zenodo at https://doi.org/10.5281/zenodo.1215581.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1800729115/-/DCSupplemental.

- Copyright © 2018 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Jha B,
- Cueto-Felgueroso L,
- Juanes R

- ↵
- Shalliker RA,
- Catchpoole HJ,
- Dennis GR,
- Guiochon G

- ↵
- ↵
- ↵
- Sahimi M

- ↵
- ↵
- Hekmatzadeh M,
- Dadvar M,
- Sahimi M

- ↵
- Coumou D,
- Drieser T,
- Geiger S,
- Heinrich C,
- Matthai S

- ↵
- Hill S

- ↵
- ↵
- ↵
- ↵
- Lindner A,
- Bonn D,
- Meunier J

- ↵
- ↵
- Holtzman R

- ↵
- Zhao B,
- MacMinn CW,
- Juanes R

- ↵
- ↵
- ↵
- Jackson S,
- Power H,
- Giddings D,
- Stevens D

- ↵
- ↵
- ↵
- ↵
- Xu B,
- Yortsos Y,
- Salin D

- ↵
- Yortsos Y,
- Xu B,
- Salin D

- ↵
- ↵
- Shokri N,
- Sahimi M,
- Or D

- ↵
- ↵
- ↵
- Armstrong R,
- Georgiadis A,
- Ott H,
- Klemin D,
- Berg S

- ↵
- Armstrong RT, et al.

- ↵
- ↵
- Zaretskiy Y,
- Geiger S,
- Sorbie K

- ↵
- Ferrari A,
- Jimenez-Martinez J,
- Borgne T,
- Méheust Y,
- Lunati I

- ↵
- Rabbani HS,
- Joekar-Niasar V,
- Pak T,
- Shokri N

- ↵
- Rabbani HS,
- Joekar-Niasar V,
- Shokri N

- ↵
- Raeini AQ,
- Bijeljic B,
- Blunt MJ

- ↵
- Deshpande S,
- Anumolu L,
- Trujillo M

- Armstrong RT,
- Berg S

- Berg S, et al.

- Osei-Bonsu K,
- Grassia P,
- Shokri N

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences