# Comparison between ”Parametrized homology via zigzag persistence” and
”Refined homology in the presence of a real-valued continuous function”

Dan Burghelea  Department of Mathematics, The Ohio State University, Columbus,
OH 43210,USA. Email<EMAIL_ADDRESS>
###### Abstract

One shows that the persistence diagrams $Dgm^{\cdots}(H\mathbb{X})$ and the
measures $\mu^{\cdots}_{H\mathbb{X}}$ defined in [5] can be derived as
restrictions of the maps $\delta^{f},$ $\gamma^{f}$ and of the measures
$\dim\mathbb{F}(\cdots),$ $\dim\mathbb{T}(\cdots)$ considered in [2] and [1]
to the points and the rectangles above and below diagonal.

###### Contents

  1. 1 Introduction
  2. 2 Definitions
  3. 3 Equalities

##  1 Introduction

In a recent paper [5] published in this journal the authors have proposed for
an $\mathbb{R}-$ space $\mathbb{X}$ (= a continuous map $f:X\to\mathbb{R}$),
under reasonable topological hypotheses, four persistence diagrams

$Dgm^{\wedge}(H\mathbb{X}),Dgm^{\lor}(H\mathbb{X}),Dgm^{\backslash\backslash}(H_{r}\mathbb{X}),Dgm^{//}(H\mathbb{X})$

regarded as the relevant invariants for parametrized homology of $\mathbb{X}.$
These persistence diagrams collect the four types of barcodes which can be
associated to $f$ via zigzag persistence, are maps with discrete support from
$\overline{R}^{2}_{+}=\\{-\infty\leq a<b\leq+\infty\\}$ to $\mathbb{Z}_{\geq
0}$ and can be interpreted as densities $d\mu$ of four integer valued measures
$\mu=\mu^{\wedge}_{H_{r}(\mathbb{X})},$
$\mu^{\lor}_{H_{r}(\mathbb{X})},\mu^{\backslash\backslash}_{H_{r}(\mathbb{X})},\mu^{//}_{H_{r}(\mathbb{X})}$
for squares $R=[a,b]\times[c.d],$ $-\infty\leq a<b<c,d\leq\infty.$

In cite [2] sections 6 and 7 (cf. also [3]published in this journal) and with
more details in the book [1] sections 5 and 6, under essentially the same
hypotheses, two maps with discrete support,
$\rho^{f}_{r}:\mathbb{R}^{2}\to\mathbb{Z}_{\geq
0},\gamma^{f}_{r}:\mathbb{R}^{2}\setminus\Delta\to\mathbb{Z}_{\geq 0}$ with
$\Delta=\\{(x,x)\in\mathbb{R}^{2}\\}$, have been defined and studied. 111 In
the above references the discussion refers mostly to $X$ compact but the
conclusions remain the same without the need of any additional hypotheses when
$X$ is locally compact and $f:X\to\mathbb{R}$ is a proper map. In this case
the support of the maps $\delta_{r}^{f}$ and $\gamma^{f}_{r}$ is discrete;
this is already used in the case of angle-valued maps. In the case $X$ is
compact the support of the maps $\delta^{f}_{r}$ and $\gamma^{f}_{r}$ is
finite, hence these maps are configurations of points with multiplicity, cf.
[2], [1] and [3].

These maps are derived from the vector space-valued maps $\hat{\rho}^{f}_{r}$
and $\hat{\gamma}^{f}_{r}$ with $\rho^{f}=\dim\hat{\gamma}^{f}_{r}$ and
$\gamma^{f}_{r}=\dim\hat{\gamma}^{f}_{r}$ and are viewed as relevant
topological invariants for a real-valued map $f.$ The support of $\rho_{r}$
consists of points in $\mathbb{R}^{2}$ which correspond to the closed $r-$bar
codes and open $(r-1)-$ barcodes while the support of $\gamma_{r}$ of points
which correspond to the closed-open and open-closed $r-$barcodes. In both [2]
and [1] these two maps 222They become four when one treats separately the
sigma algebras generated by the above diagonal squares and the below diagonal
squares are defined in two ways, one being as ”densities” of the measures
$\dim\mathbb{F}_{r}$ and $\dim\mathbb{T}_{r}$ on the sigma algebras associated
to the same type of squares. A measure theoretic formulation of both
$\delta^{f}_{r},\gamma^{f}_{r},$ and even more general of
$\hat{\delta}^{f}_{r}$ and $\hat{\gamma}^{f}_{r},$ can be explicitly found in
[1] subsection 9.2 and is implicit in [2] and [1] subsections 5.1 and 6.1.

The purpose of this note is to show that the persistence diagrams
$Dmg^{\cdots}_{\cdots}$ and of the measures $\mu^{\cdots}_{\cdots}$ are the
restrictions to the points and the rectangles above the diagonal resp. below
diagonal after composing with the map $T(x,y)=(y,x)$ of the maps
$\rho^{f}_{\cdots}$ and the measures
$\dim\mathbb{F}_{\cdots},\dim\mathbb{T}_{\cdots}$ facts which might pass
unnoticed in view of notational differences.

Precisely, one has:

###### Proposition 1.1

  * •

(a): $Dmg^{\wedge}(H_{r}(\mathbb{X}))$ equals $\rho^{f}_{r}$ restricted to
$\mathbb{R}^{2}_{+},$

  * •

(b): $Dmg^{\lor}(H_{r}(\mathbb{X}))$ equals $\rho^{f}_{r-1}\cdot T$ restricted
to $\mathbb{R}^{2}_{+},$

  * •

(c): $Dmg^{\backslash\backslash}(H_{r}(\mathbb{X}))$ equals $\rho^{f}_{r}$
restricted to $\mathbb{R}^{2}_{+},$

  * •

(d): $Dmg^{//}(H_{r}(\mathbb{X}))$ equals $\rho^{f}_{r-1}\cdot T$ restricted
to $\mathbb{R}^{2}_{+}.$

$(R^{2}_{+}=\\{a,b)\in\mathbb{R}^{2},a<b\\}$) and

###### Proposition 1.2

  * •

(A):
$\mu^{\wedge}_{H_{r}}([a,b]\times[c,d))=\dim\mathbb{F}_{r}((a,b]\times[c,d)),$

  * •

(B):
$\mu^{\backslash\backslash}_{H_{r}}([a,b]\times[c,d])=\dim\mathbb{T}_{r}((a,b]\times(c,d]),$

  * •

(C):
$\mu^{//}_{H_{r-1}}([a,b]\times[c,d])=\dim\mathbb{F}_{r}((c,d]\times[a,b)),$

  * •

(D):
$\mu^{\lor}_{H_{r}}([a,b]\times[c,d])=\dim\mathbb{F}_{r}((c,d]\times[a,b)).$

Note that the stability results as stated in [5] follows in a straightforward
manner from the stability results of [2] or [1] and the Alexander duality in
[5] can be derived without effort from the Poincaré-duality in [2] or [1] in
the same way the Alexander-duality can be derived from the Poincaré duality.
(NOTE: After the posting of the first version of this note it was brought to
my attention that the result on Alexander-duality as well as most of the
arguments in [5] were contained in the thesis of one of the author, Sara
Kalisnik, cf. http://www.matknjiz.si/doktorati/2013/Kalisnik-14521-4.pdf, and
posted on arXiv cf. Sara Kalisnik , Alexander Duality for Parametrized
Homology, arXiv:1303.1591.)

##  2 Definitions

###  2.1 CSKM-definitions

For simplicity in writing one shortens the notations in [5] by replacing the
notation : $\wedge,$ $\lor,$ $\backslash\backslash,$ $//$ by $c,$ $o,$ $co,$
$oc$ abreviatioins of closed, open, closed-open, open-closed and
$\mu^{\cdots}_{H_{r}\mathbb{X}}$ by $\mu^{\cdots}_{r}.$

Consider

$\mathcal{B}^{c}_{r}=$ the multi-set of closed $r-$bar codes,

$\mathcal{B}^{o}_{r}=$ the multi-set of open $r-$ bar codes,

$\mathcal{B}^{c,o}_{r}=$ the multi-set of closed-open $r-$bar codes,

$\mathcal{B}^{o,c}_{r}=$ the multi-set of open-closed $r-$ bar codes.

The definitions of barcodes $\mathcal{B}^{\cdots}_{r}$ in [5] (called in [5]
”intervals” and / or ”decorated pairs”) are based on the initial presentation
of zigzag persistence introduced by Carlsson, de-Silva, Morozov in 2009 . A
reformulation of these definitions in terms of ”death” and of ”observability”
is provided in [1] subsection 9.1.1.

If for a barcode $I$ one denotes by $l(I)$ resp. $r(I)$ the left end resp. the
right end, then a careful reading of the definitions in [5] shows that for a
box $R=[a,b]\times[c,d]$ with $-\infty\leq a<b<c<d\leq\infty$ one has:

$\displaystyle\mu_{H_{r}}^{c}(R)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{c}\ \ \mid a<l(I)\leq b,\ c\leq
r(I)<d\\}$ (1) $\displaystyle\mu_{H_{r}}^{c,o}(R)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{c,o}\mid a<l(I)\leq b,\ c\leq
r(I)<d\\}$ $\displaystyle\mu_{H_{r}}^{o}(R)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{o}\ \ \mid a<l(I)<b,\ c<r(I)<d\\}$
$\displaystyle\mu_{H_{r}}^{o,c}(R)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{o,c}\mid a<l(I)\leq b,\ c\leq
r(I)<d\\}$

and then

$\displaystyle Dmg^{c}_{r}(a,b)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{c}\mid a=l(I)\ r(I)=b\\}$ (2)
$\displaystyle Dmg^{c,o}_{r}(a,b)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{c,o}\mid a=l(I)\ r(I)=b\\}$
$\displaystyle Dmg^{o}_{r}(a,b)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{o}\mid a=l(I)\ r(I)=b\\}$
$\displaystyle Dmg^{o,c}_{r}(a,b)=$
$\displaystyle\sharp\\{I\in\mathcal{B}_{r}^{o,c}\mid a=l(I)\ r(I)=b\\}$

###  2.2 BH-definitions

Denote by :

$\mathbb{I}_{a}(r)=\mathrm{i}mg(H_{r}(f^{-1}((-\infty,a]))\to H_{r}(X)),$

$\mathbb{I}^{a}(r)=\mathrm{i}mg(H_{r}(f^{-1}([a,\infty)))\to H_{r}(X)),$

$\mathbb{F}_{r}(a,b)=\mathbb{I}_{a}(r)\cap\mathbb{I}^{b}(r),$
$F_{r}(a,b):=\dim\mathbb{F}_{r}(a,b),$

$\mathbb{T}_{r}(a,b):=\ker(H_{r}(f^{-1}((-\infty,a])\to
H_{r}(f^{-1}((-\infty,b]))$ when $a<b,$

$\mathbb{T}_{r}(a,b):=\ker(H_{r}(f^{-1}([a,\infty))\to f^{-1}([b,\infty))$
when $a>b,$

$T_{r}(a,b)=\dim\mathbb{T}_{r}(a,b).$

For a box $B=(a,b]\times[c,d),\ a<b,c<d$ one defines

$\mathbb{F}_{r}(B):=\mathbb{F}_{r}(b,c)/(\mathbb{F}_{r}(a,c)+\mathbb{F}_{r}(b,d));$

and one observes

$\dim\mathbb{F}_{r}(B)=F_{r}(b,c)+F_{r}(a,d)-F_{r}(a,c)-F_{r}(b,d).$ (3)

For a box above diagonal $B^{\prime}=(a,b]\times(c,d],a<b\leq c<d$ one defines

$\mathbb{T}_{r}(B^{\prime}):=\mathbb{T}_{r}({b,d})/j^{\prime}\mathbb{T}_{r}{(a,d)}+\mathbb{T}_{r}{(b,c)},$

with $j^{\prime}:\mathbb{T}_{r}{(a,d)}\to\mathbb{T}_{r}{(b,d)}$ the obviously
induced linear map, and one observes

$\dim\mathbb{T}_{r}(B^{\prime})=T_{r}(b,d)+T_{r}(a,c)-T_{r}(a,d)-T_{r}(b,c).$
(4)

For a box below diagonal $B^{\prime\prime}=[c,d)\times[a,b),a<b\leq c<d$ one
defines

$\mathbb{T}_{r}(B^{\prime\prime}):=\mathbb{T}_{r}(c,a)/j^{\prime\prime}\mathbb{T}_{r}(d,a)+\mathbb{T}_{r}(c,b)$

with $j^{\prime\prime}:\mathbb{T}_{r}({d,a})\to\mathbb{T}_{r}(c,a)$ the
obviously induced linear map, and one observes

$\dim\mathbb{T}_{r}(B^{\prime\prime})=T_{r}(c,a)+T_{r}(d,b)-T_{r}(c,b)-T_{r}(d,a).$
(5)

Recall from [4] Theorem 3.2 and Proposition 5.3 or from [2] Proposition 4.1
333a real valued map can be regarded as an angle valued map or from [1]
Proposition 4.3 the following equalities:

  * •

$\dim
H_{r}(f^{-1}(-\infty,a])=\begin{cases}\sharp\\{I\in\mathcal{B}_{r}^{c}\mid
l(I)\leq a\\}\\\ \sharp\\{I\in\mathcal{B}_{r-1}^{o}\mid
I\subset(-\infty,a)\\}\\\ \sharp\\{I\in\mathcal{B}_{r}^{co}\mid l(I)\leq
a<r(I)\\}\end{cases}$

and then

$\dim\mathbb{I}_{a}(r)=\begin{cases}\sharp\\{I\in\mathcal{B}_{r}^{c}\mid
l(I)\leq a\\}\\\ \sharp\\{I\in\mathcal{B}_{r-1}^{o}\mid
I\subset(-\infty,a)\\}\end{cases},$

  * •

$\dim
H_{r}(f^{-1}([a,\infty))=\begin{cases}\sharp\\{I\in\mathcal{B}_{r}^{c}\mid
r(I)\geq a\\}\\\ \sharp\\{I\in\mathcal{B}_{r-1}^{o}\mid
I\subset(a,\infty)\\}\\\ \sharp\\{I\in\mathcal{B}_{r}^{oc}\mid l(I)<a\leq
r(I)\\}\ \end{cases}$

and then

$\dim\mathbb{I}^{a}(r)=\begin{cases}\sharp\\{I\in\mathcal{B}_{r}^{c}\mid
r(I)\geq a\\}\\\ \sharp\\{I\in\mathcal{B}_{r-1}^{o}\mid I\subset(a,\infty)\\}\
\end{cases}.$

As a consequence we have

  1. 1.

for $a>b$

$\dim\mathbb{F}_{r}(a,b)=\begin{cases}\sharp\\{I\in\mathcal{B}_{r}^{c}\mid
l(I)\leq a,r(I)\geq b\\}\\\ \sharp\\{I\in\mathcal{B}_{r-1}^{o}\mid
b<l(I)<r(I)<a\\}\ \end{cases}$ (6)
$\delta^{f}_{r}(a,b)=\sharp\\{I=[a,b]\in\mathcal{B}_{r-1}^{o}\mid
l(I)=b,r(I)=a\\}$ (7)

  2. 2.

for $a\leq b$

$\dim\mathbb{F}_{r}(a,b)=\sharp\\{I\in\mathcal{B}_{r}^{c}\mid l(I)\leq a\leq
b\leq r(I)\\}$ (8)
$\delta^{f}_{r}(a,b)=\sharp\\{I=[a,b]\in\mathcal{B}_{r}^{c}\mid
l(I)=a,r(I)=b\\}$ (9)

  3. 3.

for $a<b$

$\dim\mathbb{T}_{a,b}(r)=\sharp\\{I\in\mathcal{B}_{r}^{co}\mid l(I)\leq
a<r(I)\leq b\\}$ (10)
$\gamma^{f}_{r}(a,b)=\sharp\\{I=[a,b)\in\mathcal{B}_{r-1}^{o}\mid
l(I)=a,r(I)=b\\}$ (11)

  4. 4.

for $a>b$

$\dim\mathbb{T}^{a,b}(r)=\sharp\\{I\in\mathcal{B}_{r}^{oc}\mid b\leq
l(I)<a\leq r(I)\\}$ (12)
$\gamma^{f}_{r}(a,b)=\sharp\\{I=[a,b]\in\mathcal{B}_{r-1}^{o}\mid
l(I)=b,r(I)=a\\}$ (13)

##  3 Equalities

For $R=[a,b]\times[c,d]$ with $a<b<c<d$

(1)+ (8) +(3) imply Proposition1.2 (A) and (2)+ (9) imply Proposition1.1(a),

(1)+ (6) +(3) imply Proposition 1.2(B) and (2)+ (7) imply Proposition1.1(b),

(1)+ (10) +(4) imply Proposition1.2(C) and (2)+ (11) imply Proposition1.1(c),

(1)+ (12) +(4) imply Proposition1.2(D) and (2)+ (13) imply Proposition1.1(d).

## References

  * [1] Dan Burghelea, New topological invariants for real- and angle-valued maps; an alternative to Morse-Novikov theory Word scientific Publishing Co. Pte. Ltd, 2017
  * [2] Dan Burghelea, Stefan Haller, Topology of angle valued maps, bar codes and Jordan blocks, J. Appl. and Comput. Topology, vol 1, pp 121-197, 2017, arXiv:1303.4328; Max Plank preprints
  * [3] D.Burghelea, A refinement of Betti numbers and homology in the presence of a continuous function I, Algebr. Geom. Topol. 17. pp 2051-2080, 2017, arXiv 1501.02486
  * [4] D. Burghelea and T. K. Dey, Persistence for circle valued maps. Discrete Comput. Geom.. 50 2013 pp69-98 ; arXiv:1104.5646
  * [5] Gunnar Carlsson, Vin de Silva, Sara Kalisnik, Dmitry Morozov. Parametrized homology via zigzag persistence, Algebr. Geom. Topol. 19 657-700, 2019


# Disorder-induced non-linear growth of viscously-unstable immiscible two-
phase flow fingers in porous media

Santanu Sinha<EMAIL_ADDRESS>PoreLab, Department of Physics,
University of Oslo, N–0316 Oslo, Norway  Yves Méheust  yves.meheust@univ-
rennes1.fr CNRS, Géosciences Rennes, UMR 6118, Université de Rennes 1,
F–350042 Rennes, France  Hursanay Fyhn<EMAIL_ADDRESS>SINTEF Energy
Research, Trondheim N-7465, Norway  PoreLab, Department of Physics, Norwegian
University of Science and Technology, N–7491 Trondheim, Norway  Subhadeep Roy
<EMAIL_ADDRESS>Department of Physics, Birla Institute
of Technology and Science Pilani, Hyderabad Campus, Telangana 50078, India
Alex Hansen<EMAIL_ADDRESS>PoreLab, Department of Physics, Norwegian
University of Science and Technology, N–7491 Trondheim, Norway

###### Abstract

The immiscible displacement of a fluid by another one inside a porous medium
produces different types of patterns depending on the capillary number Ca and
viscosity ratio $M$. At high Ca, viscous fingers resulting from the viscous
instability between fluid-fluid interfaces are believed to exhibit the same
Laplacian growth behavior as viscously-unstable fingers observed in Hele-Shaw
cells by Saffman and Taylor [1], or as diffusion limited aggregates (DLA) [2].
I.e., the interface velocity depends linearly on the local gradient of the
physical field that drives the growth process (for two-phase flow, the
pressure field). However, steady-state two-phase flow in porous media is known
to exhibit a regime for which the flow rate depends as a non-linear power law
on the global pressure drop, due to the disorder in the capillary barriers at
pore throats. A similar nonlinear growth regime was also evidenced
experimentally for viscously-unstable drainage in two-dimensional porous media
20 years ago [3]. Here we revisit this flow regime using dynamic pore-network
modeling, and explore the non-linearity in the growth properties. We
characterize the previously-unstudied dependencies of the statistical finger
width and nonlinear growth law’s exponent on Ca, and discuss quantitatively,
based on theoretical arguments, how disorder in the capillary barriers
controls the growth process’ non-linearity, and why the flow regime crosses
over to Laplacian growth at sufficiently high Ca. In addition, the statistical
properties of the fingering patterns are compared to those of Saffman-Taylor
fingers, DLA growth patterns, and the results from the aforementioned previous
experimental study.

###### pacs:

62.20.mm, 64.60.av, 46.50.+a, 81.40.Np

##  I Introduction

Fingering patterns are one of many unique features of two-phase flow which are
caused by hydrodynamic instabilities between two fluids [4, 5]. When one fluid
displaces another fluid inside a medium, depending on the properties of the
two fluids and the medium in which they are flowing, the displacement front
may exhibit fingering instead of a stable interface between the two fluids,
hereafter denoted the front [6, 7, 8]. A wide variety of fingering patterns
are observed for different types of multiphase flow, such as the flow of
miscible [9, 10, 11, 12] and immiscible [3, 13] fluids in continuum [1, 4] and
porous medium [14, 15, 16], reactive transport flow [17] and the flow of
frictional fluids in granular materials [18, 19]. The structural properties of
different types of fingers are controlled by the underlying physical forces,
such as the viscous, capillary, inertial and frictional forces. Furthermore,
depending on the driving condition or the geometry of the system, the
displacement front can also undergo transition from fingering patterns to
compact [20] or foam-like structures [21].

Our study here deals with non-reactive two-phase flow in non-deformable
medium, where the two fluids are immiscible and separated by interfaces
associated with a surface tension. Inside a continuum medium, for example a
Hele-Shaw cell, which consists of two parallel plates separated by a small
gap, the displaced fluid having a higher viscosity than the displacing fluid
leads to viscous instability of the front and the development of fingers with
a constant width, except in the vicinity of the finger tip. The smooth rounded
shape of these fingers in a two-dimensional $zx$ plane can be described by the
following parametric equation given by Saffman and Taylor in their pioneering
work in 1958 [1],

$\displaystyle
z(x)=\frac{W(1-\lambda)}{2\pi}\ln\left[\frac{1}{2}\left(1+\cos\frac{2\pi
x}{\lambda W}\right)\right]$ (1)

where $z$ is the direction of the overall front propagation, $W$ is the width
of the flow cell, and the parameter $\lambda$ is the ratio of the width of the
finger to the width of the channel. Note that this theory does not prescribe
$\lambda$ and that the surface tension between the two fluids does not enter
this equation, as it was not taken into account in its derivation. However,
Saffman and Taylor experimentally measured values of $\lambda$ close to $0.5$
over a wide range of flow rates, and it was later found that the finger width
was selected by surface tension [22]. When the flow takes place inside a
porous medium, on the other hand, for example a Hele-Shaw cell randomly filled
with glass beads, the fingers exhibit a fractal structure [23]. Properties of
these fractal fingers are controlled by the competition between viscous and
capillary forces, the ratio of which is called the capillary number, and the
viscosity contrast of the two fluids [24, 25, 26]. For slow injection, the
process is controlled by the disorder in the capillary forces at the pores,
and the displacement front generates capillary fingers with a fractal
dimension similar to invasion percolation clusters [27]. For fast injection of
a lower-viscous fluid into a higher-viscous fluid inside a porous medium [28],
the growth of the front is governed by viscous instabilities, and it produces
fractal viscous fingers with a lower fractal dimension [29, 16]. It was
pointed out that the statistical properties of viscous fingers in porous media
are analogous to diffusion limited aggregation (DLA) [30], which is a process
of aggregation of matter limited by the diffusion of random walkers arriving
from a far distance with a steady flux [31]. Both of Saffman-Taylor viscous
fingering and DLA follow Laplacian growth,

$\displaystyle\nabla^{2}p=0$ (2)

where $p$ is the probability density of the random walker for DLA or the
pressure drop across the interface between the two fluids for viscous
fingering [32, 2, 3]. For Saffman-Taylor viscous fingering, the flow in each
of the fluid phases is governed by Darcy’s law [33, 34],

$\displaystyle\bm{v}=\frac{\kappa}{\mu}\bm{\nabla}p\;,$ (3)

where $\bm{v}$ is the fluid velocity, while $\kappa$ and $\mu$ are
respectively the permeability of the porous medium and the fluid’s viscosity.
Taking into account the flow incompressibility, $\bm{\nabla}\cdot\bm{v}=0$ for
incompressible fluids, the Laplace equation $\nabla P^{2}$ ($P$ being the
pressure field) holds everywhere in both fluids. It follows from this that the
same equation 2 holds for the pressure drop across the interface, and thus the
interface displacement is driven by Laplacian growth if the fluids’
viscosities are different.

The linearity between the pressure drop and the velocity or the flow rate
indicated by Eq. 3 does not, however, hold for steady-state flow in a certain
range of capillary numbers. Steady state flow implies the simultaneous flow of
both fluids in the porous medium for a sufficiently long time so that the
system has reached a state when the statistical averages of macroscopic
quantities, such as the saturation or the global pressure drop, do not drift
with time anymore while both fluids are still flowing. Experiments in a porous
medium consisting of a mono-layer of glass beads in which air and water-
glycerol mixture were displaced simultaneously have shown that the total flow
rate $Q$ in steady state varied with the applied pressure drop $\Delta P$ as a
non-linear power law of exponent $\approx 1.85$ [35, 36]. Later experiments
with two incompressible fluids in the same porous medium found the exponent to
be $\approx 1.35$ or $1.5$ depending on the fractional flow [37]. Experiments
with three-dimensional porous media made of glass bead packings [38, 39, 40]
and real core samples [41, 42, 43], performed by different groups, have
further established this non-linearity with different exponent values. Various
theoretical approaches [35, 36, 44, 45], and numerical modeling with variable-
radii tubes [46], capillary bundles [47], pore-network models [44, 40] and
lattice Boltzmann simulations [48] have been carried out to understand the
origin of the non-linearity that makes the rheology to deviate from Eq. 3. In
general, the non-linear power law can be expressed by the expression

$\displaystyle Q\propto(\Delta P-P_{\rm t})^{\beta}~{},$ (4)

where $P_{\rm t}$ is a threshold pressure that may exist so that there will be
no flow in the system below $P_{\rm t}$, and $\beta>1$ is the non-linear
exponent. Fundamentally, this non-linearity is related to the distribution of
the capillary barriers at the pore throats, which result from surface tension
at the menisci between the two fluids; disorder in the pore geometry induces
disorder in the spatial distribution of the capillary barriers. With the
increase in $\Delta P$, more and more pores progressively allow for the
barrier to be overcome and flow to be conducted, which makes $Q$ increase
faster than the increase in $\Delta P$, and hence $\beta>1$ [49]. When all the
available pores along the interface, at any time, start flowing at a
sufficiently high $\Delta P$, the relation becomes linear. For a simplified
porous system such as a capillary bundle, the value of $\beta$ can be
determined analytically by integrating the individual flow rates of the tubes
over the distribution of the capillary barriers [47]. For a porous network,
the distribution of pore throat sizes [50] and the wettabilities [51, 52]
strongly control the value of $\beta$, as the distribution of capillary
barriers depends on them.

In this article, we investigate to which extent this non-linearity in the flow
rate, which has been demonstrated to exist for the steady state due to
disorder in the capillary barriers, also exists for the growth of viscously-
unstable fingers during drainage, and what are its characteristics. The Eq. 2
for DLA considers $p=0$ at the perimeter of the aggregate, which, fingering in
two-phase flow, would correspond to the absence of surface tension. The
analogy to DLA therefore only applies for capillary numbers that are
sufficiently high for capillary forces to be insignificant as compared to
viscous forces. In the intermediate regime between capillary and viscous
fingerings, we may expect that capillary forces will compete with viscous
forces, so that the disorder in capillary barriers will play a role in the
finger growth. Løvoll et al [3] and Toussaint et al [53] have investigated
this displacement regime experimentally using two-dimensional (2D) porous
media consisting of Hele-Shaw cells filled with a monolayer of glass beads,
and measured the fractal dimension of the fingers and the statistical width of
the front [30, 54, 55], comparable to $\lambda$ in Eq. 1. Based on these
measurements, they showed that the drainage fingers are characteristically
different from DLA. They further proposed a quadratic relationship between the
finger growth rate and the local pressure gradient in the intermediate flow
regime, resulting from pore-scale disorder. For DLA however, an ad-hoc type
surface tension was considered by introducing a sticking probability [32],
which showed no change in the fractal dimension [56, 57]. It was then argued
that the two-phase flow fingers are more similar to dielectric breakdown
models (DBM) (or the $\eta$ models) [58, 59, 60], which are generalizations of
DLAs with power-law relationships of exponent $\eta$ between the growth
probability and the local growth/displacement-driving field.

The experimental study on 2D porous media [3, 53] left a number of open
questions, in particular concerning the dependence of the nonlinear growth
exponent and statistical width of the viscous fingers on the capillary number.
In this paper we study viscously-unstable drainage as a function of the
capillary number using large-scale simulations in dynamic pore networks. We
inject a lower-viscosity fluid at one edge of the network filled with a higher
viscosity fluid, and characterize the statistical properties of the resulting
viscous finger in a reference frame attached to the most advanced finger tip,
as the finger grows inside the porous medium. We explore in particular the
relationship between the finger’s growth rate and the local pressure drop
across the interface, and how it depends on the capillary number. We also
focus on the statistical width parameter $\lambda$, analog to the $\lambda$ in
Eq. 1, compare it with the values obtained in the previous experimental study
as well as in studies of DLA and Saffman-Taylor viscous fingering, and
characterize its dependence on the capillary number.

In the following, we first describe the computational model and the simulation
procedure in Sec. II, which is based on dynamic pore-network modeling [61]
with a specific boundary condition. In Sec. III we present our results, where
we first characterize the statistical profiles related to the volume and
growth of the fingers in subsection III.1. In subsection III.2 we characterize
the shape profile and measure the width ratio $\lambda$. In subsection III.3
we then explore the relationship between the growth rate and the local
pressure drop across the interface. Finally, we discuss the results and
provide an overall conclusion in Sec. IV.

##  II System description and modeling

Figure 1:  A pore-network of dimension $L\times W$ consisting of $N_{W}=20$
and $N_{L}=80$ links. The length of each link is $l$ and therefore
$h=l/\sqrt{2}$. The network was initially filled with high-viscosity fluid
(represented in gray). The low-viscosity fluid (represented in blue) is
injected with a total flow rate $Q$ through $N_{I}=2$ inlet nodes at the
bottom, indicated by the black dots. All the other nodes at the bottom
boundary, and all on the two vertical boundaries are closed, which is
indicated by the thick black lines. The nodes at top boundary are open and
work as the outlets. The $(x,z)$ coordinate axes indicated by the orange lines
are attached to the most advanced fingertip and move with it. The part of the
network between the two red horizontal lines is used for the measurements of
different quantities in order to avoid finite size effects at the inlet and
the outlet.

The core of our simulation consists of a dynamic pore-network model [61],
which we tweaked to adapt to the present problem. The model has been developed
for over a decade [62] and tested against numerous experimental [63, 40],
theoretical [44, 64, 65, 66] and lattice-Boltzmann simulations [67]. In this
computational method the pore space of a porous medium is modeled by a network
of links and nodes. Our network is spread in 2-dimensions (2D) within a flat
cuboid domain similar to a Hele-Shaw cell, and consists of $N_{W}\times N_{L}$
number of links embedded in a diamond lattice. Such a network with $N_{W}=20$
and $N_{L}=80$ is shown in Fig. 1. Here the subscripts $W$ and $L$ refer to
the directions of the network orthogonal and parallel to the overall flow
respectively. Each link of the network is of length $l=1\,\text{mm}$. The
total network is therefore $W=hN_{W}\,\text{mm}$ wide and
$L=hN_{L}\,\text{mm}$ long, where $h=l/\sqrt{2}$. We assign the entire pore-
space of the network to the links, so the nodes only represent the position of
the link intersections. The correspondence between this geometry and that of a
granular quasi-two-dimensional porous medium consisting of cylinders is
illustrated in Fig. 2a, where a small part of the porous medium is shown. The
links are therefore composite in nature, which means that each link must
represent a narrow pore throat in between two wider pore bodies. Such a
geometry is modeled by an hourglass-shaped converging-diverging tube (as
illustrated in Fig. 2b) whose cross-sectional area varies along its length.
This results in a variation in the capillary pressure $p^{\prime}$ as the
interface (i.e., meniscus) between the two fluids moves along a link. We model
this variation with a modified Young-Laplace equation [68, 62, 45],

$\displaystyle|p^{\prime}_{i}(\epsilon)|=\frac{2\gamma}{r_{i}}\left[1-\cos\left(\frac{2\pi\epsilon}{l}\right)\right]\;,$
(5)

where $y$ is the position of an interface inside a link $i$. Such variation in
the capillary pressure is shown in Fig. 2b. Here
$\gamma=\gamma^{\prime}\cos\theta$, where $\gamma^{\prime}$ and $\theta$
represent the surface tension between the two fluids and the contact angle
respectively. The $r_{i}$ is the average radius of the $i$th link, the value
of which we choose from a uniform distribution of random numbers in the range
between $0.1l$ to $0.4l$, which introduces the disorder in the network.

The capillary numbers for both the viscous and capillary fingering regimes are
sufficiently small for any inertial effect to be negligible. For fully
developed laminar flow of incompressible fluids in each link, we therefore
consider the following equation for the flow rate of the fluids [69, 24],

$\displaystyle q_{i}=-\frac{k_{i}}{l\mu_{i}}\left[\Delta
p_{i}-\sum_{b}p^{\prime}_{i}(\epsilon_{b})\right]~{},$ (6)

where $\Delta p_{i}$ is the pressure drop between the two nodes across the
$i$th link. The network is assumed to be placed horizontally and no gravity is
considered. Here $\mu_{i}$ is the effective viscosity of the tho fluids inside
the link, which is given by $\mu_{i}=s_{i}\mu_{\rm n}+(1-s_{i})\mu_{\rm w}$,
where $\mu_{\rm n}$ and $\mu_{\rm w}$ are the viscosities of the non-wetting
and wetting fluids respectively and $s_{i}$ is the non-wetting saturation in
the $i$th link, i.e., the proportion of the link volume that is occupied by
the non-wetting fluid. The term $k_{i}$ represents the mobility of the link,
given by $k_{i}=a_{i}r_{i}^{2}/8$ for Hagen-Poiseuille flow in a circular
cross-section where $a_{i}=\pi r_{i}^{2}$ is the cross-sectional area of the
link. The summation over the capillary pressure $p^{\prime}_{i}$ in Eq. 6 is
over all the interfaces inside the link $i$, obtained using Eq. 5. For the
link shown in Fig. 2b, for example, the summation will be over the three
interfaces at $\epsilon_{1}$, $\epsilon_{2}$ and $\epsilon_{3}$.

To find out the pressures $p_{j}$ at every node $j$ of the network, we use the
Kirchhoff equation for incompressible flow $\sum_{i\in n_{j}}q_{i}=0$ for each
node. Here the summation is over the links $n_{j}$ connected to a node $j$.
This provides a set of linear equations which we solve by conjugate gradient
method [70] with proper boundary conditions. This is illustrated in Fig. 1
where the two nodes ($N_{I}=2$) at the bottom edge act as the inlets. All the
nodes at the top edge work as the outlets through which fluids leave the
system. The two vertical boundaries of the network parallel to the overall
flow are closed. The network is initially filled with wetting fluid (gray
colored in the figures), and we inject the non-wetting fluid (blue colored in
the figures) through inlet nodes with a total constant flow rate $Q$. We
therefore set $p_{j}=0$ for all the outlet nodes and $q_{i}=Q/N_{I}$ for all
the virtual links connected to the inlet nodes, indicated by the two arrows in
Fig. 1. These serve as the boundary conditions for solving the linear system
of equations for the node pressures.

(a)  (b)

Figure 2:  (a) Illustration of the network of pores and nodes. The pore space
is indicated by white color whereas the gray circles represent the grains in a
porous rock or the glass beads in a Hele-Shaw cell. One of the links is
colored by dark gray. The links intersect at nodes at the intersections of the
dashed lines. (b) Shape of individual pores in the mid-horizontal plane of the
network, and variation of the capillary pressure $p^{\prime}$, given by Eq. 5,
along the length of the pore. The white and gray segments inside the pore
represent respectively the wetting and non-wetting fluid blobs; there are
three interfaces between them in this case, which are in contact with the
beads at positions $\epsilon_{1}$, $\epsilon_{2}$ and $\epsilon_{3}$ along the
pore length.

The positions of the two fluids inside every link are assigned by the
positions of the interfaces between the blobs of the two immiscible fluids in
every link as shown in Fig. 2b. At every time step, the displacements of the
two fluids are performed by updating these positions by a distance,

$\displaystyle\Delta\epsilon_{i}=\Delta t\,q_{i}/a_{i}$ (7)

in the direction of the flow in the corresponding link. Here $\Delta t$ is a
time step chosen in such a way that the displacement of an interface inside
any link do not exceed more than $0.1l$.

One final detail about the modeling is how to distribute the fluids from links
to their neighboring links through the nodes. For this, first the links
carrying fluids towards and away from every node are identified, we call them
the incoming and outgoing links, respectively. Then for every node $j$, the
total volumes of the wetting and non-wetting fluids ($V^{\rm w}_{j}$ and
$V^{\rm n}_{j}$ respectively, and $V_{j}=V^{\rm w}_{j}+V^{\rm n}_{j}$)
arriving from the incoming links to the node are measured using Eqn. 7. These
volumes are then distributed to each outgoing links $i$ by placing new wetting
and non-wetting blobs of volumes $V^{\rm w}_{i}=q_{i}\Delta tV^{\rm
w}_{j}/V_{j}$ and $V^{\rm n}_{i}=q_{i}\Delta tV^{\rm n}_{j}/V_{j}$
respectively. This means that the ratio of $V^{\rm w}_{i}$ to $V^{\rm n}_{i}$
in any outgoing link is identical to the ratio between the incoming wetting
and non-wetting volumes at the corresponding node, and the ratio between the
total volumes ($V^{\rm w}_{i}+V^{\rm n}_{i}$) injected in different outgoing
links $i$ is the same as the ratio between the flow-rates $q_{i}$ in those
links. Further technical details about this algorithm can be found in [61].

##  III Results and discussions

(a) $\mathrm{Ca}=0.004$ (b) $\mathrm{Ca}=0.008$

(c) $\mathrm{Ca}=0.01$ (d) $\mathrm{Ca}=0.04$

(e) $\mathrm{Ca}=0.08$ (f) $\mathrm{Ca}=0.1$

(g) $\mathrm{Ca}=0.4$ (h) $\mathrm{Ca}=0.8$

Figure 3:  Drainage fingers generated for viscosity ratio $M=10^{-4}$, and for
different values of Ca. The blue and black colors respectively denote the
invading non-wetting and the defending wetting fluids. The inlets are at the
bottom edge of the networks and the overall flow direction is towards the top.
The gray shade indicates the non-wetting volume density $\tilde{\sigma}(x,z)$,
averaged over time and samples in the moving coordinate system.

Drainage simulations are performed at constant flow rates $Q$ corresponding to
chosen the capillary numbers Ca. This is a dimensionless global parameter of
the flow, which is defined as the average ratio of the viscous forces to the
capillary forces at the pore level. It is defined as [24],

$\displaystyle\text{Ca}=\frac{\mu Q/A}{\gamma}$ (8)

where $A$ is the cross-sectional area. Here we point out that for our pore-
network only the areas of the channels are properly defined, so $A$ is the
orthogonal cross-section of the channels alone, and therefore $Q/A$ is the
true mean longitudinal fluid velocity, rather than the filtration (or Darcy)
velocity as usually considered in experiments. Note also that experimental
studies, including those addressing the same topic [3, 53] may have used a
different definition of Ca, considering in particular that the viscous
pressure drop should involve the medium’s permeability; Eq. 8 assumes that the
permeability can be equated to the square of the typical channel cross-
section. Furthermore, the wetting angle could be taken into account in the
capillary pressure drop, thus appearing in the definition of Ca as well. When
comparing our simulations results to experimental results under identical flow
conditions, the similarity of the flow conditions must be assessed based on Ca
estimates obtained from the same definition, so a conversion of the Ca values
(either experimental, or numerical) may be necessary.

To find out the value of $Q$ that must be imposed in the simulation for a
given Ca, we first define it as [71],

$\displaystyle{\rm Ca}=\frac{\Delta p_{\rm visc}}{\Delta p_{\rm
cap}}=\frac{\Delta P_{W}h/L}{2\gamma/{\bar{r}}}$ (9)

where $\Delta p_{\rm visc}=\Delta P_{W}h/L$ and $\Delta p_{\rm
cap}=2\gamma/{\bar{r}}$ are respectively the average viscous and capillary
pressure drops across a pore. Here $\Delta P_{W}$ is the total viscous
pressure drop across the network when it is completely filled with the wetting
fluid and ${\bar{r}}$ is the average link radius. We then first solve the
viscous pressure drop $\Delta P_{{\rm visc},1}$ for a test flow rate $Q_{1}$
for the network saturated with wetting fluid, and then determine the flow rate
$Q$ by using,

$\displaystyle Q=\frac{2\gamma L}{h{\bar{r}}}\frac{Q_{1}}{\Delta P_{{\rm
visc},1}}{\rm Ca}\;.$ (10)

This is because $Q\propto\Delta P_{W}$ for single phase flow. The simulation
is then performed with $q_{i}=Q/N_{I}$ for all the virtual inlet links.

Table 1: Number of samples simulated for different values of Ca and $M$. The values of $M$ are $10^{-4}$ and $10^{-5}$ whereas the values of Ca are in the range $0.01$ to $0.9$. All the measured quantities are statistically averaged over these many samples for the respective values of Ca and $M$. Ca: | 0.01 | 0.02 | 0.03 | 0.04 | 0.05 | 0.06 | 0.07 | 0.08 | 0.09 | 0.10 | 0.20 | 0.30 | 0.40 | 0.50 | 0.60 | 0.70 | 0.80 | 0.90
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$M=10^{-4}$ | 57 | 79 | 93 | 105 | 112 | 117 | 123 | 125 | 125 | 125 | 121 | 105 | 98 | 90 | 84 | 81 | 76 | 73
$M=10^{-5}$ | 30 | 49 | 56 | 66 | 71 | 74 | 81 | 81 | 82 | 81 | 80 | 65 | 62 | 57 | 54 | 51 | 46 | 45

We have considered a network of dimension $32\times 160$ links in this study,
this is based on the maximum computational time we could spend reasonably. We
had to make additional sweeps through the entire network at every time step to
compute the volume density, the growth density and the local pressure drops,
which we will discuss later. This made the simulations considerably more time-
consuming compared to usual pore-network simulations with the same model.
Furthermore, in the analyses we discarded the regions of the network that are
distant from the inlet and the outlet by less than $N_{W}$ rows, as shown by
the red lines in Fig. 1, in order to avoid any boundary effect. In Tab. 1 we
show the numbers of different realizations of the network considered for
different values of $M$ and Ca. All the measured quantities in different plots
in this article are averaged over these many samples. The surface tension
between the fluids was chosen to $\gamma=0.03\,\text{N/m}$. Two viscosity
ratios, $M=\mu_{\rm n}/\mu_{\rm w}=10^{-4}$ with $\mu_{\rm
n}=10^{-5}\,\text{Pa.s}$ and $\mu_{\rm w}=10^{-1}\,\text{Pa.s}$, and
$M=10^{-5}$ with $\mu_{\rm n}=10^{-5}\,\text{Pa.s}$ and $\mu_{\rm
w}=1\,\text{Pa.s}$, were considered. A typical set of simulations for
different capillary numbers and viscosity ratios are shown in Fig. 3, where
the black and blue colors represent the wetting and non-wetting fluids
respectively. The direction of the overall flow in these images is from the
bottom to the top as indicated by the arrow. The two red lines indicate the
aforementioned regions near the inlets and the outlets that are not taken into
account in the analyses. The gray shades indicate the average volume density,
which we describe in the following.

Figure 4:  Plot of the dimensionless growth density $h\phi(z)$ as a function
of the row numbers from the tip, $z/h$, for the viscosity ratios (a)
$M=10^{-4}$ and (b) $M=10^{-5}$. The different sets in each plot correspond to
different values of Ca indicated in the legend of both plots. The solid line
corresponds the model function $B\exp(-\frac{z}{h\xi})$ with $\xi=4.6$ and
$B=0.18$, for both plots. The inset shows a linear-logarithmic plot of the
same data where the dashed line has a slope of $-1/\xi$ and intercept of $\ln
B$.

###  III.1 Longitudinal profiles of volume density and volumetric growth rate

The volume density of the fingers represents a statistical map of the
occupation of the two-dimensional network space by the invading non-wetting
fluid. This provides a statistical shape of the displacement structure which
is comparable with the smooth continuum-medium fingers [1]. The density maps
are functions of the distance from the most advanced fingertip, $z$, as the
displacement process can be considered statistically stationary in an $(x,z)$
reference frame attached to that tip (see figure 1), as previously shown for
experimental drainage fingers [3, 53], as well as for DLA fingers [30, 54].
With respect to this dynamic $(x,z)$ reference frame, we define the volume
density $\tilde{\sigma}(x,z)$ as,

$\displaystyle\tilde{\sigma}(x,z)=\frac{1}{h^{2}n(x,z)}\sum_{k=1}^{n(x,z)}V_{k}(x,z)\;,$
(11)

which is essentially the volume of injected fluid per unit area at a given
position in the $(x,z)$ reference plane, averaged over time, $V_{k}(x,z)$
being the non-wetting volume in the link positioned at position ($x,z$) at
time step $k$. Here we point out that the lowest unit of resolution for our
measurements is one single link, which means that we have only one data point
within an area of $h^{2}$. Furthermore, as the available window of network for
the measurements varies with time, the number of elements withing the
summation is different for different links. This is accounted by $n(x,z)$ in
the above expression, which is the number of times the volume in the link at
$(x,z)$ is added in the summation. This also means that the data near the
fingertip contain better statistical averages compared to those far behind the
tip. The map is then further averaged over a number of different samples of
the network as indicated in Tab. 1 to obtain statistically averaged
$\tilde{\sigma}(x,z)$ map. In Fig. 3, we show the spatial distribution of
$\tilde{\sigma}(x,z)$ as gray-scale maps where the dark and light shades
correspond to the smallest and largest values, respectively. For an observer
sitting at the origin of the moving coordinate system, this volume profile of
the finger will remain the same statistically during the invasion of the
fingers in an infinitely long network, and only the length of the profile will
increase in the $z$ direction.

We next measure the volumetric growth rate $\tilde{\nu}(x,z)$ of the invading
fingers, which is defined as the average increase in the volume of the finger
per unit time and unit area at $(x,z)$

$\displaystyle\tilde{\nu}(x,z)=\frac{1}{h^{2}n(x,z)}\sum_{k=1}^{n(x,z)}\frac{1}{\Delta
t_{k}}\left[V_{k}(x,z)-V_{k-1}(x,z)\right]\;,$ (12)

where $\Delta t_{k}$ is the time interval between the time steps $k$ and
$k-1$. Note that, $\tilde{\nu}(x,z)$ has the dimension of a velocity. By
integrating $\tilde{\nu}(x,z)$ in both $x$ and $z$ directions we have

$\displaystyle\int_{0}^{L}\int_{0}^{W}\nu(x,z)\operatorname{d\\!}x\operatorname{d\\!}z=Q\;,$
(13)

as the total increase in the volume of the whole finger per unit time is equal
to the volumetric rate of the steady fluid injection into the system. By
normalizing $\tilde{\nu}(x,z)$ by $Q$, we then define the growth density
$\tilde{\phi}(x,z)$,

$\displaystyle\tilde{\phi}(x,z)=\frac{1}{Q}\tilde{\nu}(x,z)\;.$ (14)

We assume that the system is statistically symmetric with respect to the
longitudinal mid-section ($x=w/2$) of the flow cell, and we define the
quantities $\sigma(z)$, $\nu(z)$ and $\phi(z)$ in the longitudinal direction
as a function of the sole distance $z$ from the tip, by integrating
$\tilde{\sigma}(x,z)$, $\tilde{\nu}(x,z)$ and $\tilde{\phi}(x,z)$ in the $x$
direction,

$\displaystyle\sigma(z)$
$\displaystyle=\int_{0}^{W}\tilde{\sigma}(x,z)\operatorname{d\\!}x\;,$ (15)
$\displaystyle\nu(z)$
$\displaystyle=\int_{0}^{W}\tilde{\nu}(x,z)\operatorname{d\\!}x\;,$
$\displaystyle\phi(z)$
$\displaystyle=\int_{0}^{W}\tilde{\phi}(x,z)\operatorname{d\\!}x\;,$

from which we can also define $\phi(z)=\nu(z)/Q$. In Fig. 4, we plot
$h\phi(z)$ as a function of the number of rows from the most advanced tip,
$z/h$, for different values of Ca and $M$. Near the fingertips, the data shows
an exponential-type decay of $\phi(z)$. This is indicated by the solid line
where we plot a model function $B\exp(-\frac{z}{h\xi})$ with $\xi=4.6$ and
$B=0.18$. Here $\xi$ is a characteristic decay length, or screening length,
which characterizes the active invasion zone. The value of $\xi$ appears to be
the same for the whole range of capillary numbers and viscosity ratios,
showing that it is a characteristic of the porous medium, which was also
observed in the drainage experiments Hele-Shaw cell filled with a mono-layer
of glass beads[3]. The same data are plotted in a linear-logarithmic scale in
the insets, which show that the exponential decay is valid for $z<W/2$ for
lower capillary numbers whereas the validity increases up to $z<W$ at high
capillary numbers. This is indicated by the linear dashed line which has a
slope $-1/\xi$ and intercept $\ln(B)$. This characteristic is also in
agreement with experiments [3] which showed similar exponential decay of
finger growth within the same range from the tip. Away from the tip, $\phi(z)$
is negligibly small as seen from the scattered data points in the insets,
which mean that the active growth of the fingers happen near the fingertip
whereas they are almost frozen far behind the tip.

Figure 5:  Plot of the position of the most advanced fingertip,
$\tilde{z}_{\text{tip}}$, as a function of the scaled time $t/t_{\mathrm{f}}$
where $t_{\mathrm{f}}$ is the time when the most advanced fingertip reaches
the final position in the simulation. The two plots correspond to the
viscosity ratios (a) $M=10^{-4}$ and (b) $M=10^{-5}$; the different symbols
correspond to the different values of Ca indicated in the plot legends.

Figure 6:  Plot of the scaled longitudinal volume density
$\sigma(z)/\sigma_{c}$ as a function of the row numbers $z/h$ behind the most
advanced fingertip. The two figures correspond the viscosity ratios (a)
$M=10^{-4}$ and (b) $M=10^{-5}$ and the different data sets correspond to
simulations with different injection rates $Q$, and thus to different
capillary numbers Ca. The values of Ca are indicated in the figure legends.
The average cumulative growth density $\Phi(z)$ is plotted as solid black
line, which shows the agreement with Eq. 19. The black, blue and red dashed
lines are plotted using Eq. 1 with $\lambda=0.4$, $0.5$ and $0.62$, which
correspond to the experimental observation for porous Hele-Shaw cell [3, 53],
the DLA [55, 72] and the Saffman-Taylor solution for low surface tension limit
[1] respectively. The insets show the variations of unscaled $\sigma(z)$ for
different values of Ca.

(a) $\mathrm{Ca}=0.04$ (b) $\mathrm{Ca}=0.08$ (c) $\mathrm{Ca}=0.1$

Figure 7: Volume density maps thesholded with
$\tilde{\sigma}(x,z)\geq\tilde{\sigma}_{\text{max}}/2$ for $M=10^{-4}$ and for
$\text{Ca}=0.04$, $0.08$, and $0.1$, as indicated below the figures. The gray
shades correspond to the same density map as shown in Fig. 3 for the
respective values of Ca, but now only the links preserved by the cut-off are
shown. The green, yellow and blue lines are drawn using the Saffman-Taylor
equation (Eq. 1) with $\lambda=0.4$, $0.5$ and $0.62$, which correspond to the
experimental observation for porous Hele-Shaw cell, the continuum-medium
Saffman-Taylor finger and the DLA finger respectively.

Note that $\sigma(z)$ and $\phi(z)$ are related to each other. Indeed, let us
now examine the displacement process in the referential $(x,\tilde{z})$
attached to the inlet boundary of the network, with $\tilde{z}$ and $z$
related to each other according to $z=\tilde{z}_{\rm tip}(t)-\tilde{z}=v_{\rm
tip}t-\tilde{z}$, where $v_{\rm tip}$ is the velocity of the most advanced tip
of the non-wetting fluid phase and $\tilde{z}_{\rm tip}(t)$ is its position
along $\tilde{z}$. That velocity is constant in time, as shown in Fig. 5,
where $\tilde{z}_{\rm tip}$ is plotted as a function of the time normalize by
the time $t$ scaled by $t_{f}$, where $t_{f}$ is the time at which the tip
reaches the outlet boundary of the porous network. At any time $t$, the volume
of non-wetting fluid in a slice of the network normal to the mean flow
direction and located between $\tilde{z}$ and $\tilde{z}+\Delta\tilde{z}$ is
$\tilde{\sigma}(\tilde{z},t)\Delta\tilde{z}$, with $\Delta\tilde{z}=\Delta z$.
It results from the cumulative growth of the fingers between the time
$t_{0}=\tilde{z}/v_{\rm tip}$ at which the most advanced tip first arrived at
$\tilde{z}$ and $t$, and can thus be expressed as

$\tilde{\sigma}(\tilde{z},t)\Delta z=\Delta z\int_{\tilde{z}/v_{\rm
tip}}^{t}\tilde{\nu}(\tilde{z},u)\operatorname{d\\!}u~{}.$ (16)

The steady volumetric density defined above in the $(x,z)$ reference frame,
$\sigma(z)$, is related to $\tilde{\sigma}(\tilde{z},t)$ through the following
relation:

$\forall t~{}~{}~{}\sigma(z)=\tilde{\sigma}\left(v_{\rm tip}t-z,t\right)~{}.$
(17)

A similar expression relates the volumetric growth rate expressed in the
moving reference frame, $\nu(z)$, and that expressed in the laboratory frame,
$\tilde{\nu}(\tilde{z},t)$. Hence,

$\displaystyle\forall t~{}~{}~{}\sigma(z)\Delta z$ $\displaystyle=\Delta
z\int_{t_{0}}^{t}\tilde{\nu}\left(v_{\rm
tip}t^{\prime}-z,t^{\prime}\right)\operatorname{d\\!}t^{\prime}~{}.$ (18)
$\displaystyle\mathrm{i.e.,~{}~{}}\sigma(z)\Delta z$ $\displaystyle=\Delta
z\int_{0}^{z}\nu(z^{\prime})\frac{\operatorname{d\\!}z^{\prime}}{v_{\rm tip}}$
$\displaystyle\mathrm{i.e.,~{}~{}}\sigma(z)\Delta z$
$\displaystyle=\frac{Q\Delta z}{v_{\rm
tip}}\int_{0}^{z}\phi(z^{\prime})\operatorname{d\\!}z^{\prime}~{}.$

Defining the cumulative growth density in the $z$ direction as
$\Phi(z)=\int_{0}^{z}\phi(z^{\prime})\operatorname{d\\!}z^{\prime}$ and
$\sigma_{\rm c}=Q/v_{\text{tip}}$ leads to

$\displaystyle\Phi(z)=\frac{\sigma(z)}{\sigma_{\rm c}}~{},$ (19)

which relates the volume density with the cumulative growth density.

Figure 8:  Plots of the average transverse volume density profile $\rho(x)$ as
a function of the transverse normalized coordinate $x/W$ are shown for (a)
$M=10^{-4}$ and (b) $M=10^{-5}$. The different symbols show the results for
different values of Ca and the solid line shows the profile averaged over all
datasets. The values of $\lambda$ measured using the two different methods,
(A) $\lambda=1/\rho_{\text{max}}$ and (B) $\lambda=(x_{+}-x_{-})/W$ where
$\rho(x_{+})=\rho(x_{-})=\lambda_{\text{max}}/2$, are shown in (c) as a
function of the capillary number.

To verify this functional form, the tip velocity is first measured from the
slopes as the straight lines in figure 5 multiplied by $t_{f}$. We then plot
$\sigma(z)$ scaled by $\sigma_{c}$, for different values of Ca and $M$ (Fig.
6). Notice that, though $\sigma(z)$ is a function of the capillary numbers as
shown in the insets, all $\sigma(z)/\sigma_{c}$ plots collapse on a single
master curve for different values of Ca. With Eq. 19, it therefore implies
that there is one unique cumulative growth density $\Phi(z)$, independent of
Ca. We therefore calculate the average cumulative growth density $\Phi(z)$
from the data in Fig. 4 and plot it in Fig. 6 as a solid black line. The line
shows an excellent agreement with the $\sigma(z)/\sigma_{c}$ plots, as
expected from Eq. 19. Such a growth property of the front is the
characteristics of the viscous dominated regime, or of a regime where viscous
forces compete with capillary forces, and therefore will not be observed in
pure capillary fingering regime where the invasion is completely controlled by
the disorder in the capillary thresholds.

###  III.2 The ‘continuum’ shape

The average volume densities of the fractal fingers can be used to compare
their shapes statistically with the smooth shapes of continuum-medium fingers,
as well as with the shapes of DLA. By mapping the volume density function
$\sigma(z)$ to the Saffman Taylor equation (Eq. 1) for continuum viscous
fingers, one can calculate the width ratio $\lambda$, and use that to compare
the statistical shape profile of different fingers. Before presenting the
detailed measurements of $\lambda$, we show a quick comparison in Fig. 7 where
only the part of the network with volume density
$\tilde{\sigma}(x,z)\geq\sigma_{\text{max}}/2$ are shown. This region
corresponds to the porous media with large invasion [73, 54, 55] and is
comparable to the Saffman-Taylor equation. For the viscous fingers in
continuum Hele-Shaw cell without any porous structure, the width ratio
$\lambda$ in Eq. 1 has a value of $0.5$ in the limit of negligible surface
tension, and it increases with the increase of surface tension [22, 74]. For
experiments with porous Hele-Shaw cells on the other hand, $\lambda$ was found
to be around $0.4$ [3, 53], whereas for off-lattice DLAs, $\lambda$ was found
to be around $0.62$ for linear channels [55, 72]. The solutions of Eq. 1 for
these three values are shown in Fig. 7, where the green, yellow and blue lines
correspond to $\lambda=0.4$, $0.5$ and $0.62$ respectively. The lines show
deviations from the gray-scale region and the deviations also depend on Ca. We
also show the solutions of Eq. 1 for these three values of $\lambda$ in Fig. 6
by the dashed lines, which again show some deviations from the results for the
drainage fingers here.

In order to measure the value of the width parameter $\lambda$, we measure a
transverse volume density $\rho(x)$ in a region where the growth is almost
frozen, that is, far behind the most advanced fingertip. If $L_{\rm f}$ is the
length of this frozen zone in the $z$ direction, then one has

$\displaystyle\int_{0}^{W}\int_{W}^{W+L_{\rm
f}}\tilde{\phi}(x)\operatorname{d\\!}z\operatorname{d\\!}x=\sigma_{c}L_{\rm
f}\;,$ (20)

as $L_{\rm f}$ increases with the same velocity $v_{\text{tip}}$. We therefore
define $\rho(x)$ as,

$\displaystyle\rho(x)=\frac{W}{\sigma_{c}L_{\rm f}}\int_{W}^{W+L_{\rm
f}}\tilde{\phi}(x,z)\operatorname{d\\!}z\;,$ (21)

with $\int_{0}^{W}\rho(x)\operatorname{d\\!}x=W$. In Fig. 8a and 8b, we plot
$\rho(x)$ as a function of the scaled transverse position $x/W$ for the two
different viscosity ratios. The datasets show a maximum at the middle and then
decrease on both sides. This is an expected behavior, similar to what is
observed for DLA [55] and for the experiments of viscous fingers in porous
media [3, 53]. However, the plots drift with Ca, indicating a variation of
$\lambda$ with Ca. From this data, the width ratio $\lambda$ can be estimated
in two ways: either from (A) $\lambda=(x_{+}-x_{-})/W$ where
$\rho(x_{+})=\rho(x_{-})=\lambda_{\text{max}}/2$ or alternatively from (B)
$\lambda=1/\rho_{\text{max}}$ [73, 54]. The two measurements are shown in Fig.
8c, where $\lambda$ seems to vary systematically with Ca. Interestingly, the
lowest value here is at $\approx 0.5$, which is the solution for the continuum
Saffman-Taylor fingers for negligible surface tension limit [1]. Higher values
of $\lambda$ were found for increasing surface tension for the continuum
fingers [74, 22, 75]. The value of $\lambda\approx 0.62$ for DLAs in linear
channels [55, 72] coincides with some of the data points here, however the
two-phase flow experiments in porous Hele-Shaw cells showed a much lower value
of $\lambda\approx 0.4$ [3, 53]. Note however that in contrast to the present
system, the capillary threshold distribution was not uniform in these
experiments, despite the fact that a uniform distribution was assumed for some
of the interpretations. Here for our system, $\lambda$ increases up to a much
higher value of around $0.8$. Such values are similar to the fingers in
divergent cells, which are linear cells with an angular wedge at the inlet, a
geometry that is in between the linear and circular geometries [76]. There,
$\lambda$ was found to increase when increasing the angle of the wedge, for
example $0.77$ to $0.82$ for a wedge angle of $90^{\circ}$ [55, 30]. Note
that, in our simulations, we inject fluid only through a few of the nodes at
the center of the network’s inlet boundary, and all other nodes on the two
sides of these injection nodes are blocked, see Fig. 1. This may have a wedge
effect on the value of $\lambda$, similar to the divergent inlets. However,
the fact that here $\lambda$ shows a minimum at an intermediate Ca and
increases on both sides, indicates that there is a combined effect of surface
tension and inlet geometry, and it therefore needs further in-depth study. We
leave this for future, as our main focus here is to explore the effect of
disorder on the relationship between the growth and pressure drop.

###  III.3 Growth rate versus local pressure drop

The volumetric growth of the region occupied by the invading fluid inside a
pore depends on whether the viscous pressure drop between the invading and
defending fluids is able to overcome the capillary barrier inside the pore and
thus to displace the interface between the two fluids. Assuming that the
movement of the interfaces follow Eq. 6, a theoretical methodology was
suggested in [3, 53] to find the relationship between the growth rate and the
local pressure drop across the fluid-fluid interface by integrating over the
distribution of the capillary barriers. We follow their approach in the
following and modify it to adapt our system. If $P^{\prime}$ is the threshold
pressure required to invade a link and $g(P^{\prime})$ is the normalized
distribution of the thresholds over the network, we may integrate over
$g(P^{\prime})$ to obtain the growth rate,

$\displaystyle\tilde{\nu}(x,z)$ $\displaystyle=\frac{1}{h^{2}}\int
q(x,z)g(P^{\prime})\operatorname{d\\!}P^{\prime}$ (22)
$\displaystyle=\frac{\kappa}{h^{2}}$ $\displaystyle\int(\Delta
p(x,z)-P^{\prime})\Theta(\Delta
p(x,z)-P^{\prime})g(P^{\prime})\operatorname{d\\!}P^{\prime}\;,$

where $\Delta p(x,z)$ is the pressure drop between the invading and defending
fluids inside a link $i$ at ($x,z$) and $q(x,z)$ is the corresponding flow
rate of the fluids in that link. Here we consider $\kappa=\bar{k}/\mu_{w}$,
where $\bar{k}$ is the average mobility of the pores and $\mu_{w}$ is the
viscosity of the high-viscosity defending fluid. A pore will be invaded only
when $\Delta p(x,z)>P^{\prime}$, hence the Heaviside function $\Theta(\Delta
p(x,z)-P^{\prime})$.

Figure 9:  Node pressures ($p_{j}$) obtained from conjugate gradient solver,
normalized by the maximum pressure $p_{\rm max}$ at that time step, shown by
colors for a network of $32\times 96$ links. The color scale shows the
normalized pressure values, where the darkest red corresponds to $p_{j}=0$ and
the lightest yellow correspond to the highest pressure. Here $\text{Ca}=0.05$
and $M=10^{-5}$. The injected fluid is colored blue. The box at the top right
shows a close view of a small slice of the network. Notice that the nodes
belonging to the low-viscosity finger are at much higher pressures than those
belonging to the defending fluid.

If we consider $P^{\prime}$ to be uniformly distributed between $P_{\rm t}$
and $P_{\rm m}$, in the form

$\displaystyle g(P^{\prime})=\frac{1}{G}\Theta(P^{\prime}-P_{\rm
t})\Theta(P_{\rm m}-P^{\prime})$ (23)

where $G=P_{\rm m}-P_{\rm t}$, then Eq. 22 becomes

$\displaystyle\tilde{\nu}(x,z)=\frac{\kappa}{h^{2}G}\int_{P_{\rm t}}^{P_{\rm
m}}(\Delta p-P^{\prime})\Theta(\Delta
p-P^{\prime})\operatorname{d\\!}P^{\prime}\;,$ (24)

where $\Delta p\equiv\Delta p(x,z)$. If the viscous pressure drop exceeds the
capillary threshold for all pores, i.e., if $\Delta p(x,z)>P_{\rm m}$, then
Eq. 24 reduces to

$\displaystyle\tilde{\nu}(x,z)$
$\displaystyle=\frac{\kappa}{h^{2}G}\int_{P_{\rm t}}^{P_{\rm m}}(\Delta
p-P^{\prime})\operatorname{d\\!}P^{\prime}$ (25)
$\displaystyle=\frac{\kappa}{2h^{2}(P_{\rm m}-P_{\rm t})}\left[-\left(\Delta
p-P^{\prime}\right)^{2}\right]_{P_{\rm t}}^{P_{\rm m}}$
$\displaystyle\tilde{\nu}(x,z)$
$\displaystyle=\frac{\kappa}{h^{2}}\left[\Delta p-\frac{P_{\rm m}+P_{\rm
t}}{2}\right]\;.$

The growth rate in this regime then varies linearly with the excess pressure
drop with respect to the mean capillary threshold $(P_{\rm m}+P_{\rm t})/2$.

For $\Delta p(x,z)<P_{\rm t}$ on the other hand, Eq. 24 reduces to

$\displaystyle\tilde{\nu}(x,z)$
$\displaystyle=\frac{\kappa}{h^{2}G}\int_{P_{\rm t}}^{\Delta p}(\Delta
p-P^{\prime})\operatorname{d\\!}P^{\prime}$ (26)
$\displaystyle\tilde{\nu}(x,z)$
$\displaystyle=\frac{\kappa}{2h^{2}}\frac{\left(\Delta p-P_{\rm
t}\right)^{2}}{P_{\rm m}-P_{\rm t}}\propto\left(\Delta p-P_{\rm
t}\right)^{2}\;.$

The growth rate thus varies quadratically with $[\Delta p(x,z)-P{\rm t}]$.
This regime corresponds to moderate capillary numbers where the pressure drop
$\Delta p(x,z)$ competes with the capillary thresholds, and the number of
invaded pores thus increases with an increase in $\Delta p(x,z)$.

If we further assume that $\Delta p(x,z)$ is a function of $z$ only, we may
average $\Delta p(x,z)$ over $x$, which yields

$\displaystyle\nu(z)\sim\left(\Delta P(z)-P_{\rm t}\right)^{2}\;,$ (27)

where $\Delta P(z)=\langle\Delta p(x,z)\rangle_{x}$. Hence, in this regime,
the average growth rate $\nu(z)$ depends quadratically on the excess pressure
drop at $z$, $[\Delta P(z)-P_{\rm t}]$. Such a quadratic relationship between
excess pressure drop and flow rate was also suggested for steady state two-
phase flow in porous media by mean-field calculations [44], and observed by
numerical simulations [44, 40, 48].

Figure 10:  Plots of the non-wetting pressure $P_{\rm n}(z)$, the non-wetting
pressures $P_{\rm w}(z)$ and the pressure drop $\Delta P(z)$, obtained using
Eq. 29 and 30, as a function of the distance from the most advanced fingertip
$z/h$, for different capillary numbers. The top row corresponds to $M=10^{-4}$
and the bottom row to $M=10^{-5}$.

Note however that the threshold distribution function $g(P_{\rm t})$
considered in the above derivation was assumed to be uniform. In our pore
network it is a more complex function, which is not straightforward to
determine. It not only depends on the pore properties such as the size
distribution [50] and the pore wettabilities [51, 52] but also on the
structure of the finger itself, that is, the distribution of the number of
fingertips at any $z$. We therefore may assume a generalization of Eq. 27,

$\displaystyle\nu(z)\sim\left(\Delta P(z)-P_{\rm t}\right)^{\beta}\;,$ (28)

where $\beta$ is an exponent of non-linearity. As mentioned in the
Introduction, such non-linear relationships with $\beta>1$ have been widely
observed for steady-state flow in a certain regime of capillary numbers [35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 47, 46, 44, 48].

In order to examine the form of Eq. 28, we need to measure two quantities: (a)
the average pressure difference between the defending and invading fluids,
$\Delta P(z)$, as a function of $z$, and (b) the threshold pressure $P_{\rm
t}$. At every time step, we obtain the pressures at the nodes of the network
from the conjugate gradient solver by solving the equations 6 and 5. An
example is shown in Fig. 9, where the pressure values are indicated by
different colors. Notice that the nodes connected by the invading non-wetting
low-viscosity fluid (blue links) are at much higher pressures than those that
are connected with the defending wetting fluid. The conjugate-gradient solver
does not distinguish between the wetting and non-wetting pressures; we
therefore identify, at every time step, the links filled with invading fluid
as those with a non-wetting saturation larger than a threshold value
($s_{i}>0.98$), and the nodes that are connected to them are marked as the
non-wetting nodes. The rest of the nodes are then marked as the wetting nodes.
We then measure the average excess wetting and non-wetting pressures as a
function of the distance $z$ from the most advanced fingertip as

$\displaystyle P_{\rm w}(z)$ $\displaystyle=\frac{1}{n_{\rm
w}(z)}\sum_{k=1}^{n_{\rm w}(z)}p_{k}(z)-P_{0}$ (29) $\displaystyle\text{and~{}
~{}}P_{\rm n}(z)$ $\displaystyle=\frac{1}{n_{\rm n}(z)}\sum_{k=1}^{n_{\rm
n}(z)}p_{k}(z)-P_{0}\;,$

where $n_{\rm w}(z)$ and $n_{\rm n}(z)$ are the number of wetting and non-
wetting nodes at $z$, and $P_{0}$ is the average wetting pressure at the most
advanced fingertip, $z=0$. The pressure drop $\Delta P(z)$ between the wetting
and non-wetting fluids at $z$ is then calculated as

$\displaystyle\Delta P(z)=P_{\rm n}(z)-P_{\rm w}(z)\;,$ (30)

which is then averaged over different time steps during the propagation of the
fingers, and over different samples.

Figure 11:  Dependence of the threshold pressures $P_{\rm t}$ on the capillary
numbers Ca for the two different viscosity ratios ($M$), (a) $10^{-4}$ and (b)
$10^{-5}$. The measured values of $P_{\rm t}$ represented by the red circles
are obtained by taking the averages of $\Delta P(z)$ for $z/h>76$, i.e. over
the last $20\%$ of the data points on the plateau of $\Delta P(z)$ in Fig. 10.
The estimated values of $P_{\rm t}$ represented by the blue triangles are
obtained by minimizing the least-square fit error of Eq. 28.

In Fig. 10, $P_{\rm n}(z)$, $P_{\rm w}(z)$ and $\Delta P(z)$ are plotted as
functions of the row distance $z/h$ from the fingertip for different values of
Ca and $M$. We can notice a few things here. First, the wetting pressure
$P_{\rm w}$ shows a sharp variation with $z$, whereas the non-wetting pressure
$P_{\rm n}$ do not show any noticeable variation. This is because of the low
viscosity of the invading fluid compared to the defending fluid. Secondly, the
pressure drop $\Delta P(z)$ starts from a high value at $z=0$, then falls
rapidly withing a distance $z<W/2$ and almost goes to a plateau. This behavior
is similar to the growth density $\phi(z)$ in Fig. 4, which showed a rapid
exponential fall within a similar distance, indicating that the growth
$\phi(z)$ is directly related to $\Delta P(z)$. Furthermore, note that the
plateau regimes of $\Delta P$ far behind the tip has finite values instead of
zero. This however is the stagnant zone, where the fingers do not show any
visible growth and $\phi(z)$ is essentially zero as seen in Fig. 4. Therefore
in this regime, the capillary barriers compete with the viscous pressure drop
and prevent the interfaces from moving. The plateau value thus corresponds to
the threshold pressure $P_{\rm t}$ in Eq. 28, below which there is no growth.
We can therefore measure $P_{\rm t}$ from this plateau regime, by taking the
average of $\Delta P(z)$ over its farthest data points. However, if we look
carefully on the plots, we notice that not all data sets at the plateau
exhibit a perfectly constant behavior, some of them are still decreasing
slowly and therefore would need a much larger system to saturate. We therefore
estimate the threshold pressure $P_{\rm t}$ in one more way where we minimize
the least-square fit error for Eq. 28. We try a set of trial values for
$P_{\rm t}$ around the plateau of $\Delta P(z)$ and then choose the one that
provides the lowest least-square error for the straight-line fit for
$\log[\nu(z)]$ vs $\log[\Delta P(z)-P_{\rm t}]$. These estimated values of
$P_{\rm t}$ are shown in Fig. 11 with blue triangles, where we also compare
them with the measured values of $P_{\rm t}$, which are calculated as the
average of $\Delta P$ over the farthest $20\%$ rows from the fingertip. Notice
that the values of $P_{\rm t}$ from these two different calculations are very
close to each other for the whole range of the capillary numbers. Furthermore,
$P_{\rm t}$ does not seem to depend on the viscosity ratio $M$, which confirms
that the threshold is a quantity controlled solely by the capillary barriers,
which are functions of the interfacial tension between the two fluids.

Figure 12: Plot of the local excess pressure drop
$\log_{10}[\Delta(P(z)-P_{\rm t}]$ between the invading and defending fluids
as a function of the local growth rate of the finger, $\log_{10}[\nu(z)]$, for
the viscosity ratios (a) $M=10^{-4}$ and (b) $M=10^{-5}$. Different sets
correspond to different values of Ca as indicated by the symbols. The
exponents $\beta$ defined in Eq. 28 is calculated from the slopes. The slopes
are calculated only for the range plotted with colored symbols, as the growth
far behind the fingertip is almost frozen and lead to noisy data points as
shown by the gray symbols. The blue and red straight lines are drawn with
slopes $1$ and $3$.

Figure 13: Dependence of the growth exponent $\beta$ obtained from Fig. 12 on
Ca. The error bars are of order $0.07$ towards the lower side of Ca whereas
they are $0.01$ on the higher side. They are too small compared to the scale
of the plot, and therefore not shown. The two types of symbols correspond to
the two values of $M$ as indicated in the legend of the plot.

Finally, we now set out to verify the agreement with Eq. 28 and determine the
exponent $\beta$. In Fig. 12, we plot the growth rate $\nu(z)$ as a function
of the excess local pressure drop $[\Delta P(z)-P_{\rm t}]$ in $\log$ scale,
where we use the values of $P_{\rm t}$ estimated from Fig. 11. The plots show
a linear trend over one to three decades, depending on the capillary number,
for the entire range of capillary numbers, with a linear range of the log-log
plot that is all the larger as Ca is larger. The power law exponents vary
systematically with Ca within a range of $1$ to $3$ as shown by the blue and
red straight lines in the plots. However, at very small values of the pressure
drops, the data is noisy and deviates from the linearity, which corresponds to
the gray symbols in the plots. This corresponds to the zone far behind the
fingertip, where the $\Delta P(z)$ is very close to $P_{\rm t}$ and the growth
structure is almost frozen. As mentioned before, the statistical averaging
becomes poorer as we move away from the fingertip, which therefore adds to the
statistical noise in that region. We have therefore disregarded these points
when fitting power laws to the data, and only rely on the points shown as
colored symbols in the plots. The exponent $\beta$ is thus measured from
least-square fitting a power law to the colored data points. In figure 13 we
plot the dependence of $\beta$ on Ca for the two viscosity ratios, which
exhibit very similar trend, i.e., two plateaus for $\mathrm{Ca}\lesssim
10^{-1.3}$ and $\mathrm{Ca}\gtrsim 10^{-0.2}$, and a crossover in between
them. On the lower capillary number plateau, the values fluctuate within a
range of $\approx 2.5$ to $3.0$. It then undergoes a crossover with the
increase of Ca, and then approaches another plateau towards $1$. This is very
similar to the rheological behavior of the widely studied steady-state flow,
where the total two-phase flow rate varies non-linearly with the excess
pressure drop at low capillary numbers, and then undergoes a crossover to a
linear Darcy regime at high Ca [35, 36, 37, 38, 39, 40, 41, 42, 43]. The
origin of that non-linearity and the linear crossover has been well explained
in terms of the pore-scale disorders in the capillary pressures [44, 47],
radii distribution [50] and wettabilities [51, 52]. Though the values of the
steady-state non-linear exponent has been reported to differ for different
systems, the range of them was similar to what we find here in Fig. 13.

##  IV Discussion and conclusions

We have addressed in this article how the growth rate of fingers during
immiscible two-phase flow (drainage) in porous media depends on the local
viscous pressure drops between the two fluids. In other words, we revisited
the comparison of immiscible viscous fingers with Saffman-Taylor viscous
fingering and DLA, for which a linear Laplacian growth is assumed. We
considered pore networks with a regular geometry but with disorder in the
pore/channel widths, and focused on the intermediate to higher capillary
numbers, which include the intermediate regime between capillary and viscous
fingers and the regime of pure viscous fingering. First, we looked into the
geometrical properties of the drainage fingers, and compared their statistical
shape in the reference frame of the advancing finger to the smooth Saffman-
Taylor fingers that develop in the continuum scale description of two-phase
flow in porous media. Saffman and Taylor mentioned that viscous fingering in
an empty Hele-Shaw cell without any porous structure is equivalent to what
should be expected in a porous medium [1]. Though this is true for single
phase flow, it is not necessarily true for two-phase flow, except perhaps if
the porous medium is a regular lattice (no disorder), because the continuum
description cannot account for the role of surface tension in a disordered
porous medium. Indeed, in an empty Hele-Shaw cell, the fluid-fluid interface
is continuous, and the long-range in-plane component of its curvature induces
long-range forces, whereas in a porous medium the interface is broken up in
many small menisci, and surface tension acts at the scale of these menisci,
which is the pore scale. The in-plane component of the curvature is then much
larger, and, furthermore, stochasticity in the pore sizes plays an important
role in shaping the displacement process by inducing a stochasticity in the
capillary pressure thresholds.

We thus measured the statistical longitudinal volumetric density and growth
rate profiles of the fractal drainage fingers in the dynamic reference frame
attached to the most advanced finger tip, and compared them with the
properties of smooth Saffman-Taylor fingers and of DLA growth structures. The
width ratio with respect to the medium’s width, $\lambda$, of the fingers’
statistical density map, was found to approach the $0.5$ Saffman-Taylor value
at intermediate capillary numbers Ca, and then increase to around $0.8$ when
either increasing or decreasing the capillary numbers from that intermediate
Ca range. The corresponding longitudinal volumetric density and growth rate
profiles behave accordingly. For Saffman-Taylor fingers in continuum-medium,
$\lambda$ was observed to increase with the increase of surface tension [74,
22, 75], and it also increases with the angle of inlet wedge for systems with
divergent wedges [55, 30]. In case of our present system, it may be due to a
combined effect of these two, and needs further study.

The maximum growth happens at the most advanced fingertip, and decreases
exponentially behind it. Far behind the tip, the growth is almost frozen: the
interfaces between the two fluids are held in place by the capillary barriers
and thus do not move. We show that this growth behavior is directly correlated
with the local pressure drop between the two phases, and controlled by the
distribution of the capillary forces. By computing numerically an effective
capillary threshold, we show that there exists a regime at intermediate
capillary numbers where the linear Laplacian growth property does not hold as
the local growth rate varies non-linearly with the excess local pressure drop
across the interface (i.e., the excess between the local porous pressure drop
and the capillary pressure threshold), $\Delta P$. This nonlinear regime of
the growth of the invading finger is explained by accounting for the disorder
distribution in the capillary thresholds in a theoretical assessment of the
link between $\Delta P$ and the local growth rate. Indeed, when $\Delta P$
falls within the range of the capillary pressure thresholds, increasing the
capillary number means that not only will the invasion velocity increase
linearly with $\Delta P$, but the number of pores along the front that are
invaded at any time will also increase, thus rendering the growth rate
nonlinear. For a uniform capillary threshold distribution, the resulting
growth law is expected to be quadratic. In our disorder pore networks, that
distribution is not uniform, and consequently the nonlinear growth exponent
$\beta$ depends on the capillary number. The numerical simulations provide a
range of $\beta$ values similar to that of the rheology exponent previously
reported for steady-state flow. In fact, it transitions from a plateau at
values as large as 2.5–3 at the smallest investigated capillary numbers, to
another plateau at 1 for the largest investigated capillary numbers. This
large capillary number limit of $\beta=1$ is to be expected, since it
correspond to configurations in which $\Delta P$ is larger than the upper
boundary of the medium’s capillary pressure thresholds; in such configurations
our theoretical assessment predicts a crossover of the flow regime to the
linear Laplacian growth.

This study opens many prospects for future studies. One of them is the
systematic investigation of the link between the capillary threshold
distribution of the porous medium and the nonlinear growth exponent $\beta$.
Another one would be to study whether the values of $\beta$ are comparable for
steady state flow and the present viscous fingering process in the same porous
geometry. In order for the distribution of the capillary barriers to only
depend on the pore network geometry, the fluid-fluid interface must sample the
entire barrier distribution at any time, otherwise the capillary
barrier/capillary threshold distribution will also depend on the flow
patterns, in a manner that may differ between the steady-state flow and
viscous fingering configurations.

## Acknowledgements

This work is supported by the Research Council of Norway through its Centers
of Excellence funding scheme, project number 262644.

## References

## References

  * Saffman and Taylor [1958] P. G. Saffman and G. I. Taylor, “The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid,” Proc. R. Soc. Lond. 245, 312 (1958).
  * Paterson [1984] L. Paterson, “Diffusion-limited aggregation and two-fluid displacements in porous media,” Phys. Rev. Lett. 52, 1621 (1984).
  * Løvoll  _et al._ [2004] G. Løvoll, Y. Méheust, R. Toussaint, J. Schmittbuhl, and K. J. Måløy, “Growth activity during fingering in a porous hele-shaw cell,” Phys. Rev. E 70, 026301 (2004).
  * Chuoke, van Meurs, and van der Poel [1959] R. L. Chuoke, P. van Meurs, and C. van der Poel, “The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media,” Transactions of the AIME 216, 188 (1959).
  * Yortsos and Huang [1986] Y. C. Yortsos and A. B. Huang, “Linear-stability analysis of immiscible displacement: Part 1 – simple basic flow profiles,” SPE Res. Eng. 1, 378 (1986).
  * van Batenburg  _et al._ [2015] D. W. van Batenburg, S. Berg, S. Oedai, L. L. David, A. O. N. Siemens, and K. Elewaut, “Visualisation of light oil mobilisation in asp core floods using x-ray ct imaging,” Soc. Petro. Eng. SPE-174660-MS (2015), 10.2118/174660-MS.
  * Meakin and Deutch [1986] P. Meakin and J. M. Deutch, “The formation of surfaces by diffusion limited annihilation,” J. Chem. Phys. 85, 2320 (1986).
  * Frette  _et al._ [1997] O. I. Frette, K. J. Måløy, J. Schmittbuhl, and A. Hansen, “Immiscible displacement of viscosity-matched fluids in two-dimensional porous media,” Phys. Rev. E 55, 2969 (1997).
  * Perkins, Johnston, and Hoffman [1965] T. K. Perkins, O. C. Johnston, and R. N. Hoffman, “Mechanics of viscous fingering in miscible systems,” SPE J. 5, 301 (1965).
  * Paterson [1985] L. Paterson, “Fingering with miscible fluids in a hele shaw cell,” Physics of Fluids 28, 26 (1985).
  * Jha, Cueto-Felgueroso, and Juanes [2011] B. Jha, L. Cueto-Felgueroso, and R. Juanes, “Fluid mixing from viscous fingering,” Phys. Rev. Lett. 106, 194502 (2011).
  * Fu, Cueto-Felgueroso, and Juanes [2017] X. Fu, L. Cueto-Felgueroso, and R. Juanes, “Viscous fingering with partially miscible fluids,” Phys. Rev. Fluids 2, 104001 (2017).
  * Salmo  _et al._ [2022] I. C. Salmo, K. S. Sorbie, A. Skauge, and M. A. Alzaabi, “Immiscible viscous fingering: Modelling unstable water-oil displacement experiments in porous media,” Transp. Porous Med. 145, 291 (2022).
  * Blunt [1998] M. J. Blunt, “Physically-based network modeling of multiphase flow in intermediate-wet porous media,” J. Petro. Sc. Eng. 20, 117 (1998).
  * Feder, Flekkøy, and Hansen [2022] J. Feder, E. G. Flekkøy, and A. Hansen,  _Physics of Flow in Porous Media_ (Cambridge University Press, Cambridge, 2022).
  * Pinilla, Asuaje, and Ratkovich [2021] A. Pinilla, M. Asuaje, and N. Ratkovich, “Experimental and computational advances on the study of viscous fingering: An umbrella review,” Heliyon 7, e07614 (2021).
  * Shahsavar, Fu, and Zhao [2023] N. Shahsavar, X. Fu, and B. Zhao, “Fingering instability during mixing-driven precipitation flow,” Transp. Porous. Med. (2023), 10.1007/s11242-023-02003-6.
  * Sandnes  _et al._ [2011] B. Sandnes, E. G. Flekkøy, H. Knudsen, K. J. Måløy, and H. See, “Patterns and flow in frictional fluid dynamics,” Nat. Comm. 2, 288 (2011).
  * Zhang  _et al._ [2023] D. Zhang, J. M. Campbell, J. A. Eriksen, E. G. Flekkøy, K. J. Måløy, C. W. MacMinn, and B. Sandnes, “Frictional fluid instabilities shaped by viscous forces,” Nat. Comm. 14, 3044 (2023).
  * Eriksen  _et al._ [2022] F. K. Eriksen, M. Moura, M. Jankov, A. L. Turquet, and K. J. Måløy, “Transition from viscous fingers to compact displacement during unstable drainage in porous media,” Phys. Rev. Fluids 7, 013901 (2022).
  * Lanza  _et al._ [2023] F. Lanza, S. Sinha, A. Hansen, A. Rosso, and L. Talon, “Transition from viscous fingers to foam during drainage in heterogeneous porous media,” ArXiv 2307.13451 (2023), 10.48550/arXiv.2307.13451.
  * McLean and Saffman [1981] J. W. McLean and P. G. Saffman, “The effect of surface tension on the shape of fingers in a hele shaw cell,” J. Fluid Mech. 102, 455 (1981).
  * Feder [1988] J. Feder,  _Fractals_ (Springer New York, NY, 1988).
  * Lenormand, Touboul, and Zarcone [1988] R. Lenormand, E. Touboul, and Zarcone, “Numerical models and experiments on immiscible displacements in porous media,” J. Fluid Mech. 189, 165 (1988).
  * Lenormand and Zarcone [1989] R. Lenormand and C. Zarcone, “Capillary fingering: Percolation and fractal dimension,” Transp. Porous Med. 4, 599 (1989).
  * Guo and Aryana [2019] F. Guo and S. A. Aryana, “An experimental investigation of flow regimes in imbibition and drainage using a microfluidic platform,” Energies 12, 1390 (2019).
  * Wilkinson and Willemsen [1983] D. Wilkinson and J. F. Willemsen, “Invasion percolation: a new form of percolation theory,” J. Phys. A: Math. Gen. 16, 3365 (1983).
  * Homsy [1987] G. M. Homsy, “Viscous fingering in porous media,” Annu. Rev. Fluid Mech. 19, 271 (1987).
  * Måløy, Feder, and Jøssang [1985] K. J. Måløy, J. Feder, and T. Jøssang, “Viscous fingering fractals in porous media,” Phys. Rev. Lett. 55, 2688 (1985).
  * Arnéodo  _et al._ [1989] A. Arnéodo, Y. Couder, G. Grasseau, V. Hakim, and M. Rabaud, “Uncovering the analytical saffman-taylor finger in unstable viscous fingering and diffusion-limited aggregation,” Phys. Rev. Lett. 63, 984 (1989).
  * Witten and Sander [1981] T. A. J. Witten and L. M. Sander, “Diffusion-limited aggregation, a kinetic critical phenomenon,” Phys. Rev. E. 47, 1400 (1981).
  * Witten and Sander [1983] T. A. J. Witten and L. M. Sander, “Diffusion-limited aggregation,” Phys. Rev. B. 27, 5686 (1983).
  * Darcy [1856] H. Darcy, Les Fontaines publiques de la ville de Dijon 647, Paris: Victor Dalamont (1856).
  * Whitaker [1986] S. Whitaker, “Flow in porous media i: A theoretical derivation of darcy’s law,” Transp. Porous Med. 1, 3 (1986).
  * Tallakstad  _et al._ [2009a] K. T. Tallakstad, H. A. Knudsen, T. Ramstad, G. Løvoll, K. J. Måløy, R. Toussaint, and E. G. Flekkøy, “Steady-state two-phase flow in porous media: Statistics and transport properties,” Phys. Rev. Lett. 102, 074502 (2009a).
  * Tallakstad  _et al._ [2009b] K. T. Tallakstad, G. Løvoll, H. A. Knudsen, T. Ramstad, E. G. Flekkøy, and K. J. Måløy, “Steady-state, simultaneous two-phase flow in porous media: An experimental study,” Phys. Rev. E. 80, 036308 (2009b).
  * Aursjø  _et al._ [2014] O. Aursjø, M. Erpelding, K. T. Tallakstad, E. G. Flekkøy, A. Hansen, and K. J. Måløy, “Film flow dominated simultaneous flow of two viscous incompressible fluids through a porous medium,” Front. Phys. 2, 63 (2014).
  * Rassi, Codd, and Seymour [2011] E. M. Rassi, S. L. Codd, and J. D. Seymour, “Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow,” New J. Phys. 13, 015007 (2011).
  * Rassi, Codd, and Seymour [2014] E. M. Rassi, S. L. Codd, and J. D. Seymour, “Corrigendum: Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow,” New J. Phys. 16, 039501 (2014).
  * Sinha  _et al._ [2017] S. Sinha, A. T. Bender, M. Danczyk, K. Keepseagle, C. A. Prather, J. M. Bray, L. W. Thrane, J. D. Seymour, S. L. Codd, and A. Hansen, “Effective rheology of two-phase flow in three-dimensional porous media: experiment and simulation,” Transp. Porous Med. 119, 77 (2017).
  * Gao  _et al._ [2020] Y. Gao, Q. Lin, B. Bijeljic, and M. J. Blunt, “Pore-scale dynamics and the multiphase darcy law,” Phys. Rev. Fluids. 5, 013801 (2020).
  * Zhang  _et al._ [2021] Y. Zhang, B. Bijeljic, Y. Gao, Q. Lin, and M. J. Blunt, “Quantification of non-linear multiphase flow in porous media,” Geophys. Res. Lett. 48, e2020GL090477 (2021).
  * Zhang, Bijeljic, and Blunt [2022] Y. Zhang, B. Bijeljic, and M. J. Blunt, “Nonlinear multiphase flow in hydrophobic porous media,” J. Fluid Mech. 934, R3 (2022).
  * Sinha and Hansen [2012] S. Sinha and A. Hansen, “Effective rheology of immiscible two-phase flow in porous media,” Europhys. Lett. 99, 44004 (2012).
  * Sinha  _et al._ [2013] S. Sinha, A. Hansen, D. Bedeaux, and S. Kjelstrup, “Effective rheology of bubbles moving in a capillary tube,” Phys. Rev. E 87, 025001 (2013).
  * Cheon  _et al._ [2023] H. L. Cheon, H. Fyhn, A. Hansen, Ø. Wilhelmsen, and S. Sinha, “Steady-state two-phase flow of compressible and incompressible fluids in a capillary tube of varying radius,” Transp. Porous Med. 147, 15 (2023).
  * Roy, Hansen, and Sinha [2019] S. Roy, A. Hansen, and S. Sinha, “Effective rheology of two-phase flow in a capillary fiber bundle model,” Front. Phys. 7, 92 (2019).
  * Yiotis, Talon, and Salin [2013] A. G. Yiotis, L. Talon, and D. Salin, “Blob population dynamics during immiscible two-phase flows in reconstructed porous media,” Phys. Rev. E. 87, 033001 (2013).
  * Roux and Herrmann [1987] S. Roux and H. J. Herrmann, “Disorder-induced nonlinear conductivity,” Europhys. Lett. 1227, 4 (1987).
  * Roy, Sinha, and Hansen [2021] S. Roy, S. Sinha, and A. Hansen, “Role of pore-size distribution on effective rheology of two-phase flow in porous media,” Front. Water 3, 709833 (2021).
  * Fyhn  _et al._ [2021] H. Fyhn, S. Sinha, S. Roy, and A. Hansen, “Rheology of immiscible two-phase flow in mixed wet porous media: Dynamic pore network model and capillary fiber bundle model results,” Transp. Porous Med. 139, 491 (2021).
  * Fyhn, Sinha, and Hansen [2023a] H. Fyhn, S. Sinha, and A. Hansen, “Effective rheology of immiscible two-phase flow in porous media consisting of random mixtures of grains having two types of wetting properties,” Front. Phys. 11, 1175426 (2023a).
  * Toussaint  _et al._ [2005] R. Toussaint, G. Løvoll, Y. Méheust, K. J. Måløy, and J. Schmittbuhl, “Influence of pore-scale disorder on viscous fingering during drainage,” Europhys. Lett. 71, 583 (2005).
  * Arneodo  _et al._ [1991] A. Arneodo, F. Argoul, Y. Couder, and M. Rabaud, “Anisotropic laplacian growths: From diffusion-limited aggregates to dendritic fractals,” Phys. Rev. Lett. 66, 2332 (1991).
  * Arneodo  _et al._ [1996] A. Arneodo, J. Elezgaray, M. Tabard, and F. Tallet, “Statistical analysis of off-lattice diffusion-limited aggregates in channel and sector geometries,” Phys. Rev. E 53, 6200 (1996).
  * Tao, Novotny, and Kaski [1988] R. Tao, M. A. Novotny, and K. Kaski, “Diffusion-limited aggregation with surface tension,” Phys. Rev. A 38, 1019 (1988).
  * Novotny, Tao, and Landau [1990] M. A. Novotny, R. Tao, and D. P. Landau, “Relaxation in dla with surface tension,” J. Phys. A: Math. Gen. 23, 3271 (1990).
  * Niemeyer, Pietronero, and Wiesmann [1984] L. Niemeyer, L. Pietronero, and H. J. Wiesmann, “Fractal dimension of dielectric breakdown,” Phys. Rev. Lett. 52, 1033 (1984).
  * Amitrano [1989] C. Amitrano, “Fractal dimensionality for the $\eta$ model,” Phys. Rev. A 39, 6618 (1989).
  * Mathiesen and Jensen [2002] J. Mathiesen and M. H. Jensen, “Tip splittings and phase transitions in the dielectric breakdown model: Mapping to the diffusion-limited aggregation model,” Phys. Rev. Lett. 88, 235505 (2002).
  * Sinha  _et al._ [2021] S. Sinha, M. A. Gjennestad, M. Vassvik, and A. Hansen, “Fluid meniscus algorithms for dynamic pore-network modeling of immiscible two-phase flow in porous media,” Front. Phys. 9, 548497 (2021).
  * Aker  _et al._ [1998] E. Aker, K. J. Maloy, A. Hansen, and G. G. Batrouni, “A two-dimensional network simulator for two-phase flow in porous media,” Transp. Porous Med. 32, 163 (1998).
  * Erpelding  _et al._ [2013] M. Erpelding, S. Sinha, K. T. Tallakstad, A. Hansen, E. G. Flekkøy, and K. J. Måløy, “History independence of steady state in simultaneous two-phase flow through two-dimensional porous media,” Phys. Rev. E 88, 053004 (2013).
  * Roy, Sinha, and Hansen [2020] S. Roy, S. Sinha, and H. Hansen, “Flow-area relations in immiscible two-phase flow in porous media,” Front. Phys. 8, 4 (2020).
  * Roy  _et al._ [2022] S. Roy, H. Pedersen, S. Sinha, and A. Hansen, “The co-moving velocity in immiscible two-phase flow in porous media,” Transp. Porous. Med. 143, 69 (2022).
  * Fyhn, Sinha, and Hansen [2023b] H. Fyhn, S. Sinha, and A. Hansen, “Local statistics of immiscible and incompressible two-phase flow in porous media,” Physica A 616, 128626 (2023b).
  * Sinha  _et al._ [2019] S. Sinha, M. Aa. Gjennestad, M. Vassvik, M. Winkler, A. Hansen, and E. G. Flekkøy, “Rheology of high-capillary number two-phase flow in porous media,” Front. Phys. 7, 65 (2019).
  * Dullien [1992] F. A. L. Dullien,  _Porous media: Fluid, transport and pore structure_ (Academic Press, San Diego, 1992).
  * Washburn [1921] E. W. Washburn, “The dynamics of capillary flow,” Phys. Rev. 17, 273 (1921).
  * Batrouni and Hansen [1988] G. G. Batrouni and A. Hansen, “Fourier acceleration of iterative processes in disordered systems,” J. Stat. Phys. 52, 747 (1988).
  * Løvoll  _et al._ [2005] G. Løvoll, Y. Méheust, K. J. Måløy, E. Aker, and J. Schmittbuhl, “Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study,” Energy 30, 861 (2005).
  * Somfai and Ball [2003] E. Somfai and R. C. Ball, “Diffusion-limited aggregation in channel geometry,” Phys. Rev. E 68, 020401(R) (2003).
  * Couder  _et al._ [1990] Y. Couder, F. Argoul, A. Arnéodo, J. Maurer, and M. Rabaud, “Statistical properties of fractal dendrites and anisotropic diffusion-limited aggregates,” Phys. Rev. A 42, 3499 (1990).
  * Combescot  _et al._ [1986] R. Combescot, T. Dombre, V. Hakim, Y. Pomeau, and A. Pumir, “Shape selection of saffman-taylor fingers,” Phys. Rev. Lett. 56, 2036 (1986).
  * Vanden-Broeck [1983] J.-M. Vanden-Broeck, “Fingers in a hele-shaw cell with surface tension,” Phys. Fluids 26, 2033 (1983).
  * Thomé  _et al._ [1989] H. Thomé, M. Rabaud, V. Hakim, and Y. Couder, “The saffman-taylor instability: From the linear to the circular geometry,” Phys. Fluids 1, 224 (1989).