1 Introduction
A beautiful array of flow patterns arises when a low-viscosity fluid displaces a more viscous fluid in a porous medium. The problem has been extensively examined through laboratory experiments, as well as numerical simulations and theoretical models (Saffman & Taylor Reference Saffman and Taylor1958; Paterson Reference Paterson1981; Tryggvason & Aref Reference Tryggvason and Aref1983; Chen & Wilkinson Reference Chen and Wilkinson1985; Kadanoff Reference Kadanoff1985; Måløy, Feder & Jøssang Reference Måløy, Feder and Jøssang1985; Nittmann, Daccord & Stanley Reference Nittmann, Daccord and Stanley1985; Bensimon et al.
Reference Bensimon, Kadanoff, Liang, Shraiman and Tang1986; Chen Reference Chen1987; Homsy Reference Homsy1987; Arnéodo et al.
Reference Arnéodo, Couder, Grasseau, Hakim and Rabaud1989; Fernández et al.
Reference Fernández, Albarrán, Fernandez and Albarran1990; Li et al.
Reference Li, Lowengrub, Fontana and Palffy-Muhoray2009; Bischofberger, Ramachandran & Nagel Reference Bischofberger, Ramachandran and Nagel2015). The dynamics of such displacement can be characterized by two dimensionless groups: the ratio of viscous to capillary forces, or the capillary number (
$Ca$
), and the ratio of defending to invading fluid viscosities, or viscosity contrast (
$M$
). For high
$Ca$
, the resulting displacement patterns are reminiscent of diffusion-limited aggregation (Witten, Sander & Sander Reference Witten, Sander and Sander1981; Niemeyer, Pietronero & Wiesmann Reference Niemeyer, Pietronero and Wiesmann1984; Daccord, Nittmann & Stanley Reference Daccord, Nittmann and Stanley1986; Meakin, Tolman & Blumen Reference Meakin, Tolman and Blumen1989; Conti & Marconi Reference Conti and Marconi2010). For low
$Ca$
, the displacement dynamics becomes more intricate, and the emerging patterns display a strong dependence on the pore geometry (Lenormand, Zarcone & Sarr Reference Lenormand, Zarcone and Sarr1983; Lenormand & Zarcone Reference Lenormand and Zarcone1985; Lenormand, Touboul & Zarcone Reference Lenormand, Touboul and Zarcone1988; Fernandez, Rangel & Rivero Reference Fernandez, Rangel and Rivero1991; Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Furuberg, Måløy & Feder Reference Furuberg, Måløy and Feder1996; Ferer et al.
Reference Ferer, Ji, Bromhal, Cook, Ahmadi and Smith2004; Toussaint et al.
Reference Toussaint, Løvoll, Méheust, Måløy and Schmittbuhl2005; Holtzman, Szulczewski & Juanes Reference Holtzman, Szulczewski and Juanes2012) and the wettability of the medium – that is, the chemical affinity of the solid for each fluid (Stokes et al.
Reference Stokes, Weitz, Gollub, Dougherty, Robbins, Chaikin and Lindsay1986; Trojer, Szulczewski & Juanes Reference Trojer, Szulczewski and Juanes2015; Zhao, MacMinn & Juanes Reference Zhao, MacMinn and Juanes2016; Odier et al.
Reference Odier, Levaché, Santanach-Carreras and Bartolo2017). In particular, an intermittent injection-pressure signal emerges in the limit of low
$Ca$
(Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Furuberg et al.
Reference Furuberg, Måløy and Feder1996). Given that in most practical applications visualization of the flow in porous media is not possible, the pressure signal is often the only source of information. Surprisingly, no modelling approach to date has been able to capture the injection-pressure signal across different
$Ca$
and pore wettabilities. Here, we develop a new pore-network model that fills this gap, and we use it to explore the transition from viscous-dominated to capillary-dominated flow regimes by examining the connections among fluid morphology and pressure signal.
Pore-network models of flow in porous media can be broadly classified into two groups: quasi-static and dynamic models (Blunt Reference Blunt2001; Meakin & Tartakovsky Reference Meakin and Tartakovsky2009; Joekar-Niasar & Hassanizadeh Reference Joekar-Niasar and Hassanizadeh2012). Quasi-static models neglect viscous effects and advance the invading fluid through either invasion-percolation (Chandler et al.
Reference Chandler, Koplik, Lerman and Willemsen1982; Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988) or event-based algorithms (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990). Although a quasi-static approach can be effective in reproducing experimental invasion patterns at low
$Ca$
(Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018), it is unable to capture the temporal evolution of the injection-pressure signal. Dynamic network models approximate the flow channels with a network of interconnected capillary tubes. Viscous pressure drops are calculated by assuming fully developed viscous flow within each tube. Local capillary pressures within the network are calculated from either the interface position within pore throats (Aker et al.
Reference Aker, Måløy, Hansen and Batrouni1998b
; Gjennestad et al.
Reference Gjennestad, Vassvik, Kjelstrup and Hansen2018) or through mass balance of the two phases in pore bodies (Al-Gharbi & Blunt Reference Al-Gharbi and Blunt2005; Joekar-Niasar, Hassanizadeh & Dahle Reference Joekar-Niasar, Hassanizadeh and Dahle2010). Another notable class of models is invasion percolation in a gradient: a percolation model designed to incorporate buoyancy forces (Wilkinson Reference Wilkinson1984; Birovljev et al.
Reference Birovljev, Furuberg, Feder, Jøssang, Måløy and Aharony1991; Frette et al.
Reference Frette, Feder, Jøssang and Meakin1992; Meakin et al.
Reference Meakin, Feder, Frette and Jøssang1992), and then extended to model (linear) pressure gradients (Yortsos, Xu & Salin Reference Yortsos, Xu and Salin1997). None of the studies of invasion percolation in a gradient, however, incorporate any notion of wettability (they all deal exclusively with strong drainage), pore-scale dynamics, or capillary-number-dependent pressure fluctuations.
In fact, most existing pore-network models, both quasi-static and dynamic, are limited to strong drainage (or injection of non-wetting fluid) and do not include wettability-induced cooperative pore filling (Aker et al.
Reference Aker, Måløy, Hansen and Batrouni1998b
; Al-Gharbi & Blunt Reference Al-Gharbi and Blunt2005; Holtzman & Juanes Reference Holtzman and Juanes2010; Joekar-Niasar et al.
Reference Joekar-Niasar, Hassanizadeh and Dahle2010). The only dynamic pore-network model to date that includes cooperative pore-filling events (Holtzman & Segre Reference Holtzman and Segre2015) does so by combining pore-level invasion events of Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) with viscous relaxation through the pore network. This viscous-relaxation assumption is at odds with the physics of interface motion in the capillary-dominated regime and, as a result, this model is unable to capture the injection-pressure signal observed experimentally in the limit of intermediate and low
$Ca$
(Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Furuberg et al.
Reference Furuberg, Måløy and Feder1996; Zhao et al.
Reference Zhao, MacMinn and Juanes2016). We present in § 2 a consistent framework that combines viscous, capillary, and wettability effects in a single dynamic network model that builds a direct analogy between local fluid–fluid interfaces and electric capacitors. Our model reproduces, quantitatively, the fluid–fluid displacement patterns for a wide range of
$Ca$
and wettabilities (§ 3), and points to a surprising and heretofore unrecognized transition in the pressure fluctuations between the low- and high-
$Ca$
flow regimes (§ 4).
2 Moving-capacitor model

Figure 1. (a) Schematic diagram of in-plane and out-of-plane curvatures within the flow cell. Out-of-plane curvature represents the overall affinity of the porous medium to the invading fluid. It is determined by
$\unicode[STIX]{x1D703}$
and is analogous to a battery. In-plane curvature changes as the local interface evolves while pinned to a pore throat, and it is analogous to a capacitor. (b) Evolution of burst, touch, and overlap events. (c) Temporal profiles of the injection pressure bear close resemblance to similar experiments in the drainage regime at low (orange) and high (blue)
$Ca$
(Furuberg et al.
Reference Furuberg, Måløy and Feder1996; Zhao et al.
Reference Zhao, MacMinn and Juanes2016).
Consider a moving fluid–fluid interface in a micromodel (figure 1
a). Neglecting dynamic-contact-angle effects (Hoffman Reference Hoffman1975) for simplicity, the shape of the meniscus between posts is uniquely defined by the combination of Laplace pressure and substrate wettability, defined through a contact angle
$\unicode[STIX]{x1D703}$
at which the interface meets post surfaces (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990). As the interface advances, the Laplace pressure increases until the interface encounters a burst, touch or overlap event, as defined by Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990). The burst event is equivalent to a Haines jump (Haines Reference Haines1930; Berg et al.
Reference Berg, Ott, Klapp, Schwing, Neiteler, Brussee, Makurat, Leu, Enzmann and Schwarz2013), while the touch and overlap events take place when the local interface either touches the nearest opposing post or coalesces with a neighbouring interface, respectively (figure 1
b). If the interface becomes unstable due to burst or touch, a single pore is invaded and two new interfaces appear. In the case of an overlap event, two (in some cases more) pores are filled simultaneously. These pore-level events are an integral part of the model and, indeed, this sensitivity is what permits capturing wettability effects within the model. The events evolve differently at different wettabilities – burst events are most frequent in drainage, while touch and overlap are most frequent in imbibition (or injection of wetting fluid) (Cieplak & Robbins Reference Cieplak and Robbins1990; Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018).
We can explicitly calculate the critical Laplace pressure
$\unicode[STIX]{x0394}p_{crit}$
corresponding to all events from the values of the contact angle, radii and coordinates of the posts (Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018), and thus can use the analogy between electric capacitors and fluid–fluid interfaces in constructing our network model. A capacitor represents the pinning of the fluid–fluid interface at a pore throat, and is active in both drainage and imbibition: the interface moves only when a local depinning threshold (
$\unicode[STIX]{x0394}p_{crit}$
) is reached, and the fluid front moves to restart the pinning–depinning cycle from zero in-plane curvature (figure 1
b). This progression of the in-plane curvature in our model was motivated by the work of Cieplak and Robbins (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) (see also Rabbani et al.
Reference Rabbani, Zhao, Juanes and Shokri2018) and experiments on the progression of the in-plane curvature between the Hele-Shaw cell posts (Jung et al.
Reference Jung, Brinkmann, Seemann, Hiller, Sanchez de La Lama and Herminghaus2016; Lee et al.
Reference Lee, Gupta, Hatton and Doyle2017). This is what allows capturing pressure fluctuations in the limit of low
$Ca$
(figure 1
c). The battery analogy represents the overall affinity of the porous medium to the invading fluid, set by the out-of-plane curvature at the fluid front. The out-of-plane curvature is fixed throughout a single simulation, and determined by the value of the contact angle (given the constant gap between the flow-cell plates): it is positive in drainage and negative in imbibition (figure 1
a). To complete the analogy between an electric circuit and a pore network, one can think of a network of resistors being responsible for viscous effects, capacitors and batteries responsible for capillary effects, and local rules for circuit rearrangements responsible for wettability effects (figure 1
b).
Therefore, the pressure drop across an edge of the network containing a fluid–fluid interface has three components: (i) pressure drop due to viscous dissipation, (ii) Laplace pressure drop due to in-plane curvature of the interface, and (iii) Laplace pressure drop due to out-of-plane curvature of the interface. We calculate the viscous pressure drop assuming Poiseuille flow in a capillary tube, which is analogous to the potential drop across a resistor. The out-of-plane component of the Laplace pressure can be expressed as either a positive or negative pressure jump (
$\unicode[STIX]{x0394}p_{\bot }=-(2\unicode[STIX]{x1D6FE}\cos \unicode[STIX]{x1D703}/h)$
, where
$\unicode[STIX]{x1D6FE}$
is the interfacial tension, and
$h$
is the cell height) depending on the substrate wettability; this is analogous to a battery in an electric circuit. The Laplace pressure due to in-plane curvature of the interface is analogous to a capacitor which allows flow until it reaches the critical pressure (
$\unicode[STIX]{x0394}p_{crit}=\text{min}\{p_{burst},p_{touch},p_{overlap}\}$
). Since we can calculate
$\unicode[STIX]{x0394}p_{crit}$
for all edges at the invading fluid front, we use a linear estimate of the in-plane Laplace pressure drops within our network (
$\unicode[STIX]{x1D6F7}(t)\unicode[STIX]{x0394}p_{crit}$
), where
$\unicode[STIX]{x1D6F7}(t)$
stands for the filling ratio of a given throat. When
$\unicode[STIX]{x1D6F7}(t)\rightarrow 0$
, the in-plane Laplace pressure is negligible. When
$\unicode[STIX]{x1D6F7}(t)\rightarrow 1$
, the throat is nearly full and has a critical in-plane Laplace pressure
$\unicode[STIX]{x0394}p_{crit}$
. This analogy between local interfaces and capacitors allows us to incorporate local changes in Laplace pressure due to filling of pore throats. Once a node in the network reaches its maximal potential, which coincides with its filling capacity, it becomes unstable and the interface advances. We assume that the in-plane and out-of-plane Laplace pressures are decoupled, and this is done to maintain the simplicity of the overall model. With this assumption, one can run the model for either
$h/a\gg 1$
or
$h/a\ll 1$
, where these conditions would result in negligible or dominant contributions of the out-of-plane curvature in the model, respectively.
The topology of the pore network is captured through the incidence matrix
$\unicode[STIX]{x1D63C}$
by examining the adjacency of the pores (Strang Reference Strang2007). We number all pores and adopt the convention that pore connections are oriented in the direction of increasing pore numbers. Rows of
$\unicode[STIX]{x1D63C}$
represent edges, and columns of
$\unicode[STIX]{x1D63C}$
represent nodes of the network. We also make use of the diagonal conductance matrix
$\unicode[STIX]{x1D63E}$
, whose elements are hydraulic conductivities of the network edges. The elements of this matrix can be calculated as
$c=\unicode[STIX]{x03C0}r^{4}/8\unicode[STIX]{x1D707}L$
, assuming fully developed Hagen–Poiseuille flow through a rectangular tube with hydraulic radius
$r$
and length
$L$
, where
$\unicode[STIX]{x1D707}$
is the effective viscosity of the fluid in the channel.
The pressure difference across the network edges can be calculated as
$\boldsymbol{e}=\boldsymbol{b}-\unicode[STIX]{x1D63C}\boldsymbol{p}$
, where
$\boldsymbol{b}$
and
$\boldsymbol{p}$
stand for pressure change due to out-of-plane contribution to Laplace pressure (batteries) and node pressures, respectively. The network flow rates can be calculated from this pressure difference as
$\boldsymbol{q}=\unicode[STIX]{x1D63E}\boldsymbol{e}$
. At the same time, flow rates must obey mass conservation,
$\unicode[STIX]{x1D63C}^{T}\boldsymbol{q}=\boldsymbol{f}$
, where
$\boldsymbol{f}$
stands for flow sources at the nodes. After eliminating
$\boldsymbol{e}$
, the flow through the network without the in-plane contribution to Laplace pressure (capacitors) is obtained through the following system of equations:


We set constant flow boundary conditions at the inlet pores (at the centre of the flow cell) and constant-pressure boundary conditions at the outlet pores (at the edges of the flow cell). We note that
$\unicode[STIX]{x1D63C}\boldsymbol{p}$
can be decomposed into components of nodes with prescribed pressure and all other nodes (
$\unicode[STIX]{x1D63C}\boldsymbol{p}=\unicode[STIX]{x1D63C}_{outer}\boldsymbol{p}_{outer}+\tilde{\unicode[STIX]{x1D63C}}\tilde{\boldsymbol{p}}$
), and therefore (2.1)–(2.2) transform to

The solution to (2.3) provides values of both edge flow rates and node pressures for given boundary conditions.
Finally, we incorporate the pressure drop due to in-plane Laplace pressure (capacitors) within the network. Taking into account the direction of the edges (an array
$\boldsymbol{d}(t)$
consisting of 1 and
$-1$
), the total pressure drop across the network edges can be written as
$\boldsymbol{e}=\tilde{\boldsymbol{b}}-\tilde{\unicode[STIX]{x1D63C}}\tilde{\boldsymbol{p}}-\boldsymbol{d}(t)\unicode[STIX]{x1D6F7}(t)\unicode[STIX]{x0394}p_{crit}$
. In other words, the in-plane Laplace pressure is the product of the filling ratio and the critical pressure from the quasi-static model (Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018). Therefore, the equations governing two-phase flow through the network can be written as

We now discuss the mechanics of the time stepping in our two-phase flow model. After we initialize the interface locations within the network, we use an adaptive forward Euler time stepping to update the filling ratios of the network edges at the interface
$\unicode[STIX]{x1D6F7}(t)$
. We ensure that only a fraction of the edge total volume at the interface flows within the time step (Aker et al.
Reference Aker, Måløy, Hansen and Batrouni1998b
). After every time step, we use
$\unicode[STIX]{x1D6F7}(t)$
to update the conductance matrix
$\unicode[STIX]{x1D63E}(t)$
and resolve the flow through (2.4) with updated pressure drops across the fluid–fluid front.
In the spirit of the fundamental contributions from Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990), our model takes the form of an arrangement of cylindrical posts confined between the plates of a Hele-Shaw cell. The approach is simple enough to lead to universal findings, yet sufficiently complex to have direct relevance to microfluidic geometries, as well as engineered and natural porous media – much like Lenormand’s phase diagram (Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988). By doing so, we demonstrate the ability to reproduce physics – in particular, pressure fluctuations under a wide range of wetting conditions – which, until now, were inaccessible to pore-network modelling. A limitation of the model presented here is that it does not extend to contact angles below
$45^{\circ }$
, where the wetting fluid preferentially wets the corners of the pore geometry at low
$Ca$
and forms film flow at high
$Ca$
(Zhao et al.
Reference Zhao, MacMinn and Juanes2016; Odier et al.
Reference Odier, Levaché, Santanach-Carreras and Bartolo2017).
3 Invasion patterns
We simulate immiscible fluid–fluid displacement by setting a constant injection rate at the centre of the flow cell and zero pressure at the outlets. The invading and defending fluid viscosities are set to
$8.9\times 10^{-4}~\text{Pa}~\text{s}$
and 0.34 Pa s respectively. The post height
$h$
is
$100~\unicode[STIX]{x03BC}\text{m}$
, and interfacial tension
$\unicode[STIX]{x1D6FE}$
is set to
$13\times 10^{-3}~\text{N}~\text{m}^{-1}$
. These parameters as well as the pore geometry are chosen to mimic the experiments of Zhao et al. (Reference Zhao, MacMinn and Juanes2016). The flow cell has an outer diameter of 30 cm. We perform simulations for wetting conditions from strong drainage (
$\unicode[STIX]{x1D703}=160^{\circ }$
) to weak imbibition (
$\unicode[STIX]{x1D703}=46^{\circ }$
). Figure 1(c) shows the pressure profiles for
$\unicode[STIX]{x1D703}=160^{\circ }$
at
$Ca\in \{10^{-3},10^{-7}\}$
, respectively. In the limit of high
$Ca$
, the more viscous defending fluid sustains substantial spatial pressure gradients, and the injection pressure gradually drops as more of the defending fluid is displaced (Zhao et al.
Reference Zhao, MacMinn and Juanes2016). In contrast, in the limit of low
$Ca$
, the pressure field is virtually uniform in each fluid, and the injection pressure exhibits intermittent fluctuations typical of slow capillary-dominated drainage (Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Aker, Måløy & Hansen Reference Aker, Måløy and Hansen1998a
; Knudsen & Hansen Reference Knudsen and Hansen2002; Moebius & Or Reference Moebius and Or2012).

Figure 2. (a) Phase diagram of the invading fluid morphology at breakthrough; (b) fractal dimension, computed by means of the box-counting method; (c) number of fingers per unit area of injected fluid, which exhibits a maximum near
$\unicode[STIX]{x1D703}=90^{\circ }$
; (d) normalized finger width (
$w/a$
) at different
$Ca$
and wettabilities measured at breakthrough. Finger width increases as the posts become more wetting to the invading fluid.
The morphology of the invading fluid at breakthrough can be analysed by means of a binary-image representation of the invasion patterns (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990; Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018) (figure 2
a). We estimate the width and number of fingers in the invading fluid pattern following the protocol outlined in Cieplak & Robbins (Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990) and modified in Primkulov et al. (Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018). The binary image is sliced horizontally and vertically, with each slice containing clusters of invading fluid pixels. We calculate the finger width as the mean size of these clusters. Figure 2(d) shows that the finger width, normalized by the typical pore size, increases as
$\unicode[STIX]{x1D703}\rightarrow 46^{\circ }$
for all
$Ca$
, which is in agreement with experimental observations (Stokes et al.
Reference Stokes, Weitz, Gollub, Dougherty, Robbins, Chaikin and Lindsay1986; Trojer et al.
Reference Trojer, Szulczewski and Juanes2015; Zhao et al.
Reference Zhao, MacMinn and Juanes2016). While figure 2(a) demonstrates that the number of fingers increases with
$Ca$
(Lenormand et al.
Reference Lenormand, Touboul and Zarcone1988; Fernández et al.
Reference Fernández, Albarrán, Fernandez and Albarran1990; Zhao et al.
Reference Zhao, MacMinn and Juanes2016), we observe an unexpected behaviour (figure 2
b): the finger density changes with the substrate wettability, and exhibits a maximum at approximately
$\unicode[STIX]{x1D703}=90^{\circ }$
. This effect is most pronounced for
$10^{-6}<Ca<10^{-3}$
(when viscous and capillary effects are comparable).
We explain the peak in the viscous finger density at
$\unicode[STIX]{x1D703}\approx 90^{\circ }$
in figure 2(b) by considering in-plane and out-of-plane contributions to the Laplace pressure. At a fixed
$Ca$
, the ratio of viscous and capillary forces in the micromodel changes as a function of substrate wettability. The capillary forces have out-of-plane contributions, which are nominally equal to zero when
$\unicode[STIX]{x1D703}=90^{\circ }$
, so the ratio of viscous and capillary forces increases as
$\unicode[STIX]{x1D703}$
changes from
$160^{\circ }$
to
$90^{\circ }$
at fixed
$Ca$
. In addition, when
$\unicode[STIX]{x1D703}$
changes from
$90^{\circ }$
to
$46^{\circ }$
, the cooperative pore-filling mechanisms become dominant and widen the largest fingers, which in turn consume the smaller ones and reduce the number of fingers. The combination of these two effects results in the local maximum in the number of viscous fingers around
$\unicode[STIX]{x1D703}\approx 90^{\circ }$
across different
$Ca$
(figure 2
b).
For a contact angle
$\unicode[STIX]{x1D703}$
near
$160^{\circ }$
(strong drainage) and high values of
$Ca$
(
$10^{-3}$
and
$10^{-4}$
), the invading fluid front advances through viscous fingers with fractal dimension close to 1.71, typical of diffusion-limited aggregation (DLA)-type morphology (Witten et al.
Reference Witten, Sander and Sander1981). As
$Ca$
is reduced to a low value (
$10^{-7}$
), the fractal dimension increases to approximately 1.82, characteristic of invasion percolation (Wilkinson & Willemsen Reference Wilkinson and Willemsen1983) (figure 2
b). This increasing trend in fractal dimension is consistent with the decrease in finger density (figure 2
c) and the increase in finger width (figure 2
d).
As the contact angle approaches
$46^{\circ }$
, cooperative pore filling becomes the dominant flow mechanism at all values of
$Ca$
. This flow regime results in the compact displacement of the defending fluid, and thus the fractal dimension increases, approaching a value of 2 at low
$Ca$
, indicative of stable displacement.
4 Pressure signature

Figure 3. (a,b) Temporal evolution of the injection pressure at
$Ca=10^{-3}$
and
$Ca=10^{-7}$
respectively. At high
$Ca$
, the injection pressure decreases as the viscous fingers approach the outer boundary of the flow cell. At low
$Ca$
, the injection pressure is dominated by Laplace pressure fluctuations at the interface. We use wavelet decomposition (Cai Reference Cai2002; Sygouni, Tsakiroglou & Payatakes Reference Sygouni, Tsakiroglou and Payatakes2006, Reference Sygouni, Tsakiroglou and Payatakes2007) to split the pressure signal (
$Ca=10^{-5}$
and
$\unicode[STIX]{x1D703}=160^{\circ }$
here) into its (c) global trend and (d) cyclic component. (e) The standard deviation of the pressure fluctuations point at two different regimes. At low
$Ca$
, pressure fluctuations are dominated by stick–slip changes in Laplace pressure. At high
$Ca$
, pressure fluctuations are dominated by changes in the effective hydraulic conductance of dominant flow channels.
The fundamental difference in the fluid–fluid displacement process between low and high
$Ca$
is reflected in the temporal injection-pressure signals (figure 3). When the capillary number is relatively high (
$Ca=10^{-3}$
), viscous forces dominate, and the injection pressure decreases with time for all substrate wettabilities (Zhao et al.
Reference Zhao, MacMinn and Juanes2016) (figure 3
a). Here, most of the pressure drop takes place in the more viscous defending fluid. Consequently, as more of the defending fluid is displaced, the pressure required to maintain the prescribed injection flow rate decreases. In contrast, at
$Ca=10^{-7}$
, viscous dissipation is negligible, and the injection pressure is determined by the sum of outlet and Laplace pressures. As a result, the injection pressure fluctuates in a stick–slip manner around a mean value (figure 3
b), as has been documented in slow drainage experiments (Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Furuberg et al.
Reference Furuberg, Måløy and Feder1996; Moebius & Or Reference Moebius and Or2012). The pressure signals in figure 3(b) highlight the roles that in-plane and out-of-plane curvatures play in our model. Out-of-plane curvature plays the role of batteries, and thus provides additional resistance/drive (in drainage/imbibition, respectively) to the flow at the interface. The magnitude of the pressure drop/rise at the batteries is a function of wettability, which explains why the mean value of the injection-pressure signal also varies with wettability (figure 3
b). The in-plane curvature plays the role of capacitors. As the invading fluid is injected, the in-plane component of Laplace pressure grows at the interface until the meniscus near the pore with lowest critical entry pressure becomes unstable due to burst, touch or overlap. This results in the rapid advance of the local interface, which pressurizes the defending fluid ahead. This overpressure then dissipates (see movie S1 in the supplementary materials available at https://doi.org/10.1017/jfm.2019.554). The critical pressures of touch and overlap are always smaller than the critical pressures of burst events (Cieplak & Robbins Reference Cieplak and Robbins1988, Reference Cieplak and Robbins1990; Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018), so the magnitude of the pressure fluctuations decreases as the substrate becomes more wetting to the invading fluid (figure 3
b).
To gain further insight into the difference in the pressure signature between low and high
$Ca$
, we decompose the injection-pressure signal into its global trend and fluctuating components with block James–Stein wavelet decomposition (Cai Reference Cai2002) (see figure 3
c,d). We compute the standard deviation of the fluctuating component of the pressure signal for both drainage and imbibition conditions (
$\unicode[STIX]{x1D703}=160^{\circ }$
and
$46^{\circ }$
, respectively) for a wide range of
$Ca$
, and find that it exhibits two distinct regimes (figure 3
e). At low
$Ca$
, pressure fluctuations are controlled by the stick-slip-type changes in local Laplace pressures. In contrast, at high
$Ca$
, pressure fluctuations are controlled by changes in the effective hydraulic conductance of the dominant flow channels. In the limit of high
$Ca$
, the Laplace pressure drop is negligible in comparison with the viscous pressure gradient, but the dominant flow channels are rearranged slightly as the fingers grow (see movie S2 in the supplementary materials). Since the pore geometry has a heterogeneous distribution of throat sizes, shifts in the dominant flow channels result in viscosity-driven pressure fluctuations at high
$Ca$
.

Figure 4. (a) Pore-scale perspective for the scaling of pressure fluctuations. The diagram shows a typical pore being invaded. The characteristic distance between the pore centres is
$l$
(red line), the post height is
$h$
, and a characteristic throat size is
$a$
. (b) Typical configurations of the fluid–fluid interface in drainage and imbibition. Burst events are prevalent in drainage and the typical radius of out-of-plane curvature is of order
$a$
. Overlap events are prevalent in imbibition and the typical radius of out-of-plane curvature is an order of magnitude greater than
$a$
.
Scaling arguments support the findings from the model simulations. Let us take a pore-scale perspective (see figure 4). Invading a single pore involves overcoming a capillary pressure and pushing defending fluid out through a throat of width
$a$
and height
$h$
at a speed proportional to the injection rate. The capillary pressure is
$p_{cap}\approx \unicode[STIX]{x1D6FE}((1/h)+(1/af(\unicode[STIX]{x1D703})))$
, where
$f(\unicode[STIX]{x1D703})$
is a wettability-dependent function that takes a value
${\sim}1$
near drainage and
${\sim}10$
near strong imbibition (figure 4
b). Taking variations of
$p_{cap}$
with
$a$
yields

The characteristic flow velocity through a typical throat is
$u=(k(a,h)/\unicode[STIX]{x1D707})(p_{visc}/l)$
, where
$k(a,h)=R_{h}^{2}/8$
is the rectangular channel permeability and
$R_{h}=ah/2(a+h)$
the hydraulic radius. Thus the viscous pressure is
$p_{visc}\sim (32(a+h)^{2}/a^{2}h^{2})\unicode[STIX]{x1D707}ul=32\unicode[STIX]{x1D707}ul/h^{2}(1+h/a)^{2}$
. Taking variations of
$p_{visc}$
with
$a$
yields

The magnitude of the total characteristic pressure fluctuation is
$\unicode[STIX]{x1D6FF}p_{cap}+\unicode[STIX]{x1D6FF}p_{visc}$
, and its two components are comparable when
$\unicode[STIX]{x1D6FF}p_{visc}/\unicode[STIX]{x1D6FF}p_{cap}\sim 1$
. Using (4.1) and (4.2),

which implies a crossover
$Ca$
,

between flow-rate-independent and flow-rate-dependent pressure fluctuations (figure 3
e). The above argument suggests two interesting implications. First, one can potentially infer the characteristic pore size of the material from the fluctuations of the pressure signal in both viscous-dominated and capillary-dominated flow regimes. This is especially useful when visualization of the flow in pore space is not possible, which is the case in most porous materials. Second, the characteristic
$h$
,
$a$
and
$l$
used in this study yield
$Ca^{\ast }\approx 10^{-3}/f(\unicode[STIX]{x1D703})$
, which reduces to
$Ca^{\ast }\sim 10^{-3}$
for drainage and
$Ca^{\ast }\sim 10^{-4}$
for imbibition, in agreement with the data in figure 3(e). This means that one should expect the transition from capillary-dominated to viscous-dominated flow regimes at different
$Ca^{\ast }$
in drainage and imbibition. The order of magnitude of
$f(\unicode[STIX]{x1D703})$
was obtained by calculating
$\unicode[STIX]{x0394}p_{crit}$
for all pore throats at
$\unicode[STIX]{x1D703}\in \{46^{\circ },160^{\circ }\}$
with the quasi-static model (Primkulov et al.
Reference Primkulov, Talman, Khaleghi, Rangriz Shokri, Chalaturnyk, Zhao, MacMinn and Juanes2018) and taking an average of
$f(\unicode[STIX]{x1D703})=\unicode[STIX]{x1D6FE}/a\unicode[STIX]{x0394}p_{crit}$
for each contact angle. Finally, the viscous pressure fluctuation component scales as
$\unicode[STIX]{x1D6FF}p_{visc}\sim \unicode[STIX]{x1D707}u$
, which is equivalent to
$\unicode[STIX]{x1D6FF}p_{visc}\sim Ca$
when interfacial tension is kept constant. This explains the slope of the viscous-dominated portion of the graph in figure 3(e).
5 Conclusion
Overall, our moving-capacitor network model provides new fundamental insights into the dynamics of immiscible fluid–fluid displacement in porous media for a wide range of
$Ca$
and wettabilities. The model completes the picture of the displacement by covering both high and low
$Ca$
, which allows one, for the first time, to reproduce experimental observations of invading fluid patterns (Zhao et al.
Reference Zhao, MacMinn and Juanes2016), injection pressure and front velocity in drainage (Måløy et al.
Reference Måløy, Furuberg, Feder and Jossang1992; Furuberg et al.
Reference Furuberg, Måløy and Feder1996; Moebius & Or Reference Moebius and Or2012) and imbibition. Our observations and scaling arguments on the transition from a viscous-dominated to a capillary-dominated flow regime suggest that it is possible to infer the character of the multiphase-flow displacement purely from the injection-pressure signal. This poses an exciting prospect for detailed experiments.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2019.554.