1. Introduction
Intricate dendritic patterns emerge when a less viscous fluid displaces a more viscous one within a porous medium. Such fluid–fluid displacement is inherently unstable due to viscosity contrast, giving rise to the well-known Saffman–Taylor, or viscous fingering, instability (Saffman Reference Saffman1986; Homsy Reference Homsy1987; Lenormand Reference Lenormand1990). A crucial practical question is how efficiently the invading fluid sweeps the domain upon first reaching the outer boundary (i.e. at breakthrough), and whether the sweep efficiency continues to improve thereafter. A fundamental understanding of this process can inform more effective hydrocarbon recovery, soil remediation, in situ mineral extraction and chemical reactor bed cleanup. Yet, the post-breakthrough evolution of viscous fingering remains underexplored and poorly understood.
While the evolution of invading-fluid saturation after breakthrough has been reported in numerous studies of viscous fingering (Chuoke et al. Reference Chuoke, van Meurs and van der Poel1959; Odier et al. Reference Odier, Levaché, Santanach-Carreras and Bartolo2017; Suo et al. Reference Suo, Liu and Gan2020; Mora et al. Reference Mora, Morra, Yuen and Juanes2021; Salmo et al. Reference Salmo, Sorbie, Skauge and Alzaabi2022), considerably less attention has been paid to the accompanying evolution of the displacement morphology. In particular, post-breakthrough pattern development has not been examined within the framework of Lenormand’s seminal phase diagram (Lenormand & Bories Reference Lenormand and Bories1980; Lenormand et al. Reference Lenormand, Zarcone and Sarr1983; Lenormand Reference Lenormand1990; Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021), which remains the standard classification of immiscible displacement regimes in porous media. As a result, it remains unclear whether the regime identified at breakthrough persists thereafter or evolves as the invading fluid continues to occupy the pore space.
Post-breakthrough evolution of viscous fingering is also linked to a broader class of physical problems. Paterson (Reference Paterson1984) noted the mathematical analogy between viscous fingering and diffusion-limited aggregation (Meakin et al. Reference Meakin, Tolman and Blumen1989); in both cases, pattern growth within a homogeneous domain is governed by the Laplace equation. This analogy extends further: diffusion-limited aggregation has been used to model fire-front propagation (Conti & Marconi Reference Conti and Marconi2010) and dielectric breakdown in insulators (Niemeyer et al. Reference Niemeyer, Pietronero and Wiesmann1984). The closest parallel is dissolution fingering in rocks (Zukowski et al. Reference Zukowski, Cornelissen, Osselin, Douady and Szymczak2024, Reference Zukowski, Magni, Osselin, Dutka, Cooper, Ladd and Szymczak2025), where the pressure distribution is nearly identical to that in viscous fingering, so displacement and dissolution patterns are expected to evolve similarly before and after breakthrough.
In this study, we examine the post-breakthrough dynamics of viscous fingers through experiments and simulations conducted in a micromodel porous medium (figure 1). We demonstrate that, after breakthrough, the displacement pattern may either continue to evolve or become arrested, depending on specific system parameters. We propose simple scaling relationships that describe the change in displacement rates at breakthrough and predict the conditions under which growth of the displacement pattern ceases at breakthrough.
Micromodel geometry used in the viscous fingering experiments and simulations. The porous medium consists of cylindrical posts arranged between two parallel flat plates. (a) Typical image of viscous fingering simulation, where the fluid–fluid interface is highlighted in blue, while the pore-space colouring represents local fluid pressure. (b) The radii of posts and lattice size are chosen such that (c) the average throat diameter is comparable to the gap thickness
$h$
between the micromodel plates. (d) Photo of the experimental set-up, showing a silicone-oil-filled micromodel illuminated from below. The invading water phase is dyed blue for visualisation. (e) The static advancing contact angle of a water droplet on an acrylic substrate submerged in silicone oil was measured to be
$146^\circ$
.
Figure 1. Long description
The image contains five elements: one diagram, two bar graphs, one photo, and one close-up photo. Panel A shows a diagram of a viscous fingering simulation with the fluid-fluid interface highlighted in blue. The pore-space coloring represents local fluid pressure, with a color scale indicating pressure in MPa. Panel B is a bar graph showing the number of posts versus post diameters in micrometers. Panel C is another bar graph showing the number of throats versus throat diameters in micrometers. Panel D is a photo of the experimental setup, showing a silicone-oil-filled micromodel illuminated from below, with the invading water phase dyed blue for visualization. Panel E is a close-up photo showing the static advancing contact angle of a water droplet on an acrylic substrate submerged in silicone oil, measured to be 146 degrees. The diagrams and graphs illustrate the intricate dendritic patterns that emerge when a less viscous fluid displaces a more viscous one within a porous medium, highlighting the instability due to viscosity contrast.
2. Experimental set-up and observations
We conducted our experiments in a micromodel that consists of an array of cylindrical posts confined between two parallel plates (figure 1
a). Posts were arranged on a hexagonal lattice with a centre-to-centre post spacing of
$2.8$
mm. Post radii follow a uniform random distribution (figure 1
b), resulting in throats (open spaces between posts) with a mean diameter
$a=594\;\unicode{x03BC}$
m and a standard deviation of
$68\;\unicode{x03BC}$
m (figure 1
c). The micromodel has an outer radius of
$R=7.5$
cm and an inner post-free region radius of
$R_0=8$
mm (see figure 1
a). The micromodel was fabricated by laser etching the space around an array of posts on a transparent 6 mm-thick acrylic sheet and sealing it with a second flat acrylic sheet (figure 1
d). The radii of the channel and outlet tube
$r_o$
were selected so that the associated viscous pressure drops were negligibly small compared with the pressure drop through the micromodel (
$\Delta {p}_{\textit{outlet}}/\Delta {p}_{\textit{micromodel}}=\mathcal{O}(10^{-2})$
).
Initially, the micromodel was saturated with polydimethylsiloxane silicone oil (Clearco,
$\mu _2=96.5$
mPa s). Distilled water (
$\mu _1=0.89$
mPa s) was then injected through the central port using a syringe pump (New Era Pump Systems) until the injected volume matched the micromodel’s total pore-space volume. The static advancing contact angle of water phase on an acrylic substrate surrounded by silicone oil was measured to be
$146^\circ$
(figure 1
e). Throughout the experiment, fluid injection was performed at a constant flow rate while maintaining an approximately constant pressure at the outer boundary. We characterise the injection conditions using the dimensionless capillary number
$ \textit{Ca} = {\mu _1 U}/{\sigma }$
, where
$\sigma =13$
mN m–1 denotes interfacial tension, and
$U = {q}/{(h a)}$
represents the characteristic displacement velocity (Zhao et al. Reference Zhao, MacMinn and Juanes2016). Here,
$h = 600\;\unicode{x03BC}$
m is the spacing between the micromodel plates, and
$q$
is the volumetric injection flow rate.
Experimental images of post-breakthrough evolution of a single pore volume of water injected into a micromodel filled with 96.5 mPa s silicone oil (
$M=9.2\times 10^{-3}$
) at injection rates corresponding to
$3.2\times 10^{-5}\leqslant {Ca}\leqslant 9.0\times 10^{-2}$
. Pore spaces invaded by water at breakthrough are marked in red, while the regions invaded after breakthrough are marked with white. The red dashed line corresponds to that shown in figure 4. See supplementary movie is available at https://doi.org/10.1017/jfm.2026.11962 for the temporal evolution of the patterns.
We begin by fixing the viscosity ratio,
$M \equiv \mu _1 / \mu _2$
, and conducting a series of experiments, where we displace the silicone oil with water at progressively decreasing capillary numbers. These experiments reveal qualitatively distinct post-breakthrough evolution of the displacement patterns (figure 2): at high capillary numbers (
$ \textit{Ca}$
), the invading fluid continues to advance and its saturation increases, even after breakthrough; in contrast, at low
$ \textit{Ca}$
, the displacement pattern effectively ‘freezes’ following breakthrough. These contrasting outcomes prompt several questions: Under what conditions does the displacement pattern growth cease after breakthrough? How does the pattern growth rate evolve after breakthrough, and how does this depend on system parameters such as the viscosity ratio
$M$
and the domain size?
3. Pore-network model
We can make considerable progress towards answering these questions through a moving-capacitor pore-network model, which has been described in detail elsewhere (Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2019, Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021). We include a brief summary of the model here for completeness, focusing on the injection of the non-wetting fluid that is of primary interest in this study.
The moving-capacitor model describes fluid–fluid interface evolution in porous media by taking into account both viscous and Laplace pressure drops across the network. The model approximates the flow domain as a network of capillaries and assigns capillary-entry pressures to each pore throat (
$\Delta p_{\textit{crit}}$
). A Hagen–Poiseuille flow profile is used to calculate the pressure drop across each capillary. Thus, the flow of two immiscible fluids through the network is governed by a coupled system for the network edge fluxes
$\boldsymbol{q}(t)$
and node pressures
$\boldsymbol{p}(t)$
, expressed as
\begin{align} \begin{bmatrix} \boldsymbol{C}^{-1}(t) & \boldsymbol{A} \\[5pt] \boldsymbol{A}^T & \boldsymbol{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{q}(t) \\[5pt] \boldsymbol{p}(t) \end{bmatrix} = \begin{bmatrix} \boldsymbol{b}(t) \\[5pt] \boldsymbol{f} \end{bmatrix}\!, \end{align}
where
$\boldsymbol{C}(t)$
is the diagonal matrix of edge conductivities (
$\pi a^4/128\mu L$
,
$L$
is the distance between nodes),
$\boldsymbol{A}$
is the incidence matrix (which encodes the connectivity between nodes),
$\boldsymbol{b}(t)$
stores the capillary pressure drops across network edges and
$\boldsymbol{f}$
accounts for the fluxes at the inlet pores. The conductivities and
$\boldsymbol{b}(t)$
are updated at each time step using the filling ratio of pore throats
$\boldsymbol{\varPhi }(t)$
, and adaptive Euler time stepping is used to evolve the interface.
This model has been verified against experimental data and other numerical models elsewhere (Zhao et al. Reference Zhao2019). We use it here to reproduce Lenormand’s (
$M$
,
$ \textit{Ca}$
) phase diagram for our micromodel geometry and system wettability. The experimentally measured static advancing contact angle of water on acrylic submerged in silicone oil is approximately
$146^\circ$
(figure 1
e). However, because the invading interface is mobile and our pore-network model does not account for dynamic wetting, we adopt a fixed effective contact angle of
$170^\circ$
, which lies between the measured static value and the limiting value of
$180^\circ$
expected for a rapidly advancing interface. The resulting phase diagram, corresponding to strong drainage conditions (Lenormand et al. Reference Lenormand, Touboul and Zarcone1988; Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021), is plotted in figures 3(a) and 4(a). It contains the three displacement modes characteristic of non-wetting fluid invasion: compact displacement, invasion percolation and viscous fingering.
(a) Outcomes of the moving-capacitor model synthesised in (
$M,{Ca}$
) phase diagram, where viscous fingering, compact displacement and invasion–percolation region boundaries are outlined with a dashed line and background colour corresponds to the fractal dimension of the displacement pattern at breakthrough. Here, red and white points correspond to effective (
$M,{Ca}$
) values before and after breakthrough, respectively. Panels p1–p3 show corresponding simulations of displacement patterns before and after the breakthrough in red and white, respectively. Temporal evolution of the (b) injection pressure and (c) invading-fluid saturation in a typical viscous fingering simulation; breakthrough time is indicated with a vertical dashed line. (d) Numerical approximation of
$\tilde {{Ca}}/{Ca}$
obtained from rate of change in saturation after and before breakthrough, as demonstrated in (c).
Figure 3. Long description
The image contains multiple graphs depicting the outcomes of a moving-capacitor model in a phase diagram. Panel A shows a heatmap with regions labeled as viscous fingering, compact displacement, and invasion-percolation. The x-axis is labeled Log(M) and the y-axis is labeled Log(Ca). The background color corresponds to the fractal dimension of the displacement pattern at breakthrough. Red and white points correspond to effective values before and after breakthrough, respectively. Panels p1, p2, and p3 show corresponding simulations of displacement patterns before and after breakthrough in red and white, respectively. Panel B is a line graph showing the temporal evolution of injection pressure in a typical viscous fingering simulation. The x-axis is labeled t (s) and the y-axis is labeled Pressure (M P a). A vertical dashed line indicates the breakthrough time. Panel C is a line graph showing the temporal evolution of invading-fluid saturation in a typical viscous fingering simulation. The x-axis is labeled t (s) and the y-axis is labeled Saturation. A vertical dashed line indicates the breakthrough time. Panel D is a scatter plot showing the numerical approximation of obtained from the rate of change in saturation after and before breakthrough. The x-axis is labeled M and the y-axis is labeled Ca/Ca. The color of the points corresponds to the value of Ca.
4. Results and discussion
4.1. Change in finger growth rate at breakthrough
When
$M\ll 1$
and the flow is viscously dominated, the pressure drop within the invading fluid is negligibly small before breakthrough. At breakthrough, the invading-fluid pressure remains uniform in space except for the finger that just connected to the outer boundary (see figure 1
a). Therefore, the highest invading-fluid pressure gradients are within that finger, which means that the breakthrough finger carries the majority of the flow. If we take
$a$
to be the characteristic pore-throat size, and assume that the finger width and
$h$
are comparable to
$a$
, we can approximate the pressure drop along the breakthrough finger as
\begin{align} \Delta p_{\textit{finger}} \sim \frac {\mu _1 U R}{a^2}. \end{align}
In viscously dominated flow and for
$M\ll 1$
, the pressure within the invading fluid is nearly uniform in space. This is true even at breakthrough for all but the breakthrough finger (see figure 1
a).
In a homogeneous porous medium, viscous fingers reach the boundary within a short interval of one another. Consequently, when the first finger breaks through, the fronts of the remaining fingers lie approximately a distance
$a$
from the boundary. As the degree of disorder increases, the breakthrough times of different fingers become increasingly disparate (see figure 16(a) in Primkulov et al. Reference Primkulov, Pahlavan, Fu, Zhao, MacMinn and Juanes2021). Nevertheless, the characteristic distance between the remaining finger fronts and the boundary remains of order
$a$
. We therefore assume that this distance scales as
$a$
, so that the pressure drop ahead of the remaining fingers scales as
\begin{align} \Delta p_{\textit{def}} \sim \frac {\mu _2 \tilde {U}}{a}, \end{align}
where
$\tilde {U}$
is the rate of finger growth after breakthrough. Since the pressure drop within the invading fluid is negligible for non-breakthrough fingers, the two pressure drops must be equal (
$\Delta p_{\textit{finger}}\approx \Delta p_{\textit{def}}$
). Balancing these two pressure drops yields
$\tilde {U} \sim MRq/a^3$
or
\begin{align} \tilde {\textit{Ca}}/{\textit{Ca}} \sim M \, \frac {R}{a}, \end{align}
where
$ \textit{Ca}$
and
$\tilde {{Ca}} = {\mu _1 \tilde {U}}/{\sigma }$
represent the dimensionless rates of finger growth before and after breakthrough.
Equation (4.3) suggests that the rate of viscous finger growth decreases sharply at breakthrough for
$M\ll 1$
, and this decrease is more apparent for small domains due to the ratio of macroscopic to pore scale
$R/a$
. At breakthrough, the change in the effective capillary number produces a shift on Lenormand’s phase diagram (figure 3
a), and one can understand the change in the displacement mode at breakthrough as a result of this shift in effective
$\textit{Ca}$
. For example, while points p1–p3 are equidistant from the boundary of the viscous fingering regime in figure 3(a) before breakthrough (red circles), they correspond to different
$M$
and thus different magnitudes of capillary number shift after breakthrough (white circles). While point p1 remains in the viscous fingering regime after the shift, point p3 crosses over to the invasion–percolation regime, and this is evident from the change of finger sizes after breakthrough in figure 3(a) (white vs red in panel p3).
We can examine the validity of (4.3) more quantitatively via comparison with our pore-network model. We plot the evolution of the injection pressure and saturation of the invading fluid for a typical viscous fingering within our micromodel in figure 3(b,c). At breakthrough, the pressure at the injection port drops sharply (figure 3
b). This drop coincides with the change in the invading-fluid saturation rate,
$\dot {S}(t)$
. We calculate the ratio of
$\dot {S}(t)$
immediately after to just before breakthrough (represented by the two red tangents in figure 3
c). This ratio quantifies the viscous finger growth rate across breakthrough and thus corresponds to
$\tilde {{Ca}}/{Ca}$
. We plot
$\tilde {{Ca}}/{Ca}$
estimated from the viscous fingering simulations in figure 3(d) along with the scaling relation (4.3). As expected, the linear scaling between
$\tilde {{Ca}}/{Ca}$
and
$M$
holds reasonably well for small
$M$
and breaks down as
$(M,{Ca})$
values approach the boundaries of the viscous fingering regime of Lenormand’s diagram (figure 3). We emphasise that the scaling in (4.3) is intended for the early post-breakthrough regime, where only one viscous finger has connected the inlet to the outer boundary. Once multiple fingers have broken through, additional geometric effects are expected to modify the scaling.
The critical capillary number
$ \textit{Ca}_\times$
from (4.4) is overlayed as a horizontal red dashed line over (a) the Lenormand phase diagram and (b) the values of the ratio of injection pressure after breakthrough to capillary-entry pressure values from viscous fingering simulations. The initial capillary numbers of simulation panels p4–p9 are marked as red dots in (a), and the white dots correspond to experiments from figure 2.
4.2. Cross-over to arrested growth
According to (4.3) the rate of finger growth after breakthrough approaches zero as
$M \rightarrow 0$
. However, this is not the only condition that would freeze the displacement pattern after breakthrough, as is evident from experiments in figure 2. Figure 3(b) shows that the pressure at the injection port decreases sharply at breakthrough, and this pressure is maintained across all fingers that have not yet reached the outer boundary (see figure 1
a). Therefore, one can expect the displacement pattern growth to arrest if this injection pressure drops below the mean capillary-entry pressure,
$\Delta p_{\textit{cap}}\sim \sigma /a$
, needed to squeeze the fluid–fluid interface past a pore throat with characteristic size
$a$
. Balancing
$\Delta p_{\textit{cap}}$
with
$\Delta p_{\textit{finger}}$
yields the condition for arrested finger growth after breakthrough
\begin{align} {Ca}_\times \sim a/R. \end{align}
When the initial capillary number is below
$ \textit{Ca}_\times$
(red dashed line in figure 4), the displacement pattern freezes after breakthrough. This is evident from panels p4–p9 in figure 4(a) and from the experiments in figure 2, highlighted by white circles in figure 4(a). In fact, the post-breakthrough injection pressure
$\tilde {P}_{\textit{inj}}$
becomes comparable to the capillary-entry pressure as
$ \textit{Ca}\rightarrow {Ca}_\times$
from above (figure 4
b). Below
$ \textit{Ca}_\times$
, the injection port pressure hovers around the characteristic capillary-entry pressure, such that
$\tilde {P}_{\textit{inj}}/P_{\textit{cap}}=\mathcal{O}(1)$
and the fluid displacement mode changes from viscous fingering to invasion–percolation. This is the reason that scaling (4.3) breaks down for lower values of
$ \textit{Ca}$
in figure 3(d).
Importantly,
$ \textit{Ca}_\times$
is independent of
$M$
and is only a function of the micromodel geometry. The red dashed line in figure 4 would move up or down in response to changes in characteristic pore throat size
$a$
and the size of the domain
$R$
: increasing the pore throat size or decreasing the size of the micromodel would make it more likely for the viscous fingering pattern to freeze after breakthrough (higher
$ \textit{Ca}_\times$
).
5. Conclusions
Our study presents simple scaling relations that describe how viscous fingering evolves after breakthrough. Equation (4.3) captures the sharp reduction in the growth rate of the displacement pattern at breakthrough, linking post-breakthrough behaviour to the viscosity ratio, domain size and pore scale. We further show that pattern growth ceases entirely once the injection pressure drops below the capillary-entry pressure, defining a cross-over capillary number,
$ \textit{Ca}_\times$
, determined solely by geometry (see (4.4)). Together, these results may prove useful in optimising sweep efficiency in applications such as enhanced oil recovery (Van Meurs Reference Van Meurs1957; Chuoke et al. Reference Chuoke, van Meurs and van der Poel1959), subsurface remediation, cleanup of reactor beds (Hill Reference Hill1952) and in situ mineral extraction (Seredkin et al. Reference Seredkin, Zabolotsky and Jeffress2016).
All of our simulations, as well as the outcomes of (4.3)–(4.4) apply to incompressible fluid pairs. Compressibility brings another layer of complexity after the breakthrough, where interfaces and fluid fingers may disconnect and rearrange intermittently after the breakthrough.
Furthermore, we have limited our study to the injection of a non-wetting fluid (drainage), and we expect the picture to change if one injected a strongly wetting fluid instead (strong imbibition). In drainage, fluid–fluid interfaces pinned between cylindrical posts act to resist the injection; in strong imbibition, fluid–fluid interfaces pull the invading fluid into the porous medium, and the displacement is rate limited by the injection port flow rate. Therefore, the manner in which the pattern evolves after the breakthrough may very well be very different from what we studied here. In fact, Odier et al. (Reference Odier, Levaché, Santanach-Carreras and Bartolo2017) showed that when strongly wetting fluids are injected into porous media, they advance through wetting films, which coarsen after breakthrough.
Finally, while we did our simulations in the radial micromodel geometry, we expect that (4.3)–(4.4) would hold true in the linear displacement. Verifying (4.3)–(4.4) in a linear displacement geometry through experiments or simulations, as well as exploring how the invading-fluid saturation evolves long after the breakthrough of the first finger, would be a compelling direction for future work.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2026.11962.
Acknowledgements
The authors gratefully acknowledge helpful discussions with C. MacMinn, A. Pahlavan and L. Wilen.
Funding
This work was supported by the US Department of Energy (grant no. DE-SC0018357).
Declaration of interests
The authors declare that they do not have any competing interests.
h
146∘
M=9.2×10−3
3.2×10−5⩽Ca⩽9.0×10−2
M,Ca
M,Ca
Ca~/Ca
Ca×
