# A simplified pore-scale model for slow drainage including film flow effects

Paula Reis111PoreLab, The NJORD Centre, Department of Physics, University of
Oslo, Norway Marcel Moura222PoreLab, The NJORD Centre, Department of Physics,
University of Oslo, Norway Gaute Linga333PoreLab, The NJORD Centre, Department
of Physics, University of Oslo, Norway Per Arne Rikvold444PoreLab, The NJORD
Centre, Department of Physics, University of Oslo, Norway Eirik Grude
Flekkøy555PoreLab, The NJORD Centre, Department of Physics, University of
Oslo, Norway Knut Jørgen Måloy666PoreLab, The NJORD Centre, Department of
Physics, University of Oslo, Norway

###### Abstract

Film flow through networks of corners and capillary bridges can establish
connections between seemingly isolated clusters during drainage in porous
media. Coupled with drainage through the bulk of pores and throats, the flow
through these networks constitutes a secondary drainage mechanism that can
significantly affect fluid configurations and residual saturations. In order
to investigate the prevalence of this drainage mechanism, we propose a quasi-
static pore-network model based on modifying the trapped-cluster-
identification algorithm in an invasion-percolation model. With the
modification, wetting-phase connectivity could be provided by direct
successions of pores and throats, represented by sites and bonds, as well as
by chains of interconnected capillary bridges. The advancement of the fluid
interface in the porous matrix was determined by the bonds’ invasion
thresholds, calculated based on local capillary pressure values that could be
perturbed to accommodate gravitational effects. With the proposed model,
experimentally verified phenomena related to slow drainage in granular porous
media were reproduced, showing good qualitative agreement.

††journal: arXiv

##  1 Introduction

Pore-scale modeling of fluid-fluid displacements has been an active research
topic for the past few decades, with applications spanning multiple fields [1,
2, 3, 4, 5]. In particular, extensive efforts have been directed toward the
development of models related to $CO_{2}$ storage in subsurface geological
reservoirs [6, 7, 8, 9, 10], oil recovery methods [11, 12, 13, 14, 15], drying
[16, 17, 18, 19, 20], and water management in fuel cells [21, 22, 23, 24, 25].
The quality of predictions from such models relies on appropriately
identifying and representing the primary forces driving the flow, as well as
the fundamental geometrical aspects of the porous space, which can vary
greatly according to the system of interest. In this work, we narrow our focus
to effectively modeling the impact of film flow during slow drainage in
granular porous media. With this aim, first, we identify the main flow
mechanisms exhibited during slow drainage and how their prevalence can be
affected by the particular geometrical features of pores within granular
materials. Then, we propose an approach to incorporate these mechanisms’
essential effects on a computationally cost-effective pore-scale model.

Slow drainage flows are characterized by the displacement of a wetting fluid
by a non-wetting fluid under no significant influence of viscous or inertial
forces. During such flows, two main mechanisms are commonly identified,
namely, a piston-like invasion of the non-wetting phase, followed by the
removal of the wetting phase through networks of contiguous corners and/or
capillary bridges [26, 27, 28]. The piston-like displacement is characterized
by the advancement of the fluid interface through the bulk of pore bodies and
throats in fast discrete events known as Haines jumps [29, 30, 31, 32]. As
this mechanism unfolds, the competition between local capillary pressure
($P_{c}$) and pore entry pressure thresholds ($P_{t}$), shown in eq.1,
dictates the advancement of the invading fluid front. This mechanism is
generally responsible for most of the displacement of the defending fluid and
will be referred to as the primary drainage mechanism.

$P_{c}>P_{t}=\gamma\cos{\theta}\left(\frac{1}{r_{1}}+\frac{1}{r_{2}}\right)$
(1)

where $\gamma$ is the interfacial tension between the fluids, $\theta$ is the
contact angle between the fluids and the porous material, and $r_{1}$ and
$r_{2}$ are the two main interface curvature radii.

In the regions of the porous medium swept by the invading fluid front, the
primary mechanism is succeeded by the displacement of the wetting fluid
through film flow [33, 34]. At this stage, the wetting fluid can be found in
the medium either in the form of isolated clusters – structures containing one
or several defending-phase-filled pores surrounded by the invading phase – or
at angular porous features passed by the interface, where $P_{c}$ was
insufficient to mobilize the defending phase. In this unsaturated zone, the
secondary drainage mechanism takes place, provided that contiguous sets of
corners and capillary bridges constitute paths for wetting fluid flow [35, 36,
37, 28]. Given the comparatively lower permeability of the corner/bridge flow
paths, the secondary mechanism is characterized by distinctive longer time
scales for drainage than the primary mechanism[26, 28].

The efficiency of both drainage mechanisms is inherently tied to geometric
aspects of the porous space [38, 39, 40]. During piston-like displacement, the
pore’s size distribution and spatial arrangement define the pore-invasion
order, delineating the fraction of the porous medium bypassed by the front,
i.e., the amount of wetting fluid retained in clusters. As for the secondary
mechanism, the shape of the porous space – along with wettability – is the key
element dictating the hydraulic conductance and connectivity of the wetting-
phase film flow paths established in the medium. In the following passage, a
brief review of wetting phase connectivity in granular materials is presented.

Within unsaturated granular porous media, the wetting liquid not pertaining
clusters is mostly retained in the form of capillary bridges and rings [41,
42, 43, 44, 45, 28] enclosing grain-grain contact areas. Despite their minor
contribution to the wetting-phase saturation, these structures can have a
significant impact on the transport properties, given that the conditions for
their interconnectivity are met. In a study related to drying in granular
media, Cejas et al. [43] investigated the role of grain packing in the
stability of film networks able to transport water from the bulk of the porous
medium to its surface. Their results pointed to the existence of a packing
threshold in loosely packed systems where the liquid connections break and
water removal from the porous medium is undermined. In a similar set of
studies, Chen et al. [44][45] investigated the effect of geometry in the
formation of elongated liquid films using transparent quasi-2D micromodels of
cylinders arranged between two planes. According to the spatial disposition of
the cylinders, the authors could observe the formation of films consisting of
chains of capillary bridges providing hydraulic connectivity between internal
liquid clusters and the outer rim of the micromodels. Quasi-2D models were
also adopted by Moura et al. [28] to evaluate the impact of film flow on slow
drainage under varying degrees of gravitational force influence. Their
experiments using modified Hele-Shaw cells filled with a single layer of
monodispersed glass beads allowed the analysis of several phenomena related to
film flow, such as its relative impact on the residual saturation, the size
and spatial distributions of capillary bridges, and the establishment of a
film-flow active zone trailing the invasion front. Interestingly, the results
indicated the drainage of seemingly isolated wetting-fluid clusters through
dynamic networks of capillary bridges non-coincident with the pore network
invaded during piston-like displacement.

Although significant experimental evidence points to the critical role of
capillary-bridge-chain flow paths during drainage in porous media, few
computational models represent their effect on the flow explicitly. A
pioneering model in this direction was presented by Vorhauer et al. [46], with
the intent of predicting isothermal drying rates in porous media. Based on
experimental observations of drying in etched-silicon micromodels, the authors
proposed a 2D pore-network model (PNM) that explicitly takes into account the
role of capillary bridges and rings in wetting-fluid-cluster connectivity.
With this new modeling approach, better agreement with experimental results
was obtained, when compared with classic PNMs that represent film-flow paths
as interconnected corners of angular-shaped pores. An extension of this work
was later presented by Kharaghani et al. [47], by changing the network
topology from 2D regular square lattices to 3D regular cubic lattices. Drying
rates predicted by the 3D PNM showed good quantitative agreement with
experiments using monodisperse-spherical-glass-bead packs initially filled
with distilled water. Another model taking into account the role of capillary
bridges during drying in porous media was developed by Chen et al. [45],
accompanying their experiments on quasi-2D porous media formed by cylinders
arranged between two parallel planes. In their model, liquid films consisted
of successions of capillary bridges connecting neighboring cylinders, within
which flow dynamics were shaped by the effects of viscous and capillary
forces. Using experimental drying rates as model inputs, their numerical
results were able to qualitatively predict invasion patterns and the extent of
films as the water evaporated from the porous space.

In the present work, we incorporate the capillary-bridge role in wetting-phase
connectivity successfully explored in drying models on an invasion-percolation
model for slow drainage under the effects of capillary and gravitational
forces. Driven by the experimental findings presented by Moura et al. [28], we
use the proposed PNM to investigate the capillary bridge size and spatial
distributions, the contribution of film flow to residual saturations, and the
existence of film flow active zones on quasi-2D porous media.

##  2 Methodology

The proposed model is based on the invasion-percolation (IP) method with
trapping [48]. Designed to represent immiscible fluid-fluid displacements in
porous media, this method involves the description of pores and throats as
lattices of sites and bonds, to each of which an invasion threshold is
assigned. Flow in the idealized porous lattices is then predicted based on the
movement of the fluid interface through the least resistivity path, dictated
by the bonds’ invasion thresholds in the case of drainage, and the sites’
invasion thresholds in the case of imbibition [48]. This modeling approach has
been widely adopted in the past decades [49] due to its ability to sensibly
predict porous media flows in the capillary-dominated regime associated with
low computational effort.

In the direction of investigating fundamental aspects of film flow through
corners and capillary bridges in granular porous media, our modified IP model
was specially tailored to reproduce slow drainage in Hele-Shaw cells filled
with a single layer of monodispersed spherical beads. This choice of quasi-2D
model porous media allowed us to clearly identify drainage events related to
film flow and generate results that could be readily compared with available
experimental data [28].

Insights related to fluid-fluid displacements obtained via visualization of
flow in transparent quasi-2D models can be generalized to more complex media,
by acknowledging how inherently 3D features of such models affect the observed
flow phenomena. For the study of flow through films formed by corners and
capillary bridges, an important geometric feature of quasi-2D models is the
intersection of the structures representing the porous media grains, e.g.
cylinders [50, 44, 51, 52] or spheres [53, 54, 55, 28], with the planes
containing them. Following the piston-like displacement of the interface,
fluid collected in these regions can play a fundamental role in wetting-phase
continuity by linking capillary bridges formed between pairs of grains, as
indicated in Figure 1. Several recent pore-scale models incorporate aspects of
these fluid configurations stemming from 3D structures in 2D lattices
representing the porous space [45, 56, 57, 58] and we adopt a similar
approach.

Figure 1: Examples of wetting-fluid configurations around pairs of cylinders
and spheres contained between parallel planes, obtained with Surface Evolver
[59]. The blue, magenta, and gray surfaces represent the portions of fluid in
contact with the planes, the structures representing grains, and the non-
wetting fluid, respectively.

In the next sections, we introduce the network representation of the porous
space, the criteria for fluid interface advancement, and the role of
capillary-bridge flow paths in fluid connectivity.

###  2.1 Porous-media representation

The description of porous media as networks of nodes (pore bodies) connected
by edges (pore throats) is a common approach in the development of
computationally-efficient models for multiphase flow [49, 60]. In this study,
we propose a 2D regular lattice representation of the quasi-2D porous media
designed for the experimental investigation of film flow effects during
drainage performed by Moura et al. [28]. Their porous-media models were
created by randomly placing and securing spherical glass beads with a diameter
of $1mm$ between two parallel transparent plates. The resulting porous
matrices could be divided into pore bodies and pore throats located at the
vertices and edges of the Voronoi diagrams created using the center of the
beads as seeds. Considering dense homogeneous packings of spherical beads,
these porous matrices could be approximated as regular honeycomb lattices, as
indicated in Figure 2.

Figure 2: Sketch of the modified Hele-Shaw porous matrix used in [28] and its
representation by a regular honeycomb lattice of nodes (open circles) and
edges (solid lines).

In the lattice representation, the edges portray throats with both length and
height equal to the diameter of the beads, and with width related to the
distance between the beads adjacent to them. The choice of width values
assigned to the edges followed the probability density function of beads’
separation found in the experimental models. We did not address the
geometrical aspects of the lattice nodes and assumed that pores do not exert
significant resistance to the drainage flows under investigation.

###  2.2 Modified IP Algorithm

The proposed model for slow drainage can be described as an invasion-
percolation model with trapping, in which the defending cluster labeling
algorithm was modified to allow connectivity via capillary bridges. Its
implementation exploits some basic concepts of graph theory, which will be
briefly addressed next. A graph is composed of a set of nodes and a set of
edges. An edge connects a pair of nodes and can be directed, if the pair of
nodes is ordered, or undirected, otherwise. Pairs of nodes connected by an
edge are called adjacent, and the edge is incident on the pair of nodes it
connects. The number of edges incident on a node determines the node’s degree.
A sequence of adjacent nodes in a graph defines a path. A group of nodes is
connected if there exists a path between every two nodes, and is referred to
as a connected component of a graph.

By using these concepts, the model’s main idea is to create an undirected
graph in which the connected components represent the wetting fluid occupation
of the porous medium. As the drainage simulations start with a fully saturated
porous matrix, this graph’s topology initially matches the honeycomb lattice
presented in Section 2.1. The nodes corresponding to the medium’s source of
non-wetting fluid and sink of wetting fluid are identified and labeled as the
inlet and outlet nodes, respectively. Weights are assigned to the graph’s
edges, corresponding to the capillary pressure threshold for the invasion of
the pore throats ($P_{t}$). As the model takes into account the effects of
capillary and gravitational forces, the calculation of $P_{t}$ is given by eq.
2

$P_{t}=2\gamma\cos{\theta}\left(\frac{1}{t_{w}}+\frac{1}{t_{h}}\right)-\Delta\rho
gh$ (2)

where $t_{w}$ and $t_{h}$ are the throat’s width and height, respectively,
$\Delta\rho$ is the density difference between wetting and non-wetting fluids,
$g$ is the gravitational acceleration and $h$ is the elevation difference
between the throat and the wetting-phase sink.

As with IP models, our drainage simulation unfolds in discrete steps, in which
the invasion of the available bond with the lowest $P_{t}$ takes place. This
invasion is expressed by deleting the available edge with the lowest weight
from the graph representing the wetting fluid occupation of the porous medium,
as indicated in Figure 3. At the beginning of the drainage simulation, only
the edges incident on the inlet nodes are considered available for invasion
(Figure 3a). As an edge is deleted, the degree of the nodes on which it was
incident is decreased, and these nodes’ remaining edges become available for
invasion (Figure 3b-d). The decrement of a node’s degree can be understood as
the change in its occupation from wetting to non-wetting fluid. During the
edge-deleting process, the graph can be split into multiple connected
components. Connected components that do not contain outlet nodes represent
wetting-fluid trapped clusters, and their edges become unavailable for
invasion (Figure 3e). Following these rules, the drainage process continues in
the model until all available edges are deleted.

Figure 3: Representation of the modified IP model algorithm for drainage

####  2.2.1 Cluster Merging via Capillary Bridges

The drainage process illustrated in Figure 3 does not comprise the formation
of capillary bridges and could be reproduced with a regular bond invasion
percolation model. The cluster-labeling modification we propose in our model
comes into being when capillary bridges are formed in a way that continuity is
established between otherwise disconnected clusters, as suggested in Figure 4.

Figure 4: Sketch of Hele-Shaw porous-space representation into a 2D regular
honeycomb lattice of nodes and edges.

In the upper left corner of Figure 4, a macro picture of a Hele-Shaw cell
under drainage from [28] is shown. In this insert, the two arrows highlight
chains formed by capillary bridges through which fluid transport is
permissible. In a regular IP model, this type of connectivity is not taken
into account, as no paths are formed between the nodes at the endpoints of
capillary-bridge chains. This is indicated in the sketch on the right side of
Figure 4: capillary bridges appear in the graph as edges incident on two nodes
filled with the non-wetting phase, and the connection they provide for the
wetting fluid does not coincide with the network representing the porous
space.

The capillary-bridge chain connectivity was then incorporated into our model
by considering that edges belonging to the same hexagon in the honeycomb
lattice are part of the same connected component, even if they do not directly
establish a path. This rule is equivalent to acknowledging that the wetting
fluid residing in the contact points between the spheres and the planes in the
Hele-Shaw cell links all bridges formed around the same spheres, as
illustrated in Figure 1. With this cluster-labeling modification, our model
could reproduce experimentally verified aspects of drainage events related to
both primary and secondary mechanisms, as presented in Section 3

##  3 Results and Discussion

The results presented in this section were obtained with drainage simulations
on 100 random realizations of honeycomb lattices with 30703 nodes and 45903
edges, representing quasi-2D porous media of approximately 15x26 cm. The
capillary pressure threshold values of the edges representing throats followed
the distribution presented in Figure 5, based on the experiments presented by
Moura et al. [28]. These values were calculated with the separation between
the spheres reported in [28] – which varied approximately from 0.15 to 1.5 mm
–, an interfacial tension between wetting and non-wetting fluids of
$\gamma=64mN/m$, a density contrast between the fluids of
$\Delta\rho=1203.75kg/m^{3}$, and a contact angle of $\theta=0^{\circ}$
between the fluids and the medium. The influence of gravitational forces on
the drainage was evaluated by adopting the following values of inclination of
the porous media with respect to the horizontal plane:
$\beta=[0,1.25,2.5,5,7.5,10,15,20,30,45,60,75,90,]^{\circ}$.

Figure 5: Probability density function of the capillary pressure of the porous
matrix throats, generated to match the values reported in Moura et al. [28].

In the lattices, lateral boundaries ($x=0cm$ and $x=15cm$) were closed to the
flow, nodes at the upper boundary ($y=26cm$) were sources of non-wetting
fluid, and all nodes at the lower boundary ($y=0cm$) were connected to an
external wetting-fluid sink by edges with no resistance to flow. In this way,
all drainage simulations came to an end by the time the non-wetting phase
percolated the lattice.

Figure 6 illustrates the shift in drainage dynamics predicted by the proposed
model as the influence of the gravitational component increases. These results
were generated using the same lattice and with inclination angles indicated in
the captions.

(a) $\beta=0^{\circ}$

(b) $\beta=1.25^{\circ}$

(c) $\beta=2.5^{\circ}$

(d) $\beta=7.5^{\circ}$

(e) $\beta=10^{\circ}$

(f) $\beta=20^{\circ}$

(g) $\beta=30^{\circ}$

(h) $\beta=60^{\circ}$

Figure 6: Representation of the invasion order in the porous matrix during
drainage, from dark blue (early invasion) to dark red (late invasion).

From 6(a) to 6(h), only the nodes and edges occupied by the non-wetting phase
at breakthrough are represented. The color gradient from dark blue to dark red
symbolizes the order of invasion, from early to late invaded pores and
throats. Black regions of the image are, therefore, correspondent to the
trapped wetting fluid clusters. As expected, the unstable front characteristic
of capillary fingering is progressively stabilized as $\beta$ increases,
resulting in a more efficient removal of the wetting phase. The invasion
dynamics presented in 6 are a result of both primary and secondary drainage
mechanisms, which will be discussed in more detail in the next sections.

###  3.1 The impact of film flow on residual saturation

As presented in Section 1, chains of capillary bridges can establish
connectivity between seemingly trapped clusters and enable their drainage,
leading to a reduction in wetting-fluid residual saturations. The prevalence
of primary and secondary drainage mechanisms for different degrees of
gravitational forces influence is illustrated in Figure 7.

(a) $\beta=0^{\circ}$

(b) $\beta=1.25^{\circ}$

(c) $\beta=2.5^{\circ}$

(d) $\beta=7.5^{\circ}$

(e) $\beta=10^{\circ}$

(f) $\beta=20^{\circ}$

(g) $\beta=30^{\circ}$

(h) $\beta=60^{\circ}$

Figure 7: Representation of the drainage mechanism related to invasion in the
porous matrix. Nodes and edges in blue were drained via the primary drainage
mechanism, while nodes and edges in yellow were drained by the secondary
mechanism.

Complementary to the information presented in Figure 6, the images in Figure 7
introduce the distinction between drained pores and throats that pertained
directly to the main defending cluster (blue), to the ones that were connected
to it via capillary bridges (yellow). It is clear that a significant fraction
of the wetting fluid in the porous matrix could only be dislocated through the
flow paths consisting of capillary bridges and corners, especially for
intermediate inclination angles. The quantification of the contribution of
both drainage mechanisms averaged over the different lattice realizations is
presented in Figure 8.

Figure 8: Contribution of primary and secondary drainage mechanisms.

As indicated in the graph, the fraction of wetting fluid displaced due to the
primary drainage mechanism increases monotonically as the influence of
gravitational forces grows, while the maximum saturation drained due to the
secondary mechanism is observed at $\beta=7.5^{\circ}$. As for the incremental
drainage provided by the secondary drainage, the maximum mean value is
achieved at $\beta=2.5^{\circ}$. This behavior may be explained due to the
contradicting effects gravity can have on drainage through film flow. On one
side, gravity can benefit flow from clusters placed above the average vertical
position of the drainage front. Boundaries of clusters bypassed by the front
are associated with higher capillary pressure thresholds for invasion than the
regions in their immediate vicinity. However, as these clusters get further to
the front, the negative contribution of the hydrostatic pressure to their
invasion threshold may place them in a preferential position for invasion,
when compared to nodes belonging to the front. This may justify the relatively
low value of saturation, $4.39\pm 1.54\%$, drained due to film flow at
$\beta=0^{\circ}$. Unlike in the other cases, the distance of a node from the
front without the influence of gravity would not positively impact its chance
of invasion, leading to trapping. On the other side, the positive influence of
gravity on the primary mechanism is so significant that little is left for the
secondary drainage at high values of $\beta$. In these cases, the efficient
sweep of the porous matrix by the invasion front may, therefore, reduce the
saturation drained due to film flow.

It is important to mention at this point that the absence of viscous effects
in the proposed model may affect the analysis of drainage through film flow
based on the data obtained with the proposed model. While we propose to
investigate slow drainage, which is associated with negligible influence of
viscous forces, there is a distinction between the effect of viscosity on the
flow through the bulk of pores and throats and on the flow through chains of
capillary bridges. The higher impact of viscous forces on the latter penalizes
the secondary drainage mechanism, and we expect that the contribution of film
flow to drainage presented here to be overestimated to some extent. The
incorporation of viscous effects in the model could also lead to a shift in
the values of $\beta$ associated with maximum drainage via film flow, as the
negative influence of viscous forces on this mechanism should be
proportionally more relevant for low inclination. Finally, viscous effects are
partially responsible for the different time scales associated with primary
and secondary drainage mechanisms identified experimentally in [28].The joint
analysis of Figures 6 and 7 allow us to identify some contrast in the invasion
time of clusters drained due to film flow and their surroundings, drained
through the bulk of pores and throats, but this distinction should be more
prominent. Still, the results shown in Figure 8 exhibit a reasonable agreement
with experimental data obtained by [28]. In that study, an increase of $8.3\%$
in non-wetting phase saturation was reported for a drainage experiment with
$\beta=20^{\circ}$, an up to $10\%$ in analyzed subsets of the porous matrix
for a case with $\beta=30^{\circ}$. The discrepancy between experimental and
simulated results is in line with the expected due to the lack of
representation of viscous effects in the model, and we intend to address this
shortcoming in a future work.

###  3.2 The film flow active zone

In this section, we analyze the formation region trailing the invasion front
where drainage events related to film flow are more likely to occur. First
identified experimentally by Moura et al. [28], this region was termed "active
zone", and its detailed characterization can strengthen our understanding of
the secondary drainage mechanism.

The extent of the active zone is directly related to the size and spatial
distribution of the portions of the porous matrix bypassed by the invasion
front. First, there must be a series of capillary bridges positioned in a
contiguous form in order to form a chain. This imposes a constraint on how
wide the active zone can be, as long sequences of connected capillary bridges
are less likely to occur than short ones. Second, the frequency and size of
wetting-phase clusters behind the front contribute to the width of the active
zone, as clusters connected to the front by capillary bridges can then lead to
the direct connection of further regions in the porous matrix. Thus, the
existence of large wetting-phase clusters can broaden the zone where secondary
drainage is possible, even if long chains of capillary bridges are not formed.

Figure 9 illustrated the nodes connected to the main defending cluster via
films, halfway through drainage, under different degrees of influence of
gravitational forces. There is a clear indication from these images that the
active zone gets progressively more compact as $\beta$ increases. For the
cases where the gravitational component is zero or relatively small,
$0^{\circ}\leq\beta\leq 2.5^{\circ}$, the frequency of large wetting fluid
clusters behind the front – as noticeable in Figures 6 and 7 – lead to very
extensive portions of the porous matrix connected to the front via films. As
the inclination of the porous matrix becomes steeper, there seems to be a
short range of angles where the active zone rapidly shrinks,
$2.5^{\circ}\leq\beta<10^{\circ}$, followed by a long inclination interval,
$10^{\circ}\leq\beta\leq 90^{\circ}$, where a narrow active zone is curtailed
slowly.

(a) $\beta=0^{\circ}$

(b) $\beta=1.25^{\circ}$

(c) $\beta=2.5^{\circ}$

(d) $\beta=7.5^{\circ}$

(e) $\beta=10^{\circ}$

(f) $\beta=20^{\circ}$

(g) $\beta=30^{\circ}$

(h) $\beta=60^{\circ}$

Figure 9: Representation of nodes connected to the front via films (pink),
nodes at the front (yellow), and nodes inside the main defending cluster
(blue) halfway through the drainage simulations.

The normalized probability density function (PDF) of distance from the front
to the nodes connected and drained via film flow is shown in Figure 10(a). As
suggested before, the probability of finding nodes pertaining to the active
zone closer to the front increases significantly as gravitational forces
become more effective during drainage. The same trend was observed
experimentally in [28] when comparing the PDF of the distance from the
invasion front of film flow events in experiments with
$\beta=10^{\circ},20^{\circ}$ and $30^{\circ}$. Interestingly, the normalized
PDFs presented in 10(a) and 10(b) are remarkably similar. This indicates that,
among all nodes that could be drained by film flow at a given moment, there is
no clear preference that further nodes would be invaded instead of nearer
ones.

(a)

(b)

Figure 10: Normalized PDF of the distance to the front of nodes connected via
film flow (a), and of the nodes effectively drained via film flow (b).

Another interesting fact is that wider active zones are not necessarily linked
to more significant effects of film flow events in drainage. Figure 11
indicates the average fraction of all nodes belonging to the lattice connected
to the main defending cluster via films, at any point during drainage. The
joint analysis of data presented in Figures 11 and 8 indicates that the large
amount of nodes belonging to the networks through which the secondary drainage
mechanism takes place at very low inclination angles do no translate into
greater transport of wetting fluid to the main defending cluster.

Figure 11: Average fraction of the total number nodes belonging to the lattice
connected to the invasion front during the drainage under different levels of
influence of gravitational forces

The characterization of the active zone using the size of the secondary-
drainage-mechanism networks should also take into account the fact that the
topology of such networks is under constant evolution as the porous matrix is
drained. This phenomenon was acknowledged in [28], when the authors described
the networks associated with primary drainage as "static" – formed by the
lattice of pores and throats representing the porous space –, and the ones
associated with film flow events as "dynamic" – formed by unsteady subsets of
the lattice following the invasion front. The extent of the fluctuation in the
dynamic networks’ sizes during drainage can be inferred by the error bars in
Figure 11, and an insight into their topological variations is presented in
Figure 12

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 12: Representation of nodes connected to the front via films for
$\beta=7.5^{\circ}$ at six stages during drainage

This figure illustrates the secondary drainage networks at eight equidistant
stages of a drainage simulation with $\beta=7.5^{\circ}$. We can observe that
even though their width is somewhat constant, their topology and size can vary
significantly. A reason pointed out in Moura et al. [28] for the frequent
substantial transformations in these dynamic networks is the potentially
global effect of single capillary bridges in the wetting-phase connectivity.
As shown in detail in Figure 13, which represents the region around the
invasion front in Figure 12(d), the five capillary bridges indicated by white
arrows are responsible for most of the communication between the dynamic film-
flow networks with the main defending cluster. Any drainage event in the
throats hosting these bridges, or in their immediate vicinity, could lead to
instant disconnection of large portions of the secondary drainage networks. In
order to expand our understanding of the likelihood of film flow events in a
given porous matrix, an analysis of the capillary bridge frequency in the used
lattices, as well as the probability of their occurrence given a throat’s size
is presented in the next Section.

Figure 13: Detail of the secondary drainage network shown in 12(d).
Illustration of the spherical beads surrounded by the networks is included for
clarity.

###  3.3 The frequency and size distribution of capillary bridges

Given the central role of capillary bridges in film flow events in granular
porous media, in this section we examine how likely it is for a capillary
bridge to be formed in the quasi-2D porous matrices under investigation, and
where they are mostly expected to occur. Figure 14(a) shows the fraction of
throats hosting capillary bridges, among all throats in the lattices, at the
moment of the invading fluid breakthrough. We can see that from
$\beta=0^{\circ}$ to $\beta=20^{\circ}$ there is a sharp increase in the
occurrence of bridges, which henceforth stabilizes at approximately $18.5\%$.
The comparison of these results with the ones presented in Figures 8 and 11
indicates that the pervasiveness of capillary bridges alone cannot be
correlated with a greater capacity of connecting, nor draining, extensive
regions in the lattice to the main defending cluster. In fact, the increase in
formation of capillary bridges at higher inclination angles can be linked to
the enhanced displacement of wetting-fluid via the primary drainage mechanism,
leaving smaller – and most likely disconnected – clusters behind the front. An
increase in the frequency of capillary bridges as inclination angles get
steeper was also observed experimentally by Moura et al. [28], when comparing
cases with $\beta=20^{\circ}$ and $30^{\circ}$. In that study, a fraction of
$15\%$ of throats occupied by bridges at $\beta=20^{\circ}$ was also reported,
which is approximately $17\%$ lower than the value predicted by the model.

(a)

(b)

Figure 14: Frequency of throats containing capillary bridges (a) and their
size distributions (b).

Complementary to the information presented in Figure 14(a), Figure 14(b)
presents the PDFs of the sizes of throats ($t_{w}$) in which capillary bridges
were formed during drainage at different inclination angles. In this graph we
can see that, even though the porous matrices comprised throats with sizes up
to 1.5 mm, all the bridges were found in throats smaller than 0.45 mm. It is
also noticeable that the same range of throat sizes hosting capillary bridges
was found for all tested $\beta$ values. As the inclination angles get
steeper, however, there is a tendency of the bridges to concentrate at the
lower end of this range.

Figure 15: Probability density function of the sizes of throats containing
capillary bridges experimentally and with the proposed model, for
$\beta=20^{\circ}$.

This heterogeneous occupation of throats by bridges was also identified
experimentally in [28], as shown in Figure 15 (where the number of counts of
the experimental results were rescaled to match the size of the lattices used
in the simulations). For the drainage experiment conducted at
$\beta=20^{\circ}$, the reported maximum size of throat hosting a bridge was
approximately 0.5mm, only $10\%$ larger than the size predicted by the model.
In that study, this maximum size was justified as a result of two main
factors: the Plateau-Rayleigh instability in the bridges – which would lead to
the snap-off of bridges accommodated between spheres set at wider separations
– and the nature of the drainage process – which leads to the invasion of
larger throats during the primary drainage mechanism. Given the low frequency
of snap-off events reported experimentally, and the fact that our model does
not incorporate a criterion for capillary bridge snap-off due to hydrodynamic
instabilities, we have a strong indication that the nature of the drainage
process is the predominant factor associated with the size of throats expected
to host bridges.

##  4 Conclusions

In this study, a simple quasi-static pore-network model for slow drainage in
granular porous media was presented. The model was based on the bond-invasion-
percolation method, and incorporated a modification in the trapped-cluster-
identification algorithm in order to acknowledge wetting-phase connectivity
provided by chains of capillary bridges. Porous media were represented in the
model as regular honeycomb lattices of nodes and edges, symbolizing the pores
and throats found in Hele-Shaw cells filled with a single layer of
monodispersed spherical beads. This choice of idealized quasi-2D porous media
allowed us to clearly characterize drainage events due to film flow, and the
obtained insights about the phenomenon can mostly be generalized to drainage
in 3D granular materials.

With the proposed model, the prevalence of primary and secondary drainage
mechanisms under different levels of influence of gravitational forces was
investigated, and the results were compared with experimental data from Moura
et al. [28]. Despite the model’s simple representation of the porous space and
drainage dynamics, qualitative agreement with experimental data was obtained
in the analyses regarding the impact of film flow on residual saturations, the
establishment of a film flow active zone and the occurrence of capillary
bridges in the medium. With the demonstrated ability to represent fundamental
aspects of complex physical phenomena, the proposed model could work as a
complementary tool to experiments, in the quest for a comprehensive
understanding of drainage in granular porous media.

## References

## References

  * Blunt et al. [2013] Martin J Blunt, Branko Bijeljic, Hu Dong, Oussama Gharbi, Stefan Iglauer, Peyman Mostaghimi, Adriana Paluszny, and Christopher Pentland.  Pore-scale imaging and modelling.  _Advances in Water resources_ , 51:197–216, 2013.
  * Bultreys et al. [2016] Tom Bultreys, Wesley De Boever, and Veerle Cnudde.  Imaging and image-based fluid transport modeling at the pore scale in geological materials: A practical introduction to the current state-of-the-art.  _Earth-Science Reviews_ , 155:93–128, 2016.
  * Golparvar et al. [2018] Amir Golparvar, Yingfang Zhou, Kejian Wu, Jingsheng Ma, and Zhixin Yu.  A comprehensive review of pore scale modeling methodologies for multiphase flow in porous media.  _Advances in Geo-Energy Research_ , 2(4):418–440, 2018.
  * Zhao et al. [2019] Benzhong Zhao, Christopher W MacMinn, Bauyrzhan K Primkulov, Yu Chen, Albert J Valocchi, Jianlin Zhao, Qinjun Kang, Kelsey Bruning, James E McClure, Cass T Miller, et al.  Comprehensive comparison of pore-scale models for multiphase flow in porous media.  _Proceedings of the National Academy of Sciences_ , 116(28):13799–13806, 2019.
  * Chen et al. [2022] Li Chen, An He, Jianlin Zhao, Qinjun Kang, Zeng-Yao Li, Jan Carmeliet, Naoki Shikazono, and Wen-Quan Tao.  Pore-scale modeling of complex transport phenomena in porous media.  _Progress in Energy and Combustion Science_ , 88:100968, 2022.  ISSN 0360-1285.
  * Gao et al. [2017] Jinfang Gao, Huilin Xing, Zhiwei Tian, Julie K Pearce, Mohamed Sedek, Suzanne D Golding, and Victor Rudolph.  Reactive transport in porous media for co2 sequestration: Pore scale modeling using the lattice boltzmann method.  _Computers & Geosciences_, 98:9–20, 2017.
  * Tahmasebi et al. [2017] Pejman Tahmasebi, Muhammad Sahimi, Amir H Kohanpur, and Albert Valocchi.  Pore-scale simulation of flow of co2 and brine in reconstructed and actual 3d rock cores.  _Journal of Petroleum Science and Engineering_ , 155:21–33, 2017.
  * Basirat et al. [2017] Farzad Basirat, Zhibing Yang, and Auli Niemi.  Pore-scale modeling of wettability effects on co2–brine displacement during geological storage.  _Advances in Water Resources_ , 109:181–195, 2017.  ISSN 0309-1708.
  * Masoudi et al. [2021] Mohammad Masoudi, Hossein Fazeli, Rohaldin Miri, and Helge Hellevang.  Pore scale modeling and evaluation of clogging behavior of salt crystal aggregates in co2-rich phase during carbon storage.  _International Journal of Greenhouse Gas Control_ , 111:103475, 2021.
  * Payton et al. [2022] Ryan L Payton, Yizhuo Sun, Domenico Chiarella, and Andrew Kingdon.  Pore scale numerical modelling of geological carbon storage through mineral trapping using true pore geometries.  _Transport in Porous Media_ , 141(3):667–693, 2022.
  * Zhao et al. [2010] Xiucai Zhao, Martin J Blunt, and Jun Yao.  Pore-scale modeling: Effects of wettability on waterflood oil recovery.  _Journal of Petroleum Science and Engineering_ , 71(3-4):169–178, 2010.
  * Kallel et al. [2017] Wissem Kallel, Marinus Izaak Jan Van Dijke, Kenneth Stuart Sorbie, and Rachel Wood.  Pore-scale modeling of wettability alteration during primary drainage.  _Water Resources Research_ , 53(3):1891–1907, 2017.
  * Zhu et al. [2017] Guangpu Zhu, Jun Yao, Aifen Li, Hai Sun, and Lei Zhang.  Pore-scale investigation of carbon dioxide-enhanced oil recovery.  _Energy & Fuels_, 31(5):5324–5332, 2017.
  * Raeini et al. [2018] Ali Q Raeini, Branko Bijeljic, and Martin J Blunt.  Generalized network modeling of capillary-dominated two-phase flow.  _Physical Review E_ , 97(2):023308, 2018.
  * Su et al. [2018] Junwei Su, Le Wang, Zhaolin Gu, Yunwei Zhang, and Chungang Chen.  Advances in pore-scale simulation of oil reservoirs.  _Energies_ , 11(5):1132, 2018.
  * Prat [2002] Marc Prat.  Recent advances in pore-scale models for drying of porous media.  _Chemical engineering journal_ , 86(1-2):153–164, 2002.
  * Surasani et al. [2008] VK Surasani, T Metzger, and E Tsotsas.  Influence of heating mode on drying behavior of capillary porous media: Pore scale modeling.  _Chemical engineering science_ , 63(21):5218–5228, 2008.
  * Panda et al. [2022] Debashis Panda, Supriya Bhaskaran, Shubhani Paliwal, Abdolreza Kharaghani, Evangelos Tsotsas, and Vikranth Kumar Surasani.  Pore-scale physics of drying porous media revealed by lattice boltzmann simulations.  _Drying Technology_ , 40(6):1114–1129, 2022.
  * Fei et al. [2022] Linlin Fei, Feifei Qin, Jianlin Zhao, Dominique Derome, and Jan Carmeliet.  Pore-scale study on convective drying of porous media.  _Langmuir_ , 38(19):6023–6035, 2022.
  * Zhao et al. [2022] Jianlin Zhao, Feifei Qin, Qinjun Kang, Dominique Derome, and Jan Carmeliet.  Pore-scale simulation of drying in porous media using a hybrid lattice boltzmann: pore network model.  _Drying Technology_ , 40(4):719–734, 2022.
  * Mukherjee et al. [2011] Partha P Mukherjee, Qinjun Kang, and Chao-Yang Wang.  Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells—progress and perspective.  _Energy & Environmental Science_, 4(2):346–369, 2011.
  * Molaeimanesh and Akbari [2014] GR Molaeimanesh and MH Akbari.  A three-dimensional pore-scale model of the cathode electrode in polymer-electrolyte membrane fuel cell by lattice boltzmann method.  _Journal of Power Sources_ , 258:89–97, 2014.
  * Zhu et al. [2021] Lijun Zhu, Heng Zhang, Liusheng Xiao, Aimy Bazylak, Xin Gao, and Pang-Chieh Sui.  Pore-scale modeling of gas diffusion layers: Effects of compression on transport properties.  _Journal of Power Sources_ , 496:229822, 2021.
  * Fu et al. [2022] Ya-lu Fu, Biao Zhang, Xun Zhu, Ding-ding Ye, Pang-Chieh Sui, Ned Djilali, and Qiang Liao.  Pore-scale modeling of mass transport in the air-breathing cathode of membraneless microfluidic fuel cells.  _International Journal of Heat and Mass Transfer_ , 188:122590, 2022.
  * Guo et al. [2022] Lingyi Guo, Li Chen, Ruiyuan Zhang, Ming Peng, and Wen-Quan Tao.  Pore-scale simulation of two-phase flow and oxygen reactive transport in gas diffusion layer of proton exchange membrane fuel cells: Effects of nonuniform wettability and porosity.  _Energy_ , 253:124101, 2022.
  * Hoogland et al. [2016a] Frouke Hoogland, Peter Lehmann, Rajmund Mokso, and Dani Or.  Drainage mechanisms in porous media: From piston-like invasion to formation of corner flow networks.  _Water Resources Research_ , 52(11):8413–8436, 2016a.
  * Hoogland et al. [2016b] Frouke Hoogland, Peter Lehmann, and Dani Or.  Drainage dynamics controlled by corner flow: Application of the foam drainage equation.  _Water Resources Research_ , 52(11):8402–8412, 2016b.
  * Moura et al. [2019] Marcel Moura, Eirik Grude Flekkøy, Knut Jørgen Måløy, Gerhard Schäfer, and Renaud Toussaint.  Connectivity enhancement due to film flow in porous media.  _Physical Review Fluids_ , 4(9):094102, 2019.
  * Armstrong and Berg [2013] Ryan T Armstrong and Steffen Berg.  Interfacial velocities and capillary pressure gradients during haines jumps.  _Physical Review E_ , 88(4):043010, 2013.
  * Moebius and Or [2014] Franziska Moebius and Dani Or.  Pore scale dynamics underlying the motion of drainage fronts in porous media.  _Water Resources Research_ , 50(11):8441–8457, 2014.
  * Jørgen Måløy et al. [2021] Knut Jørgen Måløy, Marcel Moura, Alex Hansen, Eirik Grude Flekkøy, and Renaud Toussaint.  Burst dynamics, up-scaling and dissipation of slow drainage in porous media.  _arXiv e-prints_ , pages arXiv–2110, 2021.
  * Mansouri-Boroujeni et al. [2023] Mahdi Mansouri-Boroujeni, Cyprien Soulaine, Mohamed Azaroual, and Sophie Roman.  How interfacial dynamics controls drainage pore-invasion patterns in porous media.  _Advances in Water Resources_ , 171:104353, 2023.
  * Zhou et al. [1997] Dengen Zhou, Martin Blunt, and FM Orr Jr.  Hydrocarbon drainage along corners of noncircular capillaries.  _Journal of colloid and interface science_ , 187(1):11–21, 1997.
  * Tuller and Or [2001] Markus Tuller and Dani Or.  Hydraulic conductivity of variably saturated porous media: Film and corner flow in angular pore space.  _Water resources research_ , 37(5):1257–1276, 2001.
  * Flekkøy et al. [2002] Eirik G Flekkøy, Jean Schmittbuhl, Finn Løvholt, Unni Oxaal, Knut Jørgen Måløy, and Per Aagaard.  Flow paths in wetting unsaturated flow: Experiments and simulations.  _Physical review E_ , 65(3):036312, 2002.
  * Ryazanov et al. [2009] Andrey V Ryazanov, Marinus Izaak Jan Van Dijke, and Kenneth Stuart Sorbie.  Two-phase pore-network modelling: existence of oil layers during water invasion.  _Transport in Porous Media_ , 80:79–99, 2009.
  * Han et al. [2009] M Han, S Youssef, E Rosenberg, M Fleury, and P Levitz.  Deviation from archie’s law in partially saturated porous media: Wetting film versus disconnectedness of the conducting phase.  _Physical Review E_ , 79(3):031127, 2009.
  * Romano et al. [2011] Marta Romano, Max Chabert, Amandine Cuenca, and Hugues Bodiguel.  Strong influence of geometrical heterogeneity on drainage in porous media.  _Physical Review E_ , 84(6):065302, 2011.
  * Xu et al. [2014] Wei Xu, Jeong Tae Ok, Feng Xiao, Keith B Neeves, and Xiaolong Yin.  Effect of pore geometry and interfacial tension on water-oil displacement efficiency in oil-wet microfluidic porous media analogs.  _Physics of Fluids_ , 26(9):093102, 2014.
  * Vahid Dastjerdi et al. [2022] Samaneh Vahid Dastjerdi, Nikolaos Karadimitriou, S Majid Hassanizadeh, and Holger Steeb.  Experimental evaluation of fluid connectivity in two-phase flow in porous media during drainage.  _Water Resources Research_ , 58(11):e2022WR033451, 2022.
  * Herminghaus [2005] S. Herminghaus.  Dynamics of wet granular matter.  _Advances in Physics_ , 54(3):221–261, 2005.  doi: 10.1080/00018730500167855.
  * Rieser et al. [2015] Jennifer M. Rieser, P. E. Arratia, A. G. Yodh, J. P. Gollub, and D. J. Durian.  Tunable capillary-induced attraction between vertical cylinders.  _Langmuir_ , 31(8):2421–2429, 2015.  doi: 10.1021/la5046139.  PMID: 25646573.
  * Cejas et al. [2018] Cesare M Cejas, Lawrence A Hough, Christian Frétigny, and Rémi Dreyfus.  Effect of geometry on the dewetting of granular chains by evaporation.  _Soft Matter_ , 14(34):6994–7002, 2018.
  * Chen et al. [2017] Chen Chen, Paul Duru, Pierre Joseph, Sandrine Geoffroy, and Marc Prat.  Control of evaporation by geometry in capillary structures. from confined pillar arrays in a gap radial gradient to phyllotaxy-inspired geometry.  _Scientific reports_ , 7(1):15110, 2017.
  * Chen et al. [2018] C. Chen, P. Joseph, S. Geoffroy, M. Prat, and P. Duru.  Evaporation with the formation of chains of liquid bridges.  _Journal of Fluid Mechanics_ , 837:703–728, 2018.  doi: 10.1017/jfm.2017.827.
  * Vorhauer et al. [2015] Nicole Vorhauer, YJ Wang, A Kharaghani, Evangelos Tsotsas, and Marc Prat.  Drying with formation of capillary rings in a model porous medium.  _Transport in Porous Media_ , 110:197–223, 2015.
  * Kharaghani et al. [2021] Abdolreza Kharaghani, Hafiz Tariq Mahmood, Yujing Wang, and Evangelos Tsotsas.  Three-dimensional visualization and modeling of capillary liquid rings observed during drying of dense particle packings.  _International Journal of Heat and Mass Transfer_ , 177:121505, 2021.
  * Wilkinson and Willemsen [1983] David Wilkinson and Jorge F Willemsen.  Invasion percolation: a new form of percolation theory.  _Journal of physics A: Mathematical and general_ , 16(14):3365, 1983.
  * Blunt [2001] Martin J Blunt.  Flow in porous media—pore-network models and multiphase flow.  _Current opinion in colloid & interface science_, 6(3):197–207, 2001.
  * Zhao et al. [2016] Benzhong Zhao, Christopher W MacMinn, and Ruben Juanes.  Wettability control on multiphase flow in patterned microfluidics.  _Proceedings of the National Academy of Sciences_ , 113(37):10251–10256, 2016.
  * Måløy et al. [2021] Knut Jørgen Måløy, Marcel Moura, Alex Hansen, Eirik Grude Flekkøy, and Renaud Toussaint.  Burst dynamics, upscaling and dissipation of slow drainage in porous media.  _Frontiers in Physics_ , page 718, 2021.
  * Vincent-Dospital et al. [2022] Tom Vincent-Dospital, Marcel Moura, Renaud Toussaint, and Knut Jørgen Måløy.  Stable and unstable capillary fingering in porous media with a gradient in grain size.  _Communications Physics_ , 5(1):306, 2022.
  * Méheust et al. [2002] Yves Méheust, Grunde Løvoll, Knut Jørgen Måløy, and Jean Schmittbuhl.  Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects.  _Physical Review E_ , 66(5):051603, 2002.
  * Løvoll et al. [2005] Grunde Løvoll, Yves Méheust, Knut Jørgen Måløy, Eyvind Aker, and Jean Schmittbuhl.  Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study.  _Energy_ , 30(6):861–872, 2005.
  * Moura et al. [2015] M Moura, E-A Fiorentino, Knut Jørgen Måløy, Gerhard Schäfer, and Renaud Toussaint.  Impact of sample geometry on the measurement of pressure-saturation curves: Experiments and simulations.  _Water Resources Research_ , 51(11):8900–8926, 2015.
  * Primkulov et al. [2018] Bauyrzhan K Primkulov, Stephen Talman, Keivan Khaleghi, Alireza Rangriz Shokri, Rick Chalaturnyk, Benzhong Zhao, Christopher W MacMinn, and Ruben Juanes.  Quasistatic fluid-fluid displacement in porous media: Invasion-percolation through a wetting transition.  _Physical Review Fluids_ , 3(10):104001, 2018\.
  * Primkulov et al. [2021] Bauyrzhan K Primkulov, Amir A Pahlavan, Xiaojing Fu, Benzhong Zhao, Christopher W MacMinn, and Ruben Juanes.  Wettability and lenormand’s diagram.  _Journal of Fluid Mechanics_ , 923:A34, 2021.
  * Primkulov et al. [2022] Bauyrzhan K Primkulov, Benzhong Zhao, Christopher W MacMinn, and Ruben Juanes.  Avalanches in strong imbibition.  _Communications Physics_ , 5(1):52, 2022.
  * Brakke [1992] Kenneth A Brakke.  The surface evolver.  _Experimental mathematics_ , 1(2):141–165, 1992\.
  * Joekar-Niasar and Hassanizadeh [2012] V Joekar-Niasar and SM Hassanizadeh.  Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review.  _Critical reviews in environmental science and technology_ , 42(18):1895–1976, 2012.