1 Introduction
Fluid invasion into porous media to displace another immiscible fluid is an important process in many practical applications, such as enhanced oil recovery (Morrow & Mason Reference Morrow and Mason2001), geological carbon sequestration (Benson & Cole Reference Benson and Cole2008), groundwater contamination by non-aqueous liquids (Dawson & Roberts Reference Dawson and Roberts1997) and the design of fuel cells and microfluidic devices (Chapuis et al.
Reference Chapuis, Prat, Quintard, Chane-Kane, Guillot and Mayer2008; Anderson, Zhang & Ding Reference Anderson, Zhang and Ding2010; Lee et al.
Reference Lee, Gupta, Hatton and Doyle2017). For such flow behaviour of multiple fluids with an interface in porous media, the instability of the displacement front directly impacts the oil recovery efficiency, the
$\text{CO}_{2}$
storage capacity and the rate of mass transfer between phases. In enhanced oil recovery, due to the unstable water–oil displacement front, water-flooding in natural reservoirs can produce only 10 %–20 % of the initial oil in place (Van’t Veld & Phillips Reference Van’t Veld and Phillips2010). In geological
$\text{CO}_{2}$
sequestration, when the injection of liquid
$\text{CO}_{2}$
stops and brine flows back to displace supercritical
$\text{CO}_{2}$
, the occurrence of fingering flow can significantly increase the efficiency of storage by increasing interfacial area and improving capillary trapping as well as dissolution trapping (Wang et al.
Reference Wang, Zhang, Wei, Oostrom, Wietsma, Li and Bonneville2012; Bachu Reference Bachu2015). Understanding and controlling the immiscible fluid–fluid displacement in porous media is therefore critical for optimizing fluid management.
The instability of the displacement front is a classic problem and continues to be the focus of an enormous number of experimental, theoretical and numerical studies over the past five decades (Saffman & Taylor Reference Saffman and Taylor1958; Paterson Reference Paterson1981; Måløy, Feder & Jøssang Reference Måløy, Feder and Jøssang1985; Lenormand, Touboul & Zarcone Reference Lenormand, Touboul and Zarcone1988; Dvraam & Payatakes Reference Dvraam and Payatakes1995; Babchin et al.
Reference Babchin, Brailovsky, Gordon and Sivashinsky2008; Cottin, Bodiguel & Colin Reference Cottin, Bodiguel and Colin2011; Armstrong & Berg Reference Armstrong and Berg2013; Bischofberger, Ramachandran & Nagel Reference Bischofberger, Ramachandran and Nagel2014; Hu et al.
Reference Hu, Wan, Kim and Tokunaga2017b
, Reference Hu, Wu, Yang and Chen2018b
; Singh et al.
Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017; Rabbani et al.
Reference Rabbani, Or, Liu, Lai, Lu, Datta, Stone and Shokri2018). When gravity can be neglected, the competition between capillarity and viscous force controls the instability of the displacement front, which results in the patterns ranging from capillary fingering to viscous fingering to compact displacement (Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988; Zhang et al.
Reference Zhang, Oostrom, Wietsma, Grate and Warner2011; Chen et al.
Reference Chen, Fang, Wu and Hu2017, Reference Chen, Wu, Fang and Hu2018). This competing effect becomes more complicated when the wettability and the disorder of porous media are both involved (Alava, Dubé & Rost Reference Alava, Dubé and Rost2004; Singh et al.
Reference Singh, Jung, Brinkmann and Seemann2019). The wettability, denoted by the invading fluid contact angle
$\unicode[STIX]{x1D703}$
, represents the affinity of fluid to the pore surface, which directly modifies the local pore-filling events via changing capillary force governed by the Young–Laplace law (Lenormand, Zarcone & Sarr Reference Lenormand, Zarcone and Sarr1983; Zacharoudiou et al.
Reference Zacharoudiou, Chapman, Boek and Crawshaw2017) and thus impacts the overall displacement patterns (Cieplak & Robbins Reference Cieplak and Robbins1988; Holtzman & Segre Reference Holtzman and Segre2015). On the other hand, the pore-scale disorder,
$\unicode[STIX]{x1D706}$
, which represents the degree of randomness of pore size (Chen & Wilkinson Reference Chen and Wilkinson1985), would also impact the local pore-filling paths via changing the threshold capillary force that needs to be overcome by the pressure drop. Therefore, the challenge in characterizing the competition between disorder and wettability on the displacement patterns in porous media is how to link the pore-filling events that highly depend upon capillary force to the non-local invasion behaviour.
To unravel the fundamentals of fluid–fluid displacement processes in disordered porous media under various wetting conditions, theoretical and experimental studies have been conducted and are still on-going (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990; Hecht & Taitelbaum Reference Hecht and Taitelbaum2004; Cottin et al.
Reference Cottin, Bodiguel and Colin2011; Holtzman & Segre Reference Holtzman and Segre2015; Trojer, Szulczewski & Juanes Reference Trojer, Szulczewski and Juanes2015; Jung et al.
Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Singh et al.
Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017; Hu et al.
Reference Hu, Wan, Yang, Chen and Tokunaga2018a
). The pioneering work that links the pore-filling events to the non-local fluid invasion was proposed by Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990). They proposed three meniscus-motion modes (burst, touch and overlap), able to satisfactorily capture the pore-scale fluid displacement, to elucidate the wettability effect in disordered media, which are now widely recognized (Jung et al.
Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Singh et al.
Reference Singh, Scholl, Brinkmann, Michiel, Scheel, Herminghaus and Seemann2017). These works provide a basic understanding of the wettability effect: namely that, as the pore surface becomes more wetting to the invading fluid within the range of
$45^{\circ }<\unicode[STIX]{x1D703}<180^{\circ }$
, the displacement pattern becomes more stable. The critical contact angle
$\unicode[STIX]{x1D703}_{c}$
that separates the unstable from the stable flow regimes is found to depend upon the porosity of the porous medium and the capillary number (Cieplak & Robbins Reference Cieplak and Robbins1988; Hu et al.
Reference Hu, Wan, Yang, Chen and Tokunaga2018a
). Given the coexisting influence of wettability and disorder on local fluid displacement, this wettability effect should vary with disorder; for instance, the critical contact angle
$\unicode[STIX]{x1D703}_{c}$
varies with disorder
$\unicode[STIX]{x1D706}$
.
Although intensive research has focused on the effect of wettability on fluid displacement, less progress has been made in examining how the competition between pore-scale disorder and wettability controls the multiphase flow. In cases where the wettability effect is isolated, previous studies have shown that decreasing disorder stabilizes the displacement for both drainage and imbibition (Chen & Wilkinson Reference Chen and Wilkinson1985; King Reference King1987; Toussaint et al.
Reference Toussaint, Løvoll, Méheust, Måløy and Schmittbuhl2005) and that the disorder also modifies the critical capillary number
$Ca_{c}$
corresponding to the cross-over from capillary to viscous fingering (Yortsos, Xu & Salin Reference Yortsos, Xu and Salin1997; Holtzman & Juanes Reference Holtzman and Juanes2010; Xu et al.
Reference Xu, Ok, Xiao, Neeves and Yin2014; Liu, Zhang & Valocchi Reference Liu, Zhang and Valocchi2015b
). Recently, a systematic study on the impact of disorder and its coupling with wettability was conducted by Holtzman (Reference Holtzman2016). Based on scaling analysis, Holtzman (Reference Holtzman2016) introduced a capillary number
$Ca^{\ast }$
dependent on disorder and contact angle to consider the competition between disorder and wettability effect, showing that the critical capillary number
$Ca^{\ast }$
decreases with increasing
$\unicode[STIX]{x1D706}$
for the range of
$Ca>10^{-5}$
. Holtzman (Reference Holtzman2016) also showed that in the slowest flow rate considered
$(Ca\approx 10^{-5})$
, higher disorder of the flow geometry would enhance overlap events that smooth the local fluid–fluid interface, which is counter to the widely recognized effect that higher disorder destabilizes displacement patterns (Koiller, Ji & Robbins Reference Koiller, Ji and Robbins1992). Given that, in the limit of quasi-static conditions, multiphase flow is significantly affected by the flow geometry, these counteracting effects that stabilize or destabilize would be amplified. Therefore, how the disorder impacts the quasi-static displacement pattern depends on which one of these counteracting effects is dominant. These effects can be clearly represented in a phase diagram in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane, which has not yet been reported.
Here, we aim to propose a phase diagram of quasi-static fluid displacement to elucidate how the competing effects between disorder and wettability control the fluid invasion pattern in porous media. It remains a great challenge to propose such a phase diagram because the number of arrangements of posts (or flow geometry) within disordered porous media can be infinite for any given disorder
$\unicode[STIX]{x1D706}$
, because the radius of posts
$r_{i}$
can take any values following an assumed probability distribution, such as
$r_{i}\sim U[(1-\unicode[STIX]{x1D706})\bar{r},(1+\unicode[STIX]{x1D706})\bar{r}]$
, where
$U$
denotes the uniform distribution,
$\bar{r}$
is a constant and
$\unicode[STIX]{x1D706}$
is the disorder of the porous medium (Chen & Wilkinson Reference Chen and Wilkinson1985; Holtzman & Segre Reference Holtzman and Segre2015). For the classic and widely used phase diagram (Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988), the boundaries that separate the different flow regimes can be uniquely determined from the given flow rate and fluid properties within a specific flow geometry. However, because of the infinite flow geometries, the lower and upper bounds, rather than the specific boundaries, are established to separate the different flow regimes in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane. The zone bounded by the lower and upper limits provides the first identification of displacement patterns. This is very important because the statistical parameters, such as porosity, average particle/pore sizes and disorder, etc., are always easy to obtain for porous media. The flow regimes can then be easily determined via our phase diagram using these statistical parameters rather than the details of pore structure, which are difficult to obtain.
We combine microfluidic experiments, numerical simulations and theoretical analysis to describe the displacement patterns as a function of disorder and contact angle. We fabricate microfluidic chips with four different disorders and image the fluid displacement under extremely slow flow conditions with a microscope and a complementary metal oxide semiconductor (CMOS) camera (§ 2.1). We develop a pore-scale numerical procedure originally proposed by Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) to consider the rectangular geometry and constant-flow-rate conditions (§ 2.2). Based on experimental observations and numerical simulations (§ 3.1), we determine the flow regimes in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane for the flow geometries of the microfluidic chips (§ 3.2). Further, from the viewpoint of probability (§ 3.3.1), we derive a theoretical model that corresponds to the lower and upper bounds of flow regimes, and then propose a phase diagram in terms of disorder and wettability (§ 3.3.2). Finally, the proposed phase diagram is evaluated by extensive numerical simulations, our microfluidic experiments and other existing experimental results (§ 3.3.3).
Figure 1. Microfluidic experimental set-up. (a) Schematic diagram of the microfluidic–microscopy system. (b,c) The flow geometry of the microfluidics is constructed by placing non-overlapping and variable-sized posts (
$r_{i}$
) on a triangular lattice with a spacing
$a=250~\unicode[STIX]{x03BC}\text{m}$
, composed of 987 posts. The radii of the posts follow a uniform distribution. The length (
$L$
), width (
$W$
) and depth (
$D_{0}$
) of the microfluidics are, respectively, 10 mm, 5 mm and
$40~\unicode[STIX]{x03BC}\text{m}$
. (d–g) Four microfluidics set-ups are fabricated with disorder
$\unicode[STIX]{x1D706}=0$
, 0.1, 0.2 and 0.3. The throat size distributions of the flow geometries with four disorders are, respectively, presented on the right side of panels (d–g). The porosity
$\unicode[STIX]{x1D719}$
and the average throat size
$\bar{d}$
for all of the microfluidics are nearly the same, i.e.
$\unicode[STIX]{x1D719}=0.55$
and
$\bar{d}=75~\unicode[STIX]{x03BC}\text{m}$
.
2 Materials and methods
2.1 Experimental section
2.1.1 Microfluidic visualization system
We designed a microfluidic visualization system to perform fluid displacement experiments (figure 1
a). The experimental set-up consists of a microfluidic flow cell, an imaging system and a syringe pump (Harvard Apparatus 70C3007). The syringe pump is used to control the flow rate. The imaging system, connecting to a computer, consists of an inverted microscope (Carl Zeiss, Observer Z1.m) that records images at the pore scale via a charge-coupled device (CCD) camera (Carl Zeiss, AxioCam MRc5), and a CMOS camera (Manta G-1236C, AVT) to record images of the entire domain of the microfluidics with a spatial resolution of
$39.0~\unicode[STIX]{x03BC}\text{m}~\text{pixel}^{-1}$
. The microscope includes a reflecting light objective (Epiplan
$5\times /0.13\,\text{W}0.8^{\prime \prime }$
, working distance = 20.5 mm) to visualize the fluid distribution within the range of
$1.5~\text{mm}\times 1.2~\text{mm}$
with a spatial resolution of
$2.0~\unicode[STIX]{x03BC}\text{m}~\text{pixel}^{-1}$
. The entire experimental set-up is maintained at room temperature (
$20\pm 0.5\,^{\circ }\text{C}$
).
2.1.2 Microfluidic fabrication
We fabricate the microfluidics with polydimethylsiloxane (PDMS) (Sylgard 184, Dow-Corning, USA). First, we construct the flow geometry of the microfluidics by placing non-overlapping and variable-sized posts on a triangular lattice with a spacing
$a=250~\unicode[STIX]{x03BC}\text{m}$
(figure 1
b,c). The post radius
$r_{i}$
follows a uniform distribution
$r_{i}\sim U[(1-\unicode[STIX]{x1D706})\bar{r},(1+\unicode[STIX]{x1D706})\bar{r}]$
, where
$\bar{r}$
is the average radius,
$\bar{r}=87.5~\unicode[STIX]{x03BC}\text{m}$
. The manufacturing resolution is
$1~\unicode[STIX]{x03BC}\text{m}$
, and thus the roughness of post surfaces is approximately
$1~\unicode[STIX]{x03BC}\text{m}$
. Four flow geometries are generated with
$\unicode[STIX]{x1D706}=0$
, 0.1, 0.2 and 0.3 (figure 1
d–g). Then, we generate the corresponding silicon masters using conventional photolithography techniques, including (1) spin-coating a negative photoresist (SU8-2035, MicroChem, USA) to achieve a film thickness of
$40~\unicode[STIX]{x03BC}\text{m}$
onto a
$4^{\prime \prime }$
silicon wafer at 2000 revolutions per minute for 30 s, and (2) exposing the photoresist to ultraviolet (UV) light with a photomask to generate the patterned silicon master. Finally, we obtain a PDMS patterned plate by pouring PDMS onto the silicon master, and then create a microfluidic flow cell by bonding this PDMS patterned plate with a smooth PDMS plate using a Plasma Cleaner (PDC-002, Harrick Plasma, USA). Four different patterns of microfluidics (figure 1
d–g), composed of 987 posts, are generated. The pore volumes (
$PV$
) of the four microfluidics are almost the same,
$PV=1.1~\unicode[STIX]{x03BC}\text{l}$
, and the porosity is
$\unicode[STIX]{x1D719}=0.55$
. We calculate the spatial correlation length of the flow geometry. Since the posts are located on triangular lattices, the correlation lengths for the microfluidics are the same, i.e.
$l_{c}=1.3a$
.
2.1.3 Experimental procedure
Given that the pore surface of microfluidics is hydrophobic, for the purpose of imbibition experiments using the water–air fluid pair, the invading fluid of air is the wetting phase and the degassed water is the defending (non-wetting) phase. We measure the contact angle between water and air on the smooth PDMS plate with a Drop Shape Analyzer (DSA25; Krüss), and the average invading fluid (air) contact angle is
$\unicode[STIX]{x1D703}=67^{\circ }$
. The water is dyed with a light-absorbing dye (Carmine, Wilton) at 0.2 wt.% concentration to increase the signal intensity between the two phases. To conduct a weak imbibition experiment, we first fully saturate the microfluidic with the withdraw mode of the syringe pump at a constant flow rate of
$1~\unicode[STIX]{x03BC}\text{l}~\text{min}^{-1}$
via the flow path
$\text{D}\rightarrow \text{B}\rightarrow \text{A}$
for 5 min (figure 1
a). Afterwards, the three-way valve B is switched to C, and the system is equilibrated at
$20\,^{\circ }\text{C}$
for 15 min. After that, imbibition processes are initiated. Air is injected into the microfluidic at the extremely slow flow rate of
$0.025~\unicode[STIX]{x03BC}\text{l}~\text{min}^{-1}$
with the withdraw mode of the pump via the flow path
$\text{C}\rightarrow \text{B}\rightarrow \text{A}$
. Images of invasion morphology for the entire domain of the microfluidic are recorded at
$10~\text{frames}~\text{min}^{-1}$
with the CMOS camera until the invading fluid reaches the outlet. Higher-resolution images of fluid distributions at the pore scale are also recorded with the microscope. The above procedures are applied to the microfluidics with
$\unicode[STIX]{x1D706}=0$
, 0.1, 0.2 and 0.3, and each disorder condition is repeated four times, with a total of 16 imbibition experiments. Each microfluidic is used only for one time. The procedure for image post-processing has been reported in the previous study (Hu et al.
Reference Hu, Wan, Kim and Tokunaga2017b
). To confirm that the flow-rate condition corresponds to the quasi-static state, we calculate the capillary number
$Ca$
, which is widely used to quantify the relative effect of viscous to capillary forces (Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988; Zhang et al.
Reference Zhang, Oostrom, Wietsma, Grate and Warner2011; Chaudhary et al.
Reference Chaudhary, Cardenas, Wolfe, Maisano, Ketcham and Bennett2013; Geistlinger et al.
Reference Geistlinger, Ataei-Dadavi, Mohammadian and Vogel2015). Here,
$Ca$
is calculated by
$Ca=v_{i}\unicode[STIX]{x1D707}_{i}/\unicode[STIX]{x1D70E}$
, where
$v_{i}=Q/A_{c}$
,
$Q=0.025~\unicode[STIX]{x03BC}\text{l}~\text{min}^{-1}$
,
$A_{c}$
is the cross-sectional area of the inlet,
$A_{c}=W\times D_{0}=5~\text{mm}\times 40~\unicode[STIX]{x03BC}\text{m}$
(figure 1
b),
$\unicode[STIX]{x1D707}_{i}=1.83\times 10^{-5}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$
and
$\unicode[STIX]{x1D70E}=7.25\times 10^{-2}~\text{N}~\text{m}^{-1}$
. Thus, the capillary number under such flow-rate conditions is
$Ca=4.2\times 10^{-10}$
, which generally corresponds to the quasi-static condition.
2.2 Numerical simulation
We develop the model of Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) to simulate the quasi-static fluid–fluid displacement processes, starting with the basic modes for menisci motion in § 2.2.1 and the numerical implementation in § 2.2.2.
2.2.1 Modes for menisci motion
Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) introduced three basic modes for menisci motion, i.e. burst, touch and overlap. These modes, able to capture the local fluid displacement, are fundamental to the pore-scale numerical simulation for fluid inversion into porous media, and also the basis for the derivation of the phase diagram in this work. From geometry analysis, we can obtain the onsets of burst, touch and overlap for non-uniform post arrangements.
The burst event occurs when no stable arc can exist in the current post arrangement with contact angle and local pressure drop. As the local pressure drop increases, the radius of curvature decreases. Once the radius of curvature reaches its minimum, an infinitesimal increment of pressure drop would cause the invading fluid to burst into the adjacent pores, also known as Haines jump (Berg et al.
Reference Berg, Ott, Klapp, Schwing, Neiteler, Brussee, Makurat, Leu, Enzmann and Schwarz2013). Thus, the critical radius of curvature,
$R_{b}$
, for the occurrence of the burst mode corresponds to its minimum. As shown in figure 2(a), based on geometry analysis (see § A.1), we have
$$\begin{eqnarray}\text{burst:}\quad R_{b}=\frac{-b_{1}+\sqrt{b_{1}^{2}-4a_{1}c_{1}}}{2a_{1}},\end{eqnarray}$$
where
$a_{1}$
,
$b_{1}$
and
$c_{1}$
are the coefficients given in table 1 and (A 3).
Figure 2. The geometry to calculate to radii of curvature of the menisci for the three basic modes as introduced by Cieplak & Robbins (Reference Cieplak and Robbins1988). (a) Burst. The onset of a burst event corresponds to the state that for the minimum radius of curvature (the maximum capillary force) no stable meniscus can exist in such an arrangement of the posts. (b) Touch. The occurrence of touch means that the meniscus has touched the edge of the nearest post before the radius of curvature reaches its minimum. (c) The order of occurrence for burst and touch before determining overlap. (d) Overlap. Two menisci merge into one at the location of the three-phase contact line. The arrow indicates the direction of menisci motion, and the red arc indicates the fluid–fluid interface.
For the touch event, the radius of curvature continues to decrease and then touches the edge of the nearest post before reaching its minimum (figure 2 b). Again, based on geometry analysis on figure 2(b) (see § A.2), the condition for the occurrence of touch is given as
$$\begin{eqnarray}\text{touch:}\quad R_{t}=\frac{-b_{2}+\sqrt{b_{2}^{2}-4a_{2}c_{2}}}{2a_{2}},\end{eqnarray}$$
where
$a_{2}$
,
$b_{2}$
and
$c_{2}$
are the coefficients given in table 1 and (A 6).
Note that during numerical simulation, one needs to identify which mode occurs first. Hecht & Taitelbaum (Reference Hecht and Taitelbaum2004) showed that the effect of the order of occurrence for the three basic modes on the multiphase flow is not very important. Holtzman & Segre (Reference Holtzman and Segre2015) check for the burst first and then check for the touch and overlap events, which is employed in this study. As shown in figure 2(c), we consider a specific fluid reconfiguration in which the point
$A$
is located on
$O_{1}O_{2}$
and the arc
$BC$
touches the nearest post, indicating that there is a critical radius of curvature corresponding to the two modes. To determine the order of occurrence of the two modes, we determine the critical radius
$R_{c}$
,
$$\begin{eqnarray}R_{c}=\sqrt{{\textstyle \frac{1}{2}}a^{2}+{\textstyle \frac{1}{2}}(r_{1}^{2}+r_{2}^{2})+R_{b}^{2}-(r_{1}+r_{2})R_{b}\cos \unicode[STIX]{x1D703}}-r_{3},\end{eqnarray}$$
where
$R_{b}$
is the critical radius given by (2.1). If
$R_{b}<R_{c}$
, then the burst has occurred before the arc has touched the nearest post (touch). Conversely, if
$R_{b}>R_{c}$
, then the touch has occurred before the burst.
The overlap mode involves two neighbouring menisci merging into a new one at the three-phase contact point or at the fluid–fluid interface (Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018; Singh et al.
Reference Singh, Jung, Brinkmann and Seemann2019). Primkulov et al. (Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018) have shown that for larger post–post spacing with
$\unicode[STIX]{x1D703}\geqslant 110^{\circ }$
, the percentage of overlap events will be underestimated when considering the two menisci merging at the three-phase contact point. We check the two treatments, and the results show that the difference in the probability of overlap (see (3.3)) is no more than 0.05. Here, as shown in figure 2(d), we consider the former case that is widely used (Koiller et al.
Reference Koiller, Ji and Robbins1992; Holtzman & Segre Reference Holtzman and Segre2015). For this post arrangement, the overlap will occur if
$\unicode[STIX]{x1D702}_{2}+\unicode[STIX]{x1D702}_{3}\geqslant \angle O_{1}O_{2}O_{3}$
. Through geometry analysis (see § A.3), we have the condition
$$\begin{eqnarray}\text{overlap}:\quad \unicode[STIX]{x1D702}_{2}(r_{1},r_{2},R,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D702}_{3}(r_{2},r_{3},R,\unicode[STIX]{x1D703})\geqslant \unicode[STIX]{x1D719}_{0},\end{eqnarray}$$
where the expressions for
$\unicode[STIX]{x1D702}_{2}$
and
$\unicode[STIX]{x1D702}_{3}$
are given in table 1 and (A 11)–(A 12), and
$R$
is the radius of curvature for the burst or touch, depending on which mode occurs first. If the burst occurs first, i.e. if
$R_{b}<R_{c}$
, then
$R=R_{b}$
; while if
$R_{b}>R_{c}$
, then
$R=R_{t}$
, and
$\unicode[STIX]{x1D719}_{0}$
is the angle of
$\angle O_{1}O_{2}O_{3}$
. For the triangular lattice considered in this work (figure 1
c),
$\unicode[STIX]{x1D719}_{0}$
is set as
$120^{\circ }$
.
2.2.2 Numerical implementation
We develop the method of Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) to consider a rectangular domain with constant-flow-rate condition, different from the previous modes, which are applicable for a circular domain and constant-pressure conditions (Cieplak & Robbins Reference Cieplak and Robbins1990; Koiller et al.
Reference Koiller, Ji and Robbins1992). To reproduce the observed quasi-static fluid displacement, we create an irregular rectangular two-dimensional (2-D) model for the microfluidics (figure 1
d–g). The inlet (the left side of figure 1
b) is fixed with a constant flow rate of
$Q_{0}=0.025~\unicode[STIX]{x03BC}\text{l}~\text{min}^{-1}$
, and zero pressure and no-flow conditions are, respectively, applied on the outlet (the right side) and the lateral boundaries (the top and bottom sides). First, all of the pores are saturated with the defending fluid, and the two rows of pores near the inlet are filled with the invading fluid. Then, the Stokes-flow-based continuity equation is assembled for the system, i.e.
$\boldsymbol{Q}=\unicode[STIX]{x1D646}\boldsymbol{P}$
, where
$\boldsymbol{Q}$
is the vector of flux into the pores,
$\boldsymbol{P}$
is the vector of pressure for each pore, and
$\unicode[STIX]{x1D646}$
is the matrix that represents the conductance for adjacent pores of the system. After solving the equation, i.e.
$\boldsymbol{P}=\unicode[STIX]{x1D646}^{-1}\boldsymbol{Q}$
, the pressure for each pore is updated, and thus the radius of curvature for all arcs. After that, the local pore-filling events, i.e. burst, touch and overlap, are identified with (2.1), (2.2) and (2.4), to update the configuration of menisci. Once the menisci advancements have been tracked, the conductance for adjacent throats,
$\unicode[STIX]{x1D646}$
, is updated. Through repeating the above procedure until one meniscus reaches the outlet, we can simulate the quasi-static fluid invasion processes in porous media. More details of the numerical implementation can be found in the existing works (Holtzman & Segre Reference Holtzman and Segre2015; Holtzman Reference Holtzman2016).
Since the numerical method is developed for the 2-D system (figure 2), the fluid advancements with three-dimensional (3-D) nature and the related mechanisms of corner flow (Dong & Chatzis Reference Dong and Chatzis1995; Weislogel & Lichter Reference Weislogel and Lichter1998), snap-off (Roof Reference Roof1970) and spreading of thin wetting films (Levaché & Bartolo Reference Levaché and Bartolo2014; Odier et al.
Reference Odier, Levaché, Santanach-Carreras and Bartolo2017) cannot be considered in numerical simulations. These important mechanisms would dominate the displacement in the case of strong imbibition such as
$\unicode[STIX]{x1D703}<45^{\circ }$
(Concus & Finn Reference Concus and Finn1969) or in the case of relatively high flow rate (Zhao, MacMinn & Juanes Reference Zhao, MacMinn and Juanes2016; Hu et al.
Reference Hu, Wan, Yang, Chen and Tokunaga2018a
). However, for the cases of quasi-static fluid displacement in the range of
$\unicode[STIX]{x1D703}>45^{\circ }$
considered in this work, the numerical model based on Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) is adequate for capturing the fluid invasion processes. The advantage of the employed numerical model is its high computational efficiency compared with other direct simulation techniques, such as the computational fluid dynamics (CFD) method (Ferrari et al.
Reference Ferrari, Jimenez-Martinez, Borgne, Méheust and Lunati2015; Hu et al.
Reference Hu, Wan, Kim and Tokunaga2017a
) and lattice Boltzmann method (LBM) (Liu et al.
Reference Liu, Ju, Wang, Xi and Zhang2015a
). For CFD and LBM, it is extremely expensive, requiring parallel computing lasting for several weeks with 24–95 processors in clusters (Raeini, Bijeljic & Blunt Reference Raeini, Bijeljic and Blunt2015), for simulating such quasi-static fluid displacement. The high efficiency of the numerical method used in this work enables us to systematically investigate the fluid displacement over a wide range of parameters that are difficult to achieve in experiments.
3 Results and discussion
3.1 Impact of disorder on displacement patterns for weak imbibition
Figure 3 presents the observed and simulated displacement patterns in the case of weak imbibition (
$\unicode[STIX]{x1D703}=67^{\circ }$
) from the initial state to the breakthrough time. The experimental results (figure 3
a,c,e,g) demonstrate that increasing disorder
$\unicode[STIX]{x1D706}$
destabilizes the immiscible two-phase flow. For uniform porous media (figure 3
g), a compact displacement is observed, but the front gradually becomes unstable as
$\unicode[STIX]{x1D706}$
increases. The destabilizing effect of the disorder can also be confirmed in terms of the variations of the invading fluid saturation at the breakthrough time
$S_{br}$
(figure 4
a) and the specific fluid–fluid interface length
$l_{nw}$
(figure 4
b), where
$l_{nw}$
is defined by the length of displacement front divided by the average pore throat, excluding the trapped regions behind the front. As shown in figure 4, the measured saturation
$S_{br}$
decreases with the disorder
$\unicode[STIX]{x1D706}$
, while the measured specific interface length
$l_{nw}$
increases with
$\unicode[STIX]{x1D706}$
.
Figure 3. The observed (a,c,e,g) and simulated (b,d,f,h) invasion morphologies at the time of breakthrough when the invading fluid reaches the outlet of the microfluidic flow cell. The disorder
$\unicode[STIX]{x1D706}$
increases from the bottom to top rows, with
$\unicode[STIX]{x1D706}=0$
(g,h),
$\unicode[STIX]{x1D706}=0.1$
(e,f),
$\unicode[STIX]{x1D706}=0.2$
(c,d) and
$\unicode[STIX]{x1D706}=0.3$
(a,b). The red colour is the invading phase and the blue colour is the defending phase. To show the evolution of the displacement front, the dark red represents the initial stage whereas the light red indicates the late time. There are 987 posts in the microfluidics used in the experiments and simulations. We isolate the effect of porosity
$\unicode[STIX]{x1D719}$
, and the values of porosity for all microfluidics are nearly the same, i.e.
$\unicode[STIX]{x1D719}=0.55$
.
The simulated invasion morphologies, shown in figure 3(b,d,f,h), are generally consistent with the experimental results. Inspection of the local fluid distributions in the marked regions (figure 3
a,c,e,g) shows that the defending phase is trapped within single or multiple pores. Although our numerical method is able to capture trapping in fluid displacement (see figure 9
d,e in § 3.3.3), the observed trapping behaviour in microfluidics is not reproduced via simulations. The discrepancy is mainly attributed to limited fabrication resolution, which induces roughness at the edges of posts of the scale of
$1~\unicode[STIX]{x03BC}\text{m}$
. From the theory of multiphase flow, for weak imbibition (
$\unicode[STIX]{x1D703}=67^{\circ }$
) in uniform porous media under quasi-static conditions, all of the neighbouring menisci would merge and the driving front should be flat, leading to 100 % of the defending fluid being displaced, which can be well captured by the simulation (figure 3
h). For the experiment with
$\unicode[STIX]{x1D706}=0$
(figure 3
g), due to manufacturing precision limits, slight changes in post thickness or radii introduce heterogeneity into a homogenous sample (Borgman et al.
Reference Borgman, Fantinel, Lühder, Goehring and Holtzman2017), which determines the invasion path and finally induces trapping (Geistlinger et al.
Reference Geistlinger, Ataei-Dadavi, Mohammadian and Vogel2015). Therefore, in comparison with the experimental results (figure 4), higher saturation
$S_{br}$
and lower interface length
$l_{nw}$
are obtained in the simulation for
$\unicode[STIX]{x1D706}=0$
and for the relatively uniform porous medium (
$\unicode[STIX]{x1D706}=0.1$
). The relative percentage errors of saturation and interface length are, respectively, 21.3 % and
$-21.6\,\%$
for
$\unicode[STIX]{x1D706}=0$
and 21.8 % and
$-16.9\,\%$
for
$\unicode[STIX]{x1D706}=0.1$
. When
$\unicode[STIX]{x1D706}$
increases up to 0.2 and 0.3, the simulated saturation
$S_{br}$
and interface length
$l_{nw}$
are closer to the measurements, with the relative percentage errors decreasing to 4.8 % and 5.2 % (
$\unicode[STIX]{x1D706}=0.2$
) and 5.2 % and 10.3 % (
$\unicode[STIX]{x1D706}=0.3$
). The reason is that, as previously indicated, for more highly disordered porous media, the disorder dominates the fluid invasion over the roughness of post edges at small scale.
Figure 4. Comparison between experimentally measured and numerically simulated results for (a) the invading fluid saturation
$S_{br}$
and (b) the specific fluid–fluid interface length
$l_{nw}$
at the time of breakthrough for different disorders. For the measurements, the data points (solid circles) indicate the average values of four replicate experiments, and the error bars indicate standard deviations.
The trapping behaviour in disordered porous media can be reproduced with simulation in the case of
$\unicode[STIX]{x1D706}=0.25$
and
$\unicode[STIX]{x1D703}=75^{\circ }$
(see figure 9
d) but cannot be reproduced in the disordered microfluidics used in experiments (
$\unicode[STIX]{x1D706}=0.1$
, 0.2, 0.3 and
$\unicode[STIX]{x1D703}=67^{\circ }$
), which means that systematic comparison of the trapped defending fluids between experiments and simulations under various disorders and wetting conditions needs to be conducted. Nevertheless, comparison between experiments and simulations mentioned above shows that the 2-D numerical simulations can generally reproduce the invasion morphology observed in the microfluidics and can capture the impact of disorder on quasi-static displacement patterns for weak imbibition.
3.2 Link pore-filling events to displacement pattern
The pore-scale fluid displacement directly determines the displacement patterns. To obtain further insight into the disorder effect, using the microscope we randomly select the fields of view in the microfluidics with
$\unicode[STIX]{x1D706}=0$
and
$\unicode[STIX]{x1D706}=0.3$
, and record the pore-scale images for local fluid advancements. As shown in figure 5, we observe that, for uniform porous media, the menisci (labelled with 1 and 2) advance from
$t_{1}$
to
$t_{3}$
and finally merge into a single meniscus at
$t_{4}$
, known as overlap (see figure 2
d), which stabilizes the fluid–fluid interface. In disordered media, however, the menisci (1 and 2) cannot merge, because
$r_{4}$
is very large so that the meniscus 1 touches its edge and finally induces a volume of defending fluid to be trapped within the pore throat. Quantitatively, the observed coalescence of neighbouring menisci (figure 5
a) suggests that the conditions for the occurrence of overlap expressed by (2.4) are satisfied in the uniform porous media with
$\unicode[STIX]{x1D703}=67^{\circ }$
, while this condition is not satisfied in the non-uniform post arrangement in figure 5(b). As the number of overlap events in the porous media increases, the displacement pattern becomes more stable.
Figure 5. The pore-scale images for local fluid displacement with the inverted microscope in the (a) uniform and (b) disordered microfluidics at
$t_{1}$
,
$t_{2}$
,
$t_{3}$
and
$t_{4}$
, where the time interval between
$t_{2}$
and
$t_{1}$
, and between
$t_{3}$
and
$t_{2}$
is
$\unicode[STIX]{x0394}t=t_{3}-t_{2}=t_{2}-t_{1}=5$
min, and that between
$t_{4}$
and
$t_{3}$
is 50 s. The arrow indicates the direction of menisci motion and the red dashed line indicates the fluid–fluid interface. The light colour is the invading phase whereas the dark colour is the defending phase. The solid phase (posts) is also shown with dark colour but can be distinguished by its edge (circumference).
To link the advancements of menisci to the displacement patterns, we calculate the probability of overlap,
$P_{ov}$
, by traversing through all of the post arrangements shown in figure 2(d). The probability,
$P_{ov}$
, is defined as
$P_{ov}=N_{ov}/N$
, where
$N_{ov}$
is the number of post arrangements for which the overlap event occurs ((2.4) is satisfied), and
$N$
is the total number of post arrangements. If
$P_{ov}=1$
, overlap occurs for all of the post arrangements, indicating that the displacement is compact. For
$P_{ov}=0$
, no overlap occurs and burst dominates the displacement, and the pattern is capillary fingering. For the cross-over between compact displacement and capillary fingering,
$0<P_{ov}<1$
. Thus, we relate the probability of overlap,
$P_{ov}$
, to the displacement patterns. Holtzman & Segre (Reference Holtzman and Segre2015) performed scaling analysis of viscous and capillary forces and found that the compact displacement only requires 40 %–50 % of overlap events, which is different from the threshold (
$P_{ov}=1$
) in this work. This is due to the different definitions between the probability and the percentage for overlap events. In the work of Holtzman & Segre (Reference Holtzman and Segre2015), the percentage of overlap events is calculated by the ratio of the number of overlaps to all instability events determined from numerical simulations. However, in this work, the probability of overlap
$P_{ov}$
is determined from geometry analysis (without numerical simulations), which means that all of the post arrangements are considered even if some of them are not occupied by the invading fluid. Under quasi-static conditions where the fluid invasion is significantly affected by the geometry, this treatment is reasonable and the relationship between
$P_{ov}$
and displacement patterns will be evaluated in § 3.3.3. Recently, Wang et al. (Reference Wang, Chauhan, Pereira and Gan2019) introduce an indicator, i.e. capillary index
$I_{c}$
, to represent the collaborative effect of disorder and wettability. Based on the variation of fluid–fluid interface length with the invading fluid saturation, they report that
$I_{c}\geqslant 0.5$
corresponds to the stable displacement and fingering occurs for
$I_{c}\leqslant 0.5$
. Note that trapping of defending fluid is also important during the immiscible displacement and would significantly affect the displacement pattern/efficiency. However, under quasi-static displacement conditions, trapping is closely related to the probability of overlap,
$P_{ov}$
. When the dynamic effect is absent, smaller
$P_{ov}$
indicates a larger number of intermittent advancements of menisci at different locations, resulting in a larger trapped volume of the defending fluid (Holtzman & Segre Reference Holtzman and Segre2015), and then destabilizes the displacement front. Therefore, the trapping is implicitly considered in our work.
By traversing through all of the post arrangements and checking the conditions of basic modes given by (2.1), (2.2) and (2.4), we calculate the variation of
$P_{ov}$
with contact angle
$\unicode[STIX]{x1D703}$
(
$45^{\circ }<\unicode[STIX]{x1D703}<180^{\circ }$
) for the four microfluidics with
$\unicode[STIX]{x1D706}=0$
, 0.1, 0.2 and 0.3. Then, we can obtain the critical contact angle below which
$P_{ov}$
is 1.0 (diamonds in figure 6
a), and another critical contact angle above which
$P_{ov}$
is 0 (squares). Thus, as shown in figure 6(a), three zones of compact displacement (
$P_{ov}=1$
), capillary fingering (
$P_{ov}=0$
) and the cross-over between them (
$0<P_{ov}<1$
) are bounded by the two polylines constructed with these critical contact angles (diamonds and squares). The variation of
$P_{ov}$
with
$\unicode[STIX]{x1D706}$
for the experimental conditions can also be directly obtained when
$\unicode[STIX]{x1D703}$
is fixed as
$67^{\circ }$
.
Figure 6. (a) The domains of capillary fingering, compact displacement and the cross-over between them in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane for four specific flow geometries with
$\unicode[STIX]{x1D706}=0$
, 0.1, 0.2 and 0.3. For each flow geometry, the blue diamonds are the critical contact angle below which the probability of overlap
$P_{ov}$
is 1.0, whereas the blue squares correspond to the critical contact angle above which the probability of overlap is 0. The red circles indicate the experimental conditions. (b) The variation of the probability of overlap
$P_{ov}$
with the disorder
$\unicode[STIX]{x1D706}$
for the contact angle of pore surface of the microfluidics, i.e.
$\unicode[STIX]{x1D703}=67^{\circ }$
. The open red star in (a) and (b) indicates the points that separate the stable displacement and the cross-over when the disorder
$\unicode[STIX]{x1D706}$
increases.
Figure 6(a) elucidates the impact of disorder not only in weak imbibition but also in the full range of drainage. Our experimental conditions in weak imbibition are also presented in figure 6 for a comparison purposes. As
$\unicode[STIX]{x1D706}$
increases from 0 to 0.3, the displacement pattern shifts from compact displacement to the cross-over at the critical disorder (
$\unicode[STIX]{x1D706}=0.151$
). The observed destabilizing effect of disorder can be quantified by the probability
$P_{ov}$
(figure 6
b). We observe that
$P_{ov}=1$
in the zone of compact displacement (
$\unicode[STIX]{x1D706}\leqslant 0.151$
), and that
$P_{ov}$
decreases as
$\unicode[STIX]{x1D706}>0.151$
for the cross-over. Figure 6(a) also shows that, for
$\unicode[STIX]{x1D703}>\unicode[STIX]{x1D703}_{c}$
, increasing disorder would stabilize the displacement patterns, and the underlying mechanism will be discussed in § 3.3.2.
3.3 Phase diagram
Section 3.2 demonstrates that the zones of different displacement patterns in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane are bounded by the contours of
$P_{ov}=1$
and
$P_{ov}=1$
. The boundaries that separate the different flow regimes are calculated based on the geometry analysis for specific flow geometry with a given
$\unicode[STIX]{x1D706}$
(0, 0.1, 0.2 and 0.3). Given the infinite flow geometries for a given
$\unicode[STIX]{x1D706}$
, there are infinite boundaries for displacement flow regimes. In this section, we derive a theoretical model for describing the lower and upper bounds for the compact displacement and capillary fingering in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane.
3.3.1 Probability of overlap for an arbitrary geometry
We first redefine the probability of overlap
$P_{ov}$
for an arbitrary geometry with a given
$\unicode[STIX]{x1D706}$
. Since the radii of posts follow a uniform distribution, i.e.
$r_{i}\sim U[(1-\unicode[STIX]{x1D706})\bar{r},(1+\unicode[STIX]{x1D706})\bar{r}]$
, the probability density function
$f(r)$
for the uniform distribution is given by
$$\begin{eqnarray}f(r)=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{1}{2\unicode[STIX]{x1D706}\bar{r}}}, & r\in [(1-\unicode[STIX]{x1D706})\bar{r},(1+\unicode[STIX]{x1D706})\bar{r}],\\ 0, & \text{otherwise}.\end{array}\right.\end{eqnarray}$$
Based on probability theory, the probability of picking a three-post arrangement (
$r_{1}$
,
$r_{2}$
and
$r_{3}$
) out of the geometry with 987 posts for a given
$\unicode[STIX]{x1D706}$
is
$$\begin{eqnarray}f_{123}(r_{1},r_{2},r_{3})=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{1}{(2\unicode[STIX]{x1D706}\bar{r})^{3}}}, & (r_{1},r_{2},r_{3})\in G,\\ 0, & \text{otherwise},\end{array}\right.\end{eqnarray}$$
where
$G$
is a bounded closed region expressed by
$G=\{(r_{1},r_{2},r_{3});(1-\unicode[STIX]{x1D706})\bar{r}\leqslant r_{i}\leqslant (1-\unicode[STIX]{x1D706})\bar{r},~i=1,2,3.\}$
with its volume being
$(2\unicode[STIX]{x1D706}\bar{r})^{3}$
. According to the condition for the occurrence of overlap (2.4), the probability of overlap for an arbitrary geometry with a given
$\unicode[STIX]{x1D706}$
can be written as
$$\begin{eqnarray}P_{ov}=\iiint _{G\cap G_{ov}}\frac{1}{(2\unicode[STIX]{x1D706}\bar{r})^{3}}\,\text{d}r_{1}\,\text{d}r_{2}\,\text{d}r_{3},\end{eqnarray}$$
where
$G_{ov}$
is the region in which the overlap occurs, depending upon
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D703}$
, expressed by
$$\begin{eqnarray}G_{ov}(\unicode[STIX]{x1D706},\unicode[STIX]{x1D703})=\{(r_{1},r_{2},r_{3}):\unicode[STIX]{x1D702}_{2}(r_{1},r_{2},r_{3};\unicode[STIX]{x1D706},\unicode[STIX]{x1D703})+\unicode[STIX]{x1D702}_{3}(r_{1},r_{2},r_{3};\unicode[STIX]{x1D706},\unicode[STIX]{x1D703})\geqslant 2\unicode[STIX]{x03C0}/3\},\end{eqnarray}$$
where
$r_{1}$
,
$r_{2}$
,
$r_{3}$
and
$\unicode[STIX]{x1D702}_{2}$
,
$\unicode[STIX]{x1D702}_{3}$
are variables in the three-post arrangement shown in figure 2(d).
By (3.3), we calculate
$P_{ov}$
as a function of disorder and contact angle, as presented in figure 7(a). The contours of
$P_{ov}$
are employed to classify the different flow regimes for an arbitrary geometry with a given
$\unicode[STIX]{x1D706}$
. Figure 7(a) shows that the flow regimes of compact displacement and capillary fingering are separated by the boundary curves of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane. The boundary curve
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
is the lower bound of the contact angles below which the probability of overlap
$P_{ov}$
is 1.0 (compact displacement), whereas the curve
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
corresponds to the upper bound of the contact angles above which the probability of overlap is 0 (capillary fingering). The area bounded by the two curves is the cross-over zone from capillary fingering to compact displacement. Thus, a phase diagram of the quasi-static fluid displacement pattern with disorder and contact angle can be straightforwardly established when the boundary curves are determined.
Figure 7. (a) Phase diagram of quasi-static immiscible displacement in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane with
$\unicode[STIX]{x1D706}$
ranging from 0 to
$\unicode[STIX]{x1D706}_{max}=a/\bar{r}-1$
(corresponding to the state of posts overlapping), and with
$\unicode[STIX]{x1D703}$
ranging from
$45^{\circ }$
to
$180^{\circ }$
. The probability of overlap
$P_{ov}$
is presented with red colour for a higher value and with blue colour for a lower value. The probability of overlap
$P_{ov}$
as a function of
$\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D706}$
is used to classify the different displacement patterns, i.e.
$P_{ov}=1$
for stable displacement,
$P_{ov}=0$
for capillary fingering and
$0<P_{ov}<1$
for the cross-over between them. The thick curves of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
are the lower and upper bounds of the cross-over zone from capillary fingering to compact displacement. (b–d) Determination of the curves of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
. Each curve,
$\unicode[STIX]{x1D703}_{i}(\unicode[STIX]{x1D706})$
, corresponding to a post arrangement with constant ratios among the three radii, describes the critical
$\unicode[STIX]{x1D703}$
below which no overlap occurs and above which burst occurs. Based on geometry analysis of sampling points, the curves
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
, respectively, correspond to the post arrangements given in (c) and (d).
3.3.2 Phase diagram
To derive the analytical (or semi-analytical) solutions for the boundary curves of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
, we perform geometry analysis on a larger number of post arrangements according to the conditions of the occurrence of overlap (2.4). We consider
$100\times 135$
intervals in the
$\unicode[STIX]{x1D703}{-}\unicode[STIX]{x1D706}$
space of
$[0,\unicode[STIX]{x1D706}_{max}]\times [45^{\circ },180^{\circ }]$
, and consider
$100\times 100\times 100$
intervals for the bounded closed region
$G$
. Thus, a total of
$1.35\times 10^{10}$
sampling post arrangements are analysed. Figure 7(b) presents a series of curves from geometry analysis. Each curve,
$\unicode[STIX]{x1D703}_{i}(\unicode[STIX]{x1D706})$
, corresponds to the critical contact angle for a post arrangement with constant ratios among the three radii. Below
$\unicode[STIX]{x1D703}_{i}(\unicode[STIX]{x1D706})$
no overlap occurs and above it burst occurs. Thus, as shown in figure 7(b), the left envelope of these curves indicates that overlap unconditionally occurs in its left area, whereas overlap never occurs and burst dominates the displacement in the right side of the right envelope. The two envelopes have the same physical meanings as those of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
shown in figure 7(a). Interestingly, we find that the two envelopes respectively correspond to the two post arrangements given in figures 7(c) and 7(d), which can be used to derive analytical (or semi-analytical) solutions of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
.
For the boundary curve of
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
, the related post arrangement (figure 7
d) shows that the three radii reach their minimum and that the burst occurs for the two menisci. Based on geometry analysis (§ B.1), we derive the analytical solution for
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
as
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{\text{CF}}(\unicode[STIX]{x1D706}):\quad \unicode[STIX]{x1D706}=1-\bar{l}^{-1}\left(\left[\frac{f_{1}(\unicode[STIX]{x1D703})}{f_{2}(\unicode[STIX]{x1D703})}\right]^{2}+\frac{f_{1}(2\unicode[STIX]{x1D703})-3}{f_{2}(\unicode[STIX]{x1D703})}+1\right)^{-1/2},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & f_{1}(\unicode[STIX]{x1D703})=\sqrt{3}\sin \unicode[STIX]{x1D703}-3\cos \unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{2}(\unicode[STIX]{x1D703})=4\cos ^{2}\unicode[STIX]{x1D703}-1, & \displaystyle\end{eqnarray}$$
where
$\bar{l}=2\bar{r}/a$
.
Similarly, the boundary curve of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
can be determined by geometry analysis on the post arrangement in figure 7(c). It denotes that the radius of post 3 reaches its minimum while the other two posts have their maximum radii. The burst occurs for the meniscus connecting posts 2 and 3 and later the two menisci merge into a new one (overlap). Again, based on geometry analysis (§ B.2), we derive the semi-analytical solution for
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
:
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706}):\quad g_{4}\unicode[STIX]{x1D6E9}^{4}+g_{3}\unicode[STIX]{x1D6E9}^{3}+g_{2}\unicode[STIX]{x1D6E9}^{2}+g_{1}\unicode[STIX]{x1D6E9}+g_{0}=0,\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=\frac{\bar{l}\cos \unicode[STIX]{x1D703}\,(1-\unicode[STIX]{x1D706}^{2}\bar{l}^{2})+\sqrt{\unicode[STIX]{x1D706}^{2}\bar{l}^{4}f_{4}(\unicode[STIX]{x1D703})-\bar{l}^{2}[\,f_{4}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D706}^{2}]+1}}{2(1-\unicode[STIX]{x1D706}^{2}\bar{l}^{2}\cos ^{2}\unicode[STIX]{x1D703})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{4}=f_{3}^{2}(2\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{3}=-2f_{3}(2\unicode[STIX]{x1D703})f_{3}(\unicode[STIX]{x1D703})(1+\unicode[STIX]{x1D706})\bar{l}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{2}=\left[f_{3}^{2}(\unicode[STIX]{x1D703})+{\textstyle \frac{1}{2}}\,f_{3}(2\unicode[STIX]{x1D703})\right](1+\unicode[STIX]{x1D706})^{2}\bar{l}^{2}-1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{1}=(1+\unicode[STIX]{x1D706})\bar{l}\left[\cos \unicode[STIX]{x1D703}-{\textstyle \frac{1}{2}}(1+\unicode[STIX]{x1D706})^{2}\bar{l}^{2}f_{3}(\unicode[STIX]{x1D703})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{0}={\textstyle \frac{1}{4}}(1+\unicode[STIX]{x1D706})^{2}\bar{l}^{2}\left[{\textstyle \frac{1}{4}}(1+\unicode[STIX]{x1D706})^{2}\bar{l}^{2}-1\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{3}(\unicode[STIX]{x1D703})=\sqrt{3}\sin \unicode[STIX]{x1D703}+\cos \unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{4}(\unicode[STIX]{x1D703})=\sin ^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D706}^{2}\cos ^{2}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$
Finally, we derive the analytical (or semi-analytical) solutions of the upper and lower bounds (in terms of contact angle
$\unicode[STIX]{x1D703}$
) of the cross-over zone for porous media with arbitrary flow geometry at a given disorder
$\unicode[STIX]{x1D706}$
. As shown in (3.5) and (3.6), only one parameter,
$\bar{l}=2\bar{r}/a$
, is included in the theoretical model.
The phase diagram predicted by the theoretical model is also presented in figure 7(a), which describes the cross-over from capillary fingering to compact displacement in disordered porous media from weak imbibition (
$\unicode[STIX]{x1D703}=45^{\circ }$
) to strong drainage (
$\unicode[STIX]{x1D703}=180^{\circ }$
). The cross-over zone denoted by the contact angle range of
$[\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706}),\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})]$
expands as the porous medium becomes more disordered, with a single point located at the critical contact angle
$\unicode[STIX]{x1D703}_{c}$
for uniform porous media and approaching the maximum of
$[\unicode[STIX]{x1D703}_{CDM},\unicode[STIX]{x1D703}_{CFM}]$
at
$\unicode[STIX]{x1D706}_{max}$
. Here,
$\unicode[STIX]{x1D703}_{CDM}$
is the minimum of
$\unicode[STIX]{x1D703}_{CD}$
;
$\unicode[STIX]{x1D703}_{CFM}$
is the maximum of
$\unicode[STIX]{x1D703}_{CF}$
; and
$\unicode[STIX]{x1D706}_{max}$
is the maximum disorder corresponding to the state of post overlapping. The phase diagram also elucidates that the impact of disorder on the displacement pattern depends on the wettability, as further illustrated in figure 8. From figure 8,
$\unicode[STIX]{x1D703}_{CDM}$
,
$\unicode[STIX]{x1D703}_{c}$
and
$\unicode[STIX]{x1D703}_{CFM}$
vary with the compactness of the porous medium that is represented by
$\bar{l}~(=2\bar{r}/a)$
or porosity
$\unicode[STIX]{x1D719}~(=1-(\unicode[STIX]{x03C0}/(2\sqrt{3}))\bar{l}^{2})$
, and the critical contact angle
$\unicode[STIX]{x1D703}_{c}$
decreases with
$\unicode[STIX]{x1D719}$
or
$\bar{r}$
, consistent with the previous numerical results (Koiller et al.
Reference Koiller, Ji and Robbins1992).
Figure 8. The disorder effect on fluid invasion depends not only on contact angle, but also on the compactness represented by
$\bar{l}$
or porosity
$\unicode[STIX]{x1D719}=1-(\unicode[STIX]{x03C0}/(2\sqrt{3}))\bar{l}^{2}$
. (a) Four zones are defined in the
$\bar{l}$
–
$\unicode[STIX]{x1D703}$
plane. For zone I and zone IV, the displacement patterns are respectively compact and unstable, independent of disorder
$\unicode[STIX]{x1D706}$
. In zone II, increasing the disorder
$\unicode[STIX]{x1D706}$
destabilizes the displacement front, and the mechanism is well described by the pore-scale images recorded with microscopy (figure 5). For zone III, increasing the disorder
$\unicode[STIX]{x1D706}$
stabilizes the displacement front, and the mechanism is shown in panel (b).
Figure 8 further indicates that the wettability-dependent disorder impact on quasi-static displacement can be represented with four zones, i.e. zone I, zone II, zone III and zone IV. For zone I (
$\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{CDM}$
) and zone IV (
$\unicode[STIX]{x1D703}>\unicode[STIX]{x1D703}_{CFM}$
), the patterns are respectively compact and unstable, independent of disorder
$\unicode[STIX]{x1D706}$
. In zone II, increasing the disorder
$\unicode[STIX]{x1D706}$
destabilizes the displacement front, which has been well investigated in the previous studies (Koiller et al.
Reference Koiller, Ji and Robbins1992; Yortsos et al.
Reference Yortsos, Xu and Salin1997; Toussaint et al.
Reference Toussaint, Løvoll, Méheust, Måløy and Schmittbuhl2005; Holtzman & Juanes Reference Holtzman and Juanes2010; Xu et al.
Reference Xu, Ok, Xiao, Neeves and Yin2014; Holtzman Reference Holtzman2016), and the mechanism is well described by the pore-scale images recorded with microscopy (figure 4). For zone III, however, increasing the disorder
$\unicode[STIX]{x1D706}$
stabilizes the displacement front, and the mechanism is given in figure 8(b). In this zone, two competing mechanisms are responsible for the disorder effect. On the one hand, for a given contact angle
$\unicode[STIX]{x1D703}$
(
$\unicode[STIX]{x1D703}_{c}<\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{CFM}$
), burst dominates the fluid displacement, destabilizing the fluid–fluid interface. On the other hand, as the porous medium becomes more disordered, overlap may occur to stabilize the interface (figure 8
b). The latter would gradually become a dominant mechanism as
$\unicode[STIX]{x1D706}$
gradually increases, leading to a stabilizing effect.
Holtzman (Reference Holtzman2016) also investigated the effect of the disorder on displacement patterns. For the microfluidics in Holtzman (Reference Holtzman2016),
$\bar{l}=0.54$
, and the corresponding
$\unicode[STIX]{x1D703}_{c}$
is
$87^{\circ }$
(obtained from figure 8
a). The numerical results in Holtzman (Reference Holtzman2016) showed that, for
$\unicode[STIX]{x1D703}=90^{\circ }$
, i.e. in zone III, increasing disorder would destabilize the displacement front, which is inconsistent with our work. This is attributed to the dynamic effect involved in the study of Holtzman (Reference Holtzman2016) (
$Ca=1.4\times 10^{-5}$
). The dynamic effect would enhance the number of fingers, destabilize the displacement front, and finally enhance trapping, which is amplified when increasing disorder (Holtzman Reference Holtzman2016). In other words, when the capillary force plays a role in multiphase flow, increasing disorder would significantly enhance trapping, which dominates the mechanism shown in figure 8(b), and finally destabilizes the displacement front. This destabilizing effect due to the role of viscous force is not considered in this work. Under quasi-static conditions, however, the effect of disorder and wettability can be well captured by the probability of overlap, which will be evaluated via simulations and experimental results shown in § 3.3.3.
3.3.3 Evaluation
Given that the pore-scale simulations can generally capture the impact of disorder on the quasi-static displacement as discussed in § 3.1, the pore-scale simulations enable us to evaluate the proposed phase diagram of quasi-static fluid displacement with a larger range of disorder and contact angle that cannot be achieved in microfluidic experiments. We perform pore-scale numerical simulations with
$45^{\circ }\leqslant \unicode[STIX]{x1D703}\leqslant 180^{\circ }$
and
$0\leqslant \unicode[STIX]{x1D706}\leqslant 0.4$
, and consider 243 points of [
$\unicode[STIX]{x1D703},\unicode[STIX]{x1D706}$
]. For each point, 10 geometries are generated independently. In total, we simulate 2430 computational cases. We evaluate the phase diagram with invading fluid saturation
$S_{br}$
and specific fluid–fluid interface length
$l_{nw}$
at the time of breakthrough, as presented in figure 9. The dots represent the average values for
$S_{br}$
and
$l_{nw}$
from the 10 simulated cases. From figure 9(a,b) we see that the theoretical model (3.5) and (3.6) well predicts the regimes of capillary fingering (figure 9
e), compact displacement (figure 9
c) and the cross-over zone (figure 9
d) in the
$\unicode[STIX]{x1D703}$
–
$\unicode[STIX]{x1D706}$
plane. Figure 9(d,e) also shows the defending fluid trapped within single or multiple pores in disordered porous media, indicating that the numerical method is able to capture the trapping behaviour in multiphase flow.
Figure 9. Evaluation of the proposed phase diagram with simulated invading fluid saturation
$S_{br}$
(a) and specific fluid–fluid interface length
$l_{nw}$
(b) at the time of breakthrough. The dots represent the average values for
$S_{br}$
and
$l_{nw}$
from 10 simulated cases. These 10 geometries are generated independently at a given disorder
$\unicode[STIX]{x1D706}$
. The theoretical model predicts the cross-over from capillary fingering to compact displacement, and the representative fluid invasion morphologies are shown in (c) compact displacement, (d) cross-over and (e) capillary fingering. In (c–e), to show the evolution of displacement front, the dark red represents the initial stage whereas the light red indicates the late time.
Based on the simulation results in figure 9, we also present the details of disorder impact on
$S_{br}$
and
$l_{nw}$
in zone II (figure 10
a,b) and in zone III (figure 10
c,d). As in figure 9, in figure 10, the data points indicate the average values of 10 simulated cases, and the error bars indicate standard deviations. The boundaries separating different displacement patterns predicted by the theoretical model are also presented with vertical lines. In zone II, increasing
$\unicode[STIX]{x1D706}$
decreases
$S_{br}$
and increases
$l_{nw}$
, thus destabilizing the displacement front ranging from compact to cross-over at the critical disorder,
$\unicode[STIX]{x1D706}_{CD}$
. In zone III, conversely, a stabilizing effect can be seen in which, as
$\unicode[STIX]{x1D706}$
increases, the displacement shifts fingering to cross-over at the critical disorder,
$\unicode[STIX]{x1D706}_{CF}$
, which supports the results in figure 8. Given that our numerical simulation can describe the trapping during quasi-static immiscible displacement, the consistency between simulations and the theoretical model indicates that the trapping can be reasonably considered when using the probability of overlap to relate displacement patterns, as discussed in §§ 3.2 and 3.3.2. Compared with zone II, however, in zone III, we observe that the variations of
$S_{br}$
and
$l_{nw}$
are insensitive to the change in
$\unicode[STIX]{x1D706}$
, and thus we can conclude that the impact of disorder on displacement patterns is less significant for
$\unicode[STIX]{x1D703}_{c}<\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{CFM}$
. Figure 10(e,f) confirms that increasing
$\unicode[STIX]{x1D703}$
decreases
$S_{br}$
and increases
$l_{nw}$
for a porous medium with a given
$\unicode[STIX]{x1D706}$
, thus destabilizing the displacement from compact to cross-over to capillary fingering (Jung et al.
Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Hu et al.
Reference Hu, Wan, Yang, Chen and Tokunaga2018a
; Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018). The two critical contact angles
$\unicode[STIX]{x1D703}_{CD}$
and
$\unicode[STIX]{x1D703}_{CF}$
predicted by the theoretical model well capture the boundaries among these three flow regimes.
Figure 10. Variations of invading fluid saturation
$S_{br}$
(a,c,e) and specific fluid–fluid interface length
$l_{nw}$
(b,d,f) with
$\unicode[STIX]{x1D706}$
in zone II (a,b), i.e.
$\unicode[STIX]{x1D703}_{CDM}<\unicode[STIX]{x1D703}=65^{\circ }<\unicode[STIX]{x1D703}_{c}$
, with
$\unicode[STIX]{x1D706}$
in zone III (c,d), i.e.
$\unicode[STIX]{x1D703}_{c}<\unicode[STIX]{x1D703}=85^{\circ }<\unicode[STIX]{x1D703}_{CFM}$
, and with
$\unicode[STIX]{x1D703}$
for
$45^{\circ }<\unicode[STIX]{x1D703}<180^{\circ }$
at
$\unicode[STIX]{x1D706}=0.2$
(e,f). The data points (squares) indicate the average values of 10 simulated cases, and the error bars indicate standard deviations. The boundaries separating different displacement patterns predicted by the theoretical model are also presented with vertical lines. (a,b) In zone II, increasing
$\unicode[STIX]{x1D706}$
decreases
$S_{br}$
and increases
$l_{nw}$
, thus destabilizing the displacement front ranging from compact to cross-over at
$\unicode[STIX]{x1D706}_{CD}$
. (c,d) In zone III, a stabilizing effect is observed. The variations of
$S_{br}$
and
$l_{nw}$
are insensitive to the change in
$\unicode[STIX]{x1D706}$
, indicating that the impact of disorder on the displacement pattern is less significant. (e,f) Increasing
$\unicode[STIX]{x1D703}$
decreases
$S_{br}$
and increases
$l_{nw}$
at
$\unicode[STIX]{x1D706}=0.2$
, thus destabilizing the displacement front, which shifts from compact to cross-over to capillary fingering. The two critical contact angles
$\unicode[STIX]{x1D703}_{CD}$
and
$\unicode[STIX]{x1D703}_{CF}$
separate these three flow regimes.
Figure 11. Evaluation of phase diagram using the experimental results in this work and the experiments by Jung et al. (Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016). In this work,
$\bar{l}=0.7$
and
$\unicode[STIX]{x1D703}=67^{\circ }$
, and the theoretical model is then presented in (a), which is able to well capture the observed displacement patterns of cross-over (open circles) and compact displacement (solid circles), respectively. (b) In the experiments of Jung et al. (Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016) with porosity
$\unicode[STIX]{x1D719}=0.85$
, the theoretical model predicts well the cross-over from the capillary fingering to compact displacement as
$\unicode[STIX]{x1D703}$
decreases from
$150^{\circ }$
to
$45^{\circ }$
. (c) A microfluidic experiment with
$\unicode[STIX]{x1D719}=0.7$
is also adopted herein. The cross-over from fingering to compact displacement is also captured by the theoretical model. In (b,c), the horizontal bars indicate standard deviations of contact angle measurements.
We also evaluate the phase diagram using experimental results. The experiments in microfluidics with different disorders under different wetting conditions are quite limited. Here, the experimental results in this work and the experiments conducted by Jung et al. (Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016) are employed to evaluate the proposed phase diagram. Given that one parameter (
$\bar{l}=2\bar{r}/a$
) is involved in the theoretical model (
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
and
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
), the parameter
$\bar{l}$
needs to be determined. In this work,
$a=250~\unicode[STIX]{x03BC}\text{m}$
,
$\bar{r}=78.5~\unicode[STIX]{x03BC}\text{m}$
and hence
$\bar{l}=0.628$
. The theoretical model is then presented in figure 11(a), which is able to well capture the observed displacement patterns of cross-over (open circles) and compact displacement (solid circles), respectively. We also present the boundary curves for the four specific geometries (previously shown in figure 6
a) in figure 11(a), demonstrating that the cross-over zone determined via specific geometries is included in the cross-over zone predicted by the theoretical model. This is because the boundary curves of
$\unicode[STIX]{x1D703}_{CDM}$
and
$\unicode[STIX]{x1D703}_{CFM}$
represent the lower and upper bounds for the cross-over zone for an arbitrary geometry with a given disorder. Jung et al. (Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016) conducted slow immiscible displacement in disordered microfluidics to investigate the wettability effect. For the microfluidics with porosity
$\unicode[STIX]{x1D719}=0.85$
,
$r=\bar{r}=16~\unicode[STIX]{x03BC}\text{m}$
and
$21~\unicode[STIX]{x03BC}\text{m}\leqslant d\leqslant 65~\unicode[STIX]{x03BC}\text{m}$
, the parameter
$\bar{l}$
can be estimated with
$\unicode[STIX]{x1D719}=1-(\unicode[STIX]{x03C0}/2\sqrt{3})\bar{l}^{2}=0.4067$
and the disorder
$\unicode[STIX]{x1D706}$
is determined as
$\unicode[STIX]{x1D706}=(d_{max}-d_{min})/(4\bar{r})=0.6875$
. As shown in figure 11(b), the theoretical model predicts well the cross-over from capillary fingering to compact displacement as
$\unicode[STIX]{x1D703}$
decreases from
$150^{\circ }$
to
$45^{\circ }$
. For another microfluidic with
$\unicode[STIX]{x1D719}=0.7$
(figure 11
c),
$r=\bar{r}=16~\unicode[STIX]{x03BC}\text{m}$
and
$14~\unicode[STIX]{x03BC}\text{m}\leqslant d\leqslant 40~\unicode[STIX]{x03BC}\text{m}$
,
$\bar{l}$
is determined as 0.5751 and
$\unicode[STIX]{x1D706}=(d_{max}-d_{min})/(4\bar{r})=0.4063$
. The theoretical model again captures well the cross-over from fingering to compact displacement for a given disorder.
4 Conclusions
We study the competing effect of disorder and wettability on the quasi-static fluid displacement in porous media by combining microfluidic experiments and pore-scale simulations with theoretical analysis. Through microfluidic experiments, we show that increasing disorder destabilizes the invasion morphologies at weak imbibition and the underlying mechanism at the pore scale can be well recorded with the microscope. We develop the numerical method originally proposed by Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) to consider the rectangular domain with constant-flow-rate boundaries to match the experimental conditions. Numerical simulations generally agree with the experimental results, and thus the numerical method can be employed to systematically investigate the quasi-static fluid displacement in porous media with different disorders and contact angles that cannot be achieved in microfluidic experiments.
Through theoretical analysis on the probability of overlap events that stabilize the fluid–fluid interface, we propose a theoretical model to describe the lower and upper bounds of cross-over between capillary fingering and compact displacement in porous media with a given disorder. The phase diagram predicted by the theoretical model not only captures the boundaries that separate the regimes but also identifies four zones to elucidate that the impact of disorder depends on contact angle. In zone II (
$\unicode[STIX]{x1D703}_{CDM}<\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{c}$
), increasing
$\unicode[STIX]{x1D706}$
destabilizes the displacement patterns, and the underlying mechanism is well explained by the pore-scale images recorded with the microscope. For zone III (
$\unicode[STIX]{x1D703}_{c}<\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{CFM}$
), however, we show that a stabilizing effect of disorder plays a role because the non-uniform posts would stabilize the local fluid–fluid interface. We find that this destabilizing effect is less significant than that in zone II. In zone I and zone IV, displacement patterns respectively exhibit compact displacement and capillary fingering, independent of both contact angle and disorder. Finally, the proposed phase diagram is evaluated using the extensive pore-scale simulations, and the experimental results in this work and in the literature.
Our proposed phase diagram, different from the previous works that consider the flow rate, the wetting condition and the fluid properties (Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988; Holtzman & Juanes Reference Holtzman and Juanes2010; Zhang et al.
Reference Zhang, Oostrom, Wietsma, Grate and Warner2011), defines the lower and upper bounds of the cross-over zone between compact displacement and capillary fingering for an arbitrary flow geometry with a given disorder. Thus, our work extends the previous studies to consider the impact of flow geometry. The proposed phase diagram provides the first identification of displacement patterns just using the statistical parameters such as porosity, average particle/pore sizes and disorder. Thus, our work is of practical significance for engineering applications, such as geological carbon sequestration, oil recovery and shale gas production, where the identification of the displacement pattern is critical for controlling the hydraulic properties of the multiphase flow system in order to ensure recovery efficiency. Note that in this work we only consider triangular lattice packing and uniform post size distribution for a fixed porosity and a fixed correlation length. Many geologic media often have an intrinsic spatial correlation and do not follow a uniform pore size distribution (or do not exhibit triangular lattice packing). The system size of the flow geometry is also fixed as
$10~\text{mm}\times 5~\text{mm}$
(including 987 posts) in this work, and the system size may also affect the displacement, which is not considered in this work. Thus, the impacts of the spatial correlation length (Holtzman Reference Holtzman2016; Borgman et al.
Reference Borgman, Darwent, Segre, Goehring and Holtzman2019), the distribution type, the post packing and the system size on the multiphase flow need further investigation to improve the generality of the theoretical model. Another important open question is to extend the proposed phase diagram from weak imbibition to strong imbibition in which corner flow plays a key role in quasi-static displacement (Girardo et al.
Reference Girardo, Cingolani, Chibbaro, Diotallevi, Succi and Pisignano2009; Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018).
Acknowledgements
We acknowledge support from the National Natural Science Foundation of China (nos 51779188, 51579188). The numerical calculations in this paper have been done on the supercomputing system in the Supercomputing Center of Wuhan University. We thank J. Wan and T. Tokunaga for their help in the design of the microfluidic chip and the experimental set-up.
Appendix A. The occurrence of three basic modes
A.1 Burst
As shown in figure 2(a), for a three-post arrangement, the arc touches the posts
$O_{1}$
and
$O_{2}$
at
$B$
and
$C$
, and the radius of curvature is
$R=\overline{AB}=\overline{AC}$
. Note that because
$\angle DBO_{1}=\angle ABE=90^{\circ }$
, we have
$\angle ABO_{1}=\angle DBE=\unicode[STIX]{x1D703}$
and similarly
$\angle ACO_{2}=\unicode[STIX]{x1D703}$
. Applying the law of cosines in
$\triangle AO_{1}B$
and
$\triangle AO_{2}C$
, we have
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}y_{1}^{2}=r_{1}^{2}+R^{2}-2r_{1}R\cos \unicode[STIX]{x1D703},\\ y_{2}^{2}=r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703},\end{array}\right\}\end{eqnarray}$$
where
$R$
is the radius of curvature of the arc. The radius
$R$
increases with
$y_{1}$
and
$y_{2}$
. The radius
$R$
reaches its minimum (
$R_{b}$
) corresponding to the condition that
$y_{1}+y_{2}$
also reaches its minimum. In the post arrangement shown in figure 2(a), the minimum of
$y_{1}+y_{2}$
means that the point
$A$
lies on the line
$O_{1}O_{2}$
, i.e.
$y_{1}+y_{2}=a$
. Substituting
$y_{1}+y_{2}=a$
into (A 1), we can obtain the critical radius of curvature for the burst mode as
$$\begin{eqnarray}\text{burst}:\quad R_{b}=\frac{-b_{1}+\sqrt{b_{1}^{2}-4a_{1}c_{1}}}{2a_{1}},\end{eqnarray}$$
with
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle a_{1}=\frac{2}{a^{2}}(r_{1}-r_{2})^{2}\cos ^{2}\unicode[STIX]{x1D703}-2,\\ \displaystyle b_{1}=2(r_{1}+r_{2})\cos \unicode[STIX]{x1D703}-\frac{2}{a^{2}}(r_{1}-r_{2})^{2}(r_{1}+r_{2})\cos \unicode[STIX]{x1D703},\\ \displaystyle c_{1}=\frac{a^{2}}{2}+\frac{1}{2a^{2}}(r_{1}^{2}-r_{2}^{2})^{2}-r_{1}^{2}-r_{2}^{2},\end{array}\right\}\end{eqnarray}$$
where
$r_{1}$
and
$r_{2}$
are the radii of post 1 and post 2, and
$\unicode[STIX]{x1D703}$
is the invading fluid contact angle.
A.2 Touch
As shown in figure 2(b), again, applying the law of cosines in
$\triangle AO_{1}B$
,
$\triangle AO_{2}C$
,
$\triangle AO_{1}O_{3}$
and
$\triangle AO_{2}O_{3}$
, we have
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}y_{1}^{2}=r_{1}^{2}+R_{t}^{2}-2r_{1}R_{t}\cos \unicode[STIX]{x1D703},\\ y_{2}^{2}=r_{2}^{2}+R_{t}^{2}-2r_{2}R_{t}\cos \unicode[STIX]{x1D703},\\ y_{1}^{2}=a^{2}+(R_{t}+r_{3})^{2}-2a(R_{t}+r_{3})\cos \unicode[STIX]{x1D6FE},\\ y_{2}^{2}=a^{2}+(R_{t}+r_{3})^{2}-2a(R_{t}+r_{3})\cos (60^{\circ }-\unicode[STIX]{x1D6FE}).\end{array}\right\}\end{eqnarray}$$
From (A 4), we can obtain the analytical solution of the critical radius of curvature,
$R_{t}$
, for the touch mode as
$$\begin{eqnarray}\text{touch}:\quad R_{t}=\frac{-b_{2}+\sqrt{b_{2}^{2}-4a_{2}c_{2}}}{2a_{2}},\end{eqnarray}$$
with
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}a_{2}=4r_{3}^{2}+4\cos ^{2}\unicode[STIX]{x1D703}\,(r_{1}^{2}+r_{2}^{2}-4r_{1}r_{3})+4r_{3}\cos \unicode[STIX]{x1D703}\,(r_{2}+r_{1})-3a^{2},\\ b_{2}=2a^{2}(r_{1}\cos \unicode[STIX]{x1D703}+r_{2}\cos \unicode[STIX]{x1D703}-2r_{3})+2r_{1}\cos \unicode[STIX]{x1D703}\,(r_{2}^{2}+r_{3}^{2}-2r_{1}^{2})\\ \qquad +\,2r_{2}\cos \unicode[STIX]{x1D703}(r_{1}^{2}+r_{3}^{2}-2r_{2}^{2})+2r_{3}(2r_{3}^{2}-r_{1}^{2}-r_{2}^{2}),\\ c_{2}=a^{2}(a^{2}-r_{1}^{2}-r_{2}^{2}-r_{3}^{2})+{\textstyle \frac{1}{2}}(r_{1}^{2}-r_{2}^{2})^{2}+{\textstyle \frac{1}{2}}(r_{1}^{2}-r_{3}^{2})^{2}+{\textstyle \frac{1}{2}}(r_{2}^{2}-r_{3}^{2})^{2}.\end{array}\right\}\end{eqnarray}$$
A.3 Overlap
As shown in figure 2(d), the overlap will occur if
$\unicode[STIX]{x1D702}_{2}+\unicode[STIX]{x1D702}_{3}\geqslant \angle O_{1}O_{2}O_{3}$
. Similar to the procedures given in (A 2) and (A 5), applying the law of cosines in
$\triangle AO_{1}B$
and
$\triangle AO_{2}C$
, we have
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}y_{1}^{2}=r_{1}^{2}+R^{2}-2r_{1}R\cos \unicode[STIX]{x1D703},\\ y_{2}^{2}=r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703},\end{array}\right\}\end{eqnarray}$$
and the law of sines in
$\triangle AO_{2}C$
leads to
$$\begin{eqnarray}\frac{y_{2}}{\sin \unicode[STIX]{x1D703}}=\frac{R}{\sin \unicode[STIX]{x1D6FC}_{2}}.\end{eqnarray}$$
By the above two equations, the expression for
$\unicode[STIX]{x1D6FC}_{2}$
can be written as
$$\begin{eqnarray}\cos \unicode[STIX]{x1D6FC}_{2}=\frac{r_{2}-R\cos \unicode[STIX]{x1D703}}{\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}.\end{eqnarray}$$
The expression for
$\unicode[STIX]{x1D6FD}_{2}$
can also be obtained via the law of cosines in
$\triangle AO_{1}O_{2}$
,
$\triangle AO_{1}B$
and
$\triangle AO_{2}C$
,
$$\begin{eqnarray}\cos \unicode[STIX]{x1D6FD}_{2}=\frac{a^{2}+r_{2}^{2}-r_{1}^{2}+2(r_{1}-r_{2})\cos \unicode[STIX]{x1D703}}{2a\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}.\end{eqnarray}$$
By (A 9) and (A 10), the angle of
$\unicode[STIX]{x1D702}_{2}$
and
$\unicode[STIX]{x1D702}_{3}$
can be written as
$$\begin{eqnarray}\unicode[STIX]{x1D702}_{2}=\cos ^{-1}\left(\frac{r_{2}-R\cos \unicode[STIX]{x1D703}}{\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}\right)-\cos ^{-1}\left(\frac{a^{2}+r_{2}^{2}-r_{1}^{2}+2(r_{1}-r_{2})\cos \unicode[STIX]{x1D703}}{2a\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}\right).\end{eqnarray}$$
Similarly, the angle of
$\unicode[STIX]{x1D702}_{2}$
can also be expressed as
$$\begin{eqnarray}\unicode[STIX]{x1D702}_{3}=\cos ^{-1}\left(\frac{r_{2}-R\cos \unicode[STIX]{x1D703}}{\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}\right)-\cos ^{-1}\left(\frac{a^{2}+r_{2}^{2}-r_{3}^{2}+2(r_{3}-r_{2})\cos \unicode[STIX]{x1D703}}{2a\sqrt{r_{2}^{2}+R^{2}-2r_{2}R\cos \unicode[STIX]{x1D703}}}\right).\end{eqnarray}$$
Appendix B. The derivation of the theoretical model
B.1 The derivation of
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
For
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
, as shown in figure 7(d), the three radii reach their minimum, i.e.
$r_{1}=r_{2}=r_{3}=(1-\unicode[STIX]{x1D706})\bar{r}$
, and the burst occurs for the two menisci. Based on (A 9), with
$\angle AO_{2}C=60^{\circ }$
, we have
$$\begin{eqnarray}\cos 60^{\circ }=\frac{(1-\unicode[STIX]{x1D706})\bar{r}-R_{b}\cos \unicode[STIX]{x1D703}}{\sqrt{R_{b}^{2}+(1-\unicode[STIX]{x1D706})^{2}\bar{r}^{2}-2(1-\unicode[STIX]{x1D706})\bar{r}R_{b}\cos \unicode[STIX]{x1D703}}}.\end{eqnarray}$$
From (B 1),
$R_{b}$
can be solved as
$$\begin{eqnarray}R_{b}=\frac{\sqrt{3}\sin \unicode[STIX]{x1D703}-3\cos \unicode[STIX]{x1D703}}{1-4\cos ^{2}\unicode[STIX]{x1D703}}(1-\unicode[STIX]{x1D706})\bar{r}.\end{eqnarray}$$
For the burst mode, the critical radius of curvature is expressed by (A 2), and substituting
$r_{1}=r_{2}=r_{3}=(1-\unicode[STIX]{x1D706})\bar{r}$
into (A 2) leads to
$$\begin{eqnarray}R_{b}={\textstyle \frac{1}{2}}\sqrt{a^{2}-[2(1-\unicode[STIX]{x1D706})\bar{r}]^{2}\sin ^{2}\unicode[STIX]{x1D703}}+(1-\unicode[STIX]{x1D706})\bar{r}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$
Combining (B 2) with (B 3), we derive the analytical solution for
$\unicode[STIX]{x1D703}_{CF}(\unicode[STIX]{x1D706})$
, as given in (3.5).
B.2 The derivation of
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
For
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
, as shown in figure 7(c), the radius of post 3 reaches its minimum, i.e.
$r_{3}=(1-\unicode[STIX]{x1D706})\bar{r}$
, while the other two posts reach their maximum radii, i.e.
$r_{1}=r_{2}=(1+\unicode[STIX]{x1D706})\bar{r}$
. The burst occurs for the meniscus connecting posts 2 and 3 and later the two menisci merge into a single one (overlap). Again, for the burst, the critical radius of curvature is expressed by (A 2), and substituting
$r_{3}=(1-\unicode[STIX]{x1D706})\bar{r}$
and
$r_{2}=(1+\unicode[STIX]{x1D706})\bar{r}$
into (A 2) leads to
$$\begin{eqnarray}R_{b}=\frac{2\bar{r}\cos \unicode[STIX]{x1D703}\,(a^{2}-4\unicode[STIX]{x1D706}^{2}\bar{r}^{2})+a\sqrt{\left(\begin{array}{@{}l@{}}16\unicode[STIX]{x1D706}^{2}\bar{r}^{4}[\sin ^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D706}^{2}\cos ^{2}\unicode[STIX]{x1D703}]\\ \quad -\,4a^{2}\bar{r}^{2}(\sin ^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D706}^{2}+\unicode[STIX]{x1D706}^{2}\cos ^{2}\unicode[STIX]{x1D703})+a^{4}\end{array}\right)}}{2(a^{2}-4\unicode[STIX]{x1D706}^{2}\bar{r}^{2}\cos ^{2}\unicode[STIX]{x1D703})}.\end{eqnarray}$$
From the geometry shown in figure 7(c), we have
$$\begin{eqnarray}2\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}=120^{\circ },\quad \cos \unicode[STIX]{x1D6FD}=\frac{a}{2y},\quad \sin \unicode[STIX]{x1D6FC}=\frac{R_{b}\sin \unicode[STIX]{x1D703}}{y}.\end{eqnarray}$$
By (A 9),
$\cos \unicode[STIX]{x1D6FC}$
is written as
$$\begin{eqnarray}\cos \unicode[STIX]{x1D6FC}=\frac{(1+\unicode[STIX]{x1D706})\bar{r}-R_{b}\cos \unicode[STIX]{x1D703}}{y}.\end{eqnarray}$$
Using trigonometric functions, we have
$$\begin{eqnarray}\displaystyle & \sin 2\unicode[STIX]{x1D6FC}=2\sin \unicode[STIX]{x1D6FC}\cos \unicode[STIX]{x1D6FC}=2R_{b}\sin \unicode[STIX]{x1D703}[(1+\unicode[STIX]{x1D706})\bar{r}-R_{b}\cos \unicode[STIX]{x1D703}]/y^{2},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \cos 2\unicode[STIX]{x1D6FC}=\cos ^{2}\unicode[STIX]{x1D6FC}-\sin ^{2}\unicode[STIX]{x1D6FC}=\frac{R_{b}^{2}(\cos ^{2}\unicode[STIX]{x1D703}-\sin ^{2}\unicode[STIX]{x1D703})-2R_{b}(1+\unicode[STIX]{x1D706})\bar{r}\cos \unicode[STIX]{x1D703}+(1+\unicode[STIX]{x1D706})^{2}\bar{r}^{2}}{2y^{2}},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \cos \unicode[STIX]{x1D6FD}=\cos (2\unicode[STIX]{x1D6FC}-120^{\circ })=-\frac{1}{2}\cos 2\unicode[STIX]{x1D6FC}+\frac{\sqrt{3}}{2}\sin 2\unicode[STIX]{x1D6FC}.\qquad & \displaystyle\end{eqnarray}$$
Substituting (B 5), (B 7) and (B 8) into (B 9) leads to
$$\begin{eqnarray}ay=-(\sqrt{3}\sin 2\unicode[STIX]{x1D703}+\cos 2\unicode[STIX]{x1D703})R_{b}^{2}+2(\sqrt{3}\sin \unicode[STIX]{x1D703}+\cos \unicode[STIX]{x1D703})(1+\unicode[STIX]{x1D706})\bar{r}R_{b}-(1+\unicode[STIX]{x1D706})^{2}\bar{r}^{2}.\end{eqnarray}$$
Applying the law of cosines in
$\triangle AO_{2}C$
provides the relation between
$R_{b}$
and
$y$
,
$$\begin{eqnarray}y^{2}=(1+\unicode[STIX]{x1D706})^{2}\bar{r}^{2}+R_{b}^{2}-2(1+\unicode[STIX]{x1D706})\bar{r}R_{b}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$
Finally, the boundary curve
$\unicode[STIX]{x1D703}_{CD}(\unicode[STIX]{x1D706})$
can be determined via the combination of (B 4), (B 10) and (B 11), as rearranged in (3.6).
