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

Paula Reis, Marcel Moura, Gaute Linga, Per Arne Rikvold, Renaud

Toussaint, Eirik Grude Flekkøy, Knut Jørgen Måløy

► To cite this version:

Paula Reis, Marcel Moura, Gaute Linga, Per Arne Rikvold, Renaud Toussaint, et al. A simplified pore-scale model for slow drainage including film-flow effects. Advances in Water Resources, 2023, 182, pp.104580. 10.1016/j.advwatres.2023.104580. hal-04311793

HAL Id: hal-04311793

https://hal.science/hal-04311793v1

Submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Paula Reis<sup>a,*</sup>, Marcel Moura<sup>a</sup>, Gaute Linga<sup>a</sup>, Per Arne Rikvold<sup>a</sup>, Renaud Toussaint<sup>a,d</sup>, Eirik Grude Flekkøy<sup>a,b</sup>, Knut Jørgen Måløy<sup>a,c</sup>

## ARTICLE INFO

Keywords:

Slow drainage

Film flow

Pore-scale modeling

Granular porous media

Multiphase flow

## ABSTRACT

Liquid-filled corners and capillary bridges can establish networks connecting 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 is 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 is determined by the bonds' invasion thresholds and local capillary pressure values, calculated taking into consideration gravitational effects. With the proposed model, experimentally verified phenomena related to slow drainage in granular porous media are reproduced, showing good qualitative agreement.

## 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 (Blunt et al., 2013; Bultreys et al., 2016; Golparvar et al., 2018; Zhao et al., 2019; Chen et al., 2022). In particular, extensive efforts have been directed toward the development of models related to CO₂ storage in subsurface geological reservoirs (Gao et al., 2017; Tahmasebi et al., 2017; Basirat et al., 2017; Masoudi et al., 2021; Payton et al., 2022), oil recovery methods (Zhao et al., 2010; Kallel et al., 2017; Zhu et al., 2017; Raeini et al., 2018; Su et al., 2018), drying (Prat, 2002; Surasani et al., 2008; Panda et al., 2022; Fei et al., 2022; Zhao et al., 2022), and water management in fuel cells (Mukherjee et al., 2011; Molaeimanesh and Akbari, 2014; Zhu et al., 2021; Fu et al., 2022; Guo et al., 2022). 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 pore 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 in a computationally efficient 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 (Hoogland et al., 2016a,b; Moura et al., 2019). 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 (Haines, 1930; Måløy et al., 1992; Furuberg et al., 1996; Armstrong and Berg, 2013; Moebius and Or, 2014; Måløy et al., 2021; Mansouri-Boroujeni et al., 2023). These sudden interface jumps are triggered when a pore entry pressure threshold ($P_t$) is exceeded by the capillary pressure ($P_c$), as expressed by:

$$P_c > P_t = \gamma \cos \theta \left( \frac{1}{r_1} + \frac{1}{r_2} \right) \qquad (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. Haines jumps are associated with

<footnote><sup>a</sup> PoreLab, The NJORD Centre, Department of Physics, University of Oslo, Oslo, Norway</footnote>

<footnote><sup>b</sup> PoreLab, Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway</footnote>

<footnote><sup>c</sup> PoreLab, Department of Geoscience and Petroleum, Norwegian University of Science and Technology, Trondheim, Norway</footnote>

<footnote><sup>d</sup> ITES/Institut Terre et Environnement de Strasbourg, CNRS UMR7063 - Université de Strasbourg, Strasbourg, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: paula.reis@mn.uio.no (P. Reis).</footnote>
a local transient pressure response and fluid reconfiguration (Sun and Santamarina, 2019). Piston-like displacement is generally responsible for most of the defending fluid mobilization and will be referred to as the primary drainage mechanism.

In the regions of the porous medium swept by the invading fluid front, the primary mechanism is followed by the displacement of the wetting fluid through film flow (Zhou et al., 1997; Tuller and Or, 2001). 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 (Flekkøy et al., 2002; Ryazanov et al., 2009; Han et al., 2009; Moura et al., 2019). Given the comparatively lower permeability of the corner/bridge flow paths, the secondary mechanism is characterized by distinctively longer time scales for drainage than the primary mechanism (Hoogland et al., 2016a; Moura et al., 2019).

The efficiency of both drainage mechanisms is inherently tied to geometric aspects of the pore space (Romano et al., 2011; Xu et al., 2014; Vahid Dastjerdi et al., 2022). During piston-like displacement, the size distribution and spatial arrangement of pores 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 porous 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 (Herminghaus, 2005; Rieser et al., 2015; Cejas et al., 2018; Chen et al., 2017, 2018; Moura et al., 2019) 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. (2018) 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 above which the liquid connections break and water removal from the porous medium is undermined. In a similar set of studies, Chen et al. (2017) and Chen et al. (2018) investigated the effect of geometry in the formation of elongated liquid films using transparent quasi-2D micromodels of cylinders arranged between two planes. Depending on the spatial distribution 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. (2019) to evaluate the impact of film flow on slow drainage under the influence of varying gravitational forces. 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. (2015), 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. (2021), 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 monodispers 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. (2018), 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 pore space.

In the present work, we incorporate the capillary-bridge role in wetting-phase connectivity successfully explored in drying models into 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. (2019), in addition to some new analyses of their experimental data, 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 (Wilkinson and Willemse, 1983). Designed to represent immiscible fluid-fluid displacements in porous media, this method involves the description of pores and throats as networks of sites and bonds, to each of which an invasion threshold is assigned. Flow in the idealized porous networks is then predicted based on the movement of the fluid interface through the least resistivity path, dictated by the available bonds' invasion thresholds in the case of drainage, and the available sites' invasion thresholds in the case of imbibition (Wilkinson and Willemse, 1983). Only bonds and sites located at the interface of the main defending-fluid cluster are available for invasion, and disconnected finite clusters are considered trapped. In particular, a loose trapping assumption is adopted in our model, meaning that fluid can be drained through bonds connected to the same site, irrespective of the site's fluid occupation (Joekar-Niasar et al., 2008). This modeling approach has been widely adopted in the past decades (Blunt, 2001) 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 is specially tailored to reproduce slow drainage – i.e., drainage in the absence of significant viscous and inertial forces – in Hele-Shaw cells filled with a single layer of monodispersed spherical beads. This choice of quasi-2D model porous media allows us to clearly identify drainage events related to film flow and generate results that could be readily compared with available experimental data (Moura et al., 2019).

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
<figure><img src="image_4.png" /><figcaption>Fig. 1. Examples of wetting-fluid configurations around pairs of cylinders and spheres, represented as triangular meshes, contained between transparent parallel planes, obtained with Surface Evolver (Brakke, 1992). 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.</figcaption></figure>

geometric feature of quasi-2D models is the intersection of the structures representing the porous media grains, e.g. cylinders (Zhao et al., 2016; Chen et al., 2017; Måløy et al., 2021; Vincent-Dospital et al., 2022) or spheres (Méheust et al., 2002; Løvoll et al., 2005; Moura et al., 2015, 2019), 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 Fig. 1. Several recent pore-scale models incorporate aspects of these fluid configurations stemming from 3D structures in 2D lattices representing the porous space (Chen et al., 2018; Primkulov et al., 2018, 2021, 2022), and we here adopt a similar approach.

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 (Blunt, 2001; Joekar-Niasar and Hassanizadeh, 2012). 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. (2019). Their porous-media models were created by randomly placing and securing spherical glass beads with a diameter of 1 mm 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 are approximated as regular honeycomb lattices, as indicated in Fig. 2.

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 do not address the geometrical aspects of the lattice nodes and assume 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 an object composed of a set of nodes and a set of edges. Nodes can be, but are not necessarily, connected by edges. 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 a 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.

The representation of the role of capillary and gravitational forces in the drainage process is achieved by assigning weights to the graph's edges. In analogy with Eq. (1), these edge weights are calculated as the difference between the local capillary pressure (Eq. (2a)), and the local capillary pressure threshold for invasion of each throat, given by the Young-Laplace Eq. (2b). In this calculation, the local capillary pressure, $P_c(y)$, varies as a function of the vertical distance between a throat and the porous-medium outlet, to which a reference pressure $P_0$ is assigned. In this way, we map the hydrostatic pressure field into the pressure difference between the wetting and non-wetting phases, incorporating the effect of gravity during slow drainage (Wilkinson, 1984; Birovljev et al., 1991; Auradou et al., 1999; Ayaz et al., 2020).

$$P_c(y) = P_0 + \Delta\rho gy \sin \beta \quad (2a)$$

$$P_t = 2\gamma \cos \theta \left( \frac{1}{t_w} + \frac{1}{t_h} \right) \quad (2b)$$

where $P_0$ is the wetting-phase pressure at the outlet, $\Delta\rho$ is the density difference between wetting and non-wetting fluids, $g$ is the gravitational acceleration, $y$ is the perpendicular distance between a throat and the model's outlet, $\beta$ is the inclination angle between model and the horizontal plane, and $t_w$ and $t_h$ are the throat's width and height, respectively.

As with IP models, our drainage simulation unfolds in discrete steps, in each of which the available bond with the lowest resistance is invaded. This invasion is expressed by deleting the available edge with the largest weight from the graph representing the wetting-fluid occupation of the porous medium, as indicated in Fig. 3. At the beginning of the drainage simulation, only the edges incident on the inlet nodes are considered available for invasion (Fig. 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 (Fig. 3b-d). The reduction 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 (Fig. 3e). Following these rules, the drainage process continues in the model until all available edges are deleted.

### 2.2.1. Cluster merging via capillary bridges

The drainage process illustrated in Fig. 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 connects otherwise disconnected clusters, as suggested in Fig. 4.

In the upper left corner of Fig. 4, a picture of a Hele-Shaw cell under drainage from Moura et al. (2019) is shown. In this insert, the two arrows highlight chains formed by capillary bridges through which fluid transport is permissible. In a standard 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
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Sketch of the modified Hele-Shaw porous matrix used in Moura et al. (2019) and its representation by a regular honeycomb lattice of nodes (open circles) and edges (solid lines).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Representation of the modified IP model algorithm for drainage.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Sketch of Hele-Shaw porous-space representation into a 2D regular honeycomb lattice of nodes and edges.</figcaption></figure>

in the sketch on the right side of Fig. 4: capillary bridges appear in the graph as edges incident on two nodes filled with the non-wetting phase. Notice that the connection they provide for the wetting fluid, indicated by the arrow, matches the dual lattice of the network representing the porous space, as suggested in the experimental work of Moura et al. (2019).

The capillary-bridge 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
<figure><img src="image_8.png" /><figcaption>Fig. 5. Edge configurations in a lattice unit cell that can lead to cluster-merging.</figcaption></figure>

around the same spheres, as illustrated in Fig. 1. To further clarify the implemented connectivity rule, Fig. 5 displays all the possible edge configurations in a lattice unit cell that can lead to the merging of otherwise disconnected components. Note that configurations A, D, and E were exemplified in the sketch contained in Fig. 4.

As an additional modification to the standard IP model, we allow the snap-off of capillary bridges connected to the main defending cluster, conceding that the local capillary pressure exceeds the throat's $P_t$. Although the snap-off pressure of a bridge in the quasi-2D porous media under investigation should not be identical to the threshold for throat invasion, we adopt this approximation as a model-simplifying assumption. More detailed analyses regarding the hydrodynamic instability of capillary bridges held between spherical beads and parallel planes, as depicted in Fig. 1, will be incorporated in a future study.

With the proposed modifications related to the cluster-labeling algorithm and the possibility of capillary bridge snap-off, our model could reproduce experimentally verified aspects of drainage events related to both primary and secondary mechanisms, as presented in Section 3. In the next section, the properties of the porous media and fluids represented in our numerical study are described.

## 2.3. Porous-medium and fluid properties

The results presented in this work 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 15 × 26 cm. The capillary pressure threshold values of the edges representing throats follow the distribution presented in Fig. 6, based on the experiments presented by Moura et al. (2019). In that study, glass spheres with a diameter of 1 mm were placed between the Hele-Shaw cell plates, and the reported separation between the spheres varied approximately from 0 to 1.5 mm, characterizing a porous medium with properties compatible with coarse-grained soils (Zhai et al., 2020). Adopted fluid properties also proceed from the same experimental work (Moura et al., 2019), and correspond to an interfacial tension between wetting and non-wetting fluids of $\gamma = 64$ mN/m, a density contrast between the fluids of $\Delta\rho = 1204$ kg/m$^3$, and an assumption that the wetting-fluid perfectly wets the porous medium. Such properties lead to a capillary length $\lambda_c = \sqrt{\gamma/\Delta\rho g} \approx 2$ mm, much smaller than the analyzed porous-media sample size, meaning that gravitational forces are likely to impact the investigated flow. Therefore, the influence of gravity on 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}$.

In the lattices, lateral boundaries ($x = 0$ cm and $x = 15$ cm) were closed to the flow, nodes at the upper boundary ($y = 26$ cm) were sources of non-wetting fluid, and all nodes at the lower boundary ($y = 0$ cm) 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><img src="image_9.png" /><figcaption>Fig. 6. Probability density function of the capillary pressure threshold ($P_t$) of the porous matrix throats, generated to match the values reported in Moura et al. (2019).</figcaption></figure>

## 3. Results

Several experimentally reported (Moura et al., 2019) aspects related to drainage through capillary bridges and liquid rings could be reproduced with the proposed numerical model, and are described in this section. An analysis of the advantages and shortcomings of our modified-IP approach to model slow drainage in granular porous media is presented thereafter, in Section 4.

Fig. 7 illustrates the change in drainage dynamics predicted by the proposed model as the influence of the gravitational component increases. These results were generated using the same network and with inclination angles indicated in the captions. From 7(a) to 7(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, therefore, correspond to the trapped wetting fluid clusters. We can observe that as $\beta$ increases, the unstable front characteristic of capillary fingering is progressively stabilized, resulting in a more efficient removal of the wetting phase. This trend has been observed both in experimental (Méheust et al., 2002; Løvoll et al., 2005) and numerical (Birovljev et al., 1991) studies of slow drainage in porous media. The invasion dynamics presented here 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
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Representation of the invasion order in the porous matrix during the drainage simulations, from dark blue (early invasion) to dark red (late invasion), at the moment of the non-wetting phase breakthrough. In the images, the lattice's bottom boundary, where the outlet is located, corresponds to $y = 0$ cm, while the top boundary, containing the inlet, corresponds to $y = 26$ cm.</figcaption></figure>

for different degrees of gravitational forces influence is illustrated in Fig. 8.

Complementary to the information presented in Fig. 7, the images in Fig. 8 introduce the distinction between drained pores and throats that pertained directly to the main defending cluster (teal), 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 Fig. 9(a), while corresponding experimental data from Moura et al. (2019) is presented in 9(b).

As indicated in 9(a), 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 relative importance of film flow to drainage, the maximum mean value is achieved at $\beta = 2.5^\circ$. This behavior may be explained due to the opposing 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 from the front, the increase in capillary pressure due to the hydrostatic pressure may place them in a preferential position for invasion, when compared to nodes belonging to the front, as described in Eq. (2b). 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.

The results shown in Fig. 9(a) exhibit a reasonable agreement with experimental data obtained by Moura et al. (2019), shown in 9(b). The amount of fluid recovered with the primary drainage mechanism measured with both methods is almost equivalent, while the drainage related to the secondary mechanism seems to be overestimated in the numerical simulations. At low values of $\beta$, however, a one-to-one comparison between the two sets of results is not straightforward, as the simulated results with different random networks displayed considerable variations and the experiments were not repeated. Still, it is noticeable that in both cases the flow through capillary bridges is minimized at $\beta = 0^\circ$, and its contribution to the total drainage varies non-monotonically with the influence of gravitational forces. Moreover, this discrepancy between experimental and simulated results is expected due to the model's idealized assumptions, as discussed in Section 4, and we intend to address these shortcomings in future work.

### 3.2. The film flow active zone

In this section, we analyze the formation of a region trailing the invasion front where drainage events related to film flow are more likely
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Representation of the drainage mechanisms related to invasion in the porous matrix during the drainage simulations. Nodes and edges in teal were drained via the primary drainage mechanism, while nodes and edges in yellow were drained by the secondary mechanism. In the images, the lattice's bottom boundary, where the outlet is located, corresponds to y = 0 cm, while the top boundary, containing the inlet, corresponds to y = 26 cm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Wetting-fluid saturation drained via the primary (S₁, blue curve) and secondary (S₂, pink curve) mechanisms. The green curve represents the relative contribution of the secondary mechanisms to drainage, based on how much wetting fluid would be drained via the primary mechanism alone. Results presented in (a) correspond to this numerical work, and its error bars represent the standard deviation in the measurements obtained with a hundred random realizations of the lattice representing the porous space. Results presented in (b) correspond to the experimental work conducted by Moura et al. (2019), and are obtained with a single experiment for each inclination.</figcaption></figure>

to occur. First identified experimentally by Moura et al. (2019), 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, capillary bridges must form a continuous path
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** 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.</figcaption></figure>

between the invasion front and the wetting-phase clusters behind it. This imposes a constraint on how wide the active zone can be, as long sequences of connected capillary bridges are significantly less likely to occur than single bridges or short bridge chains. 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.

Fig. 10 illustrates the nodes connected to the main defending cluster via corners and capillary bridges, halfway through drainage, under different degrees of influence of gravitational forces. We will use the size of the region covered by these nodes as a proxy for the extent of the active zone, as it encloses the domain in which the secondary drainage can take place. From Figs. 10(a) to 10(h), there is a clear indication that the active zone gets progressively more compact as $β$ increases. For the cases where the gravitational component is zero or relatively small, $0° ≤ β ≤ 2.5°$, the frequency of large wetting fluid clusters behind the front – as noticeable in Figs. 7 and 8 – lead to very large parts of the porous matrix being 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° ≤ β < 10°$, followed by a long inclination interval, $10° ≤ β ≤ 90°$, where a narrow active zone is curtailed slowly. This compaction of the active zone takes place simultaneously with the reduction of wetting-phase-cluster sizes behind the front, and can also be attributed to the increased chance of capillary bridge snap-off at steeper angles due to the increased local capillary pressure values.

The probability density functions (PDFs) of the distance from the front to the nodes connected and drained via film flow are shown in Fig. 11. As suggested before, the probability of finding nodes pertaining to the active zone closer to the front increases significantly as gravitational forces become more pronounced during drainage. The same trend was observed experimentally in Moura et al. (2019) when comparing the PDF of the distance from the invasion front of film flow events in experiments with $β = 10°, 20°$ and $30°$. Interestingly, the normalized PDFs presented in 11(a) and 11(b) are remarkably similar. This resemblance 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. Besides, it supports the validity of using the dynamic networks' sizes, as shown in Fig. 10, as a representative of the active zone width.

Another interesting fact is that wider active zones are not necessarily linked to more significant effects of film flow events in drainage. Fig. 12 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 Figs. 12 and 9(a) 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. This finding substantiates the claim presented in Section 3.1 that gravitational effects exert a significant influence on the occurrence of the secondary drainage mechanism.

The characterization of the active zone using the size of the secondary-drainage-mechanism networks should also take into account the fact that the connectivity of such networks is under constant
<figure><img src="image_14.png" /><figcaption>Fig. 11. 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).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. 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 error bars represent the standard deviation in the measurements obtained with a hundred random realizations of the lattice representing the porous space.</figcaption></figure>

evolution as the porous matrix is drained. This phenomenon was acknowledged in Moura et al. (2019), where the networks associated with primary drainage were described 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 magnitude of the fluctuation in the dynamic networks’ sizes during drainage can be inferred from the error bars in Fig. 12. Additionally, an insight into these variations is presented in Fig. 13. This figure illustrates the secondary drainage networks at eight equidistant stages of a drainage simulation with $\beta = 7.5°$. We can observe that even though their width is somewhat constant, their shape and size can vary significantly. A reason pointed out in Moura et al. (2019) 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 Fig. 14, which represents the region around the invasion front in Fig. 13(d), the five capillary bridges indicated by 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.

### 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. Fig. 15(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°$ to $\beta = 20°$ there is a sharp increase in the occurrence of bridges, which then stabilizes at approximately 18.5%. The comparison of these results with the ones presented in Figs. 9(a) and 12 indicates that the pervasiveness of capillary bridges alone is not correlated with a greater capacity of connecting, or draining, extensive regions in the lattice to the main defending cluster. In fact, the increase in the 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 relative frequency of capillary bridges as inclination angles get steeper was also observed experimentally by Moura et al. (2019) when comparing cases with $\beta = 10°$ and $20°$. In that study, a fraction of 12.8% of throats was occupied by bridges at $\beta = 10°$, and 15% at $\beta = 20°$, which is approximately 20% lower than the values predicted by the model.

Complementary to the information presented in Figs. 15(a) and 15(b) presents the PDFs of the sizes of throats ($t_w$) in which capillary bridges were formed during drainage at different inclination angles. In
<figure><img src="image_16.png" /><figcaption>Fig. 13. Representation of nodes connected to the front via capillary bridges and corners for $\beta = 7.5^\circ$, at six stages during drainage.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Detail of the secondary drainage network shown in 13(d). The arrows point to the bridges responsible for the connection of the surrounded regions with matching colors. The depiction of the spherical beads surrounded by the networks is included for clarity.</figcaption></figure>

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 for the bridges to concentrate at the lower end of this range.

This heterogeneous occupation of throats by bridges was also identified experimentally in Moura et al. (2019), as shown in Fig. 16

(where the number of counts of the experimental results was re-scaled 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 a throat hosting a bridge was approximately 0.5 mm, 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 —
<figure><img src="image_18.png" /><figcaption>Fig. 15. Frequency of throats containing capillary bridges (a) and their size distributions (b).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Probability density function of the sizes of throats containing capillary bridges experimentally and with the proposed model, for $\beta = 20^{\circ}$.</figcaption></figure>

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, not the Plateau–Rayleigh instability (Rayleigh, 1878), is the predominant factor associated with the size of throats expected to host bridges.

## 4. Discussion

The results presented in Section 3 indicate that the proposed modified-IP model, despite its many idealizing assumptions, is able to reproduce recently experimentally identified (Moura et al., 2019) phenomena taking place during slow drainage in granular porous media. Intended to investigate the partition of drainage into two main mechanisms, namely a primary one — through the bulk of pore bodies and throats — and a secondary one — through connected corners and capillary bridges — the developed quasi-static network model was used to investigate drainage flows under the influence of capillary and gravitational forces. A summary of the findings is outlined next.

First, the influence of primary and secondary mechanisms was assessed by quantifying the fractions of the wetting fluid drained by each. As described in Section 3.1, the relative contribution of flow through corners and capillary bridges is maximized at low inclination angles of the porous media with respect to the horizontal plane, $1.25^{\circ} \le \beta \le 10^{\circ}$. In the absence of gravitational effects, i.e. $\beta = 0^{\circ}$, very little fluid is drained by the secondary mechanism, as connected clusters behind the front do not experience increased capillary forces, for being at the same elevation as the front itself. At higher inclination angles, $\beta > 10^{\circ}$, on the other hand, the contribution of the secondary mechanism is greatly

surpassed by the primary drainage mechanism efficacy, leading to a diminutive relative importance.

Once the relevance of each mechanism was quantified, we analyzed the fluids’ configuration in the pore space originated by drainage through the bulk of pores, and its link to the flow through connected capillary bridges and corners. In Section 3.2, the concept of an active zone where the secondary mechanism takes place was introduced. Analogously to experimental observations presented in Moura et al. (2019), we identified in our simulations a finite region trailing the front where clusters could be drained by the secondary mechanism. It was observed that the width of the active zone decreased monotonically with the influence of gravity and could not, therefore, be correlated to the amount of fluid drained through bridges and corners. This mismatch between active zone size and secondary drainage mechanism efficacy is attributed to the combination of large clusters formed behind the front at low $\beta$ values, and their difficulty of being drained subjected to weak gravitational forces.

The frequency and size of capillary bridges in the pore throats were also investigated, in Section 3.3. The results indicated that the percentage of throats hosting bridges undergo a threefold increase from $\beta = 0^{\circ}$ to $\beta = 45^{\circ}$, followed by a slow decrease at steeper angles. This increase can be explained by the gain in efficiency of the primary drainage with $\beta$ (see Fig. 9(a)), which enlarges the air-invaded portions of the porous matrix, where bridges can be formed. For $\beta > 45^{\circ}$, however, this effect may be overshadowed by the increased likelihood of invasion and snap-off of capillary bridges, as suggested by the decrease in the mean size of throats hosting a bridge (see Fig. 15(b)).

The proposed numeric model and drainage mechanisms analysis summed up above are deemed pertinent amidst a series of studies published in the past few years acknowledging the role of capillary bridges and liquid rings in quasi-static fluid displacements in porous media. Zhao et al. (2016) highlighted the importance of chains of
pendular rings in fluid transport during strong imbibition in patterned microfluidics. The wetting-fluid connectivity provided by these chains in imbibition flows was later incorporated in 2D pore-network models proposed by Primkulov et al. (2018, 2021) and Primkulov et al. (2022), departing from the typical concept of representing film flow in porous media as flow at corners of angular capillaries. The significance of fluid transport through capillary bridges and liquid rings during drying in porous media was also pointed out in a series of experimental and numerical (Vorhauer et al., 2015; Cejas et al., 2018; Chen et al., 2017, 2018; Kharaghani et al., 2021) studies. In this context, chains of bridges and rings provide continuous channels for liquid transport from saturated zones within the porous medium to the evaporation front, allowing the maintenance of fast drying rates. In the same direction as the studies regarding imbibition and drying, experimental observations presented in Moura et al. (2019) made a case for the relevance of capillary bridges and liquid rings in slow drainage in porous media. To the best of our knowledge, this work is the first to explicitly incorporate these extra paths for connectivity in a simple network model for drainage in porous media subjected to capillary and gravitational forces.

Beyond the relevance of its targeted flow phenomena, we consider that the significance of our model also lies in its approach to representing film flow, network conceptualization, and implementation method. Similarly to Primkulov et al. (2018, 2021) and Primkulov et al. (2022), we incorporated film flow effects in a PNM without portraying pore throats as angular capillaries with wetting fluid flowing at corners. By assuming that all solid-solid contacts in the porous medium (plane-sphere contacts in the case of the modified Hele-Shaw cell) adjacent to liquid-filled throats must also be filled with liquid, we could identify connected and disconnected fluid segments at a low-level description of the porous medium geometric features. With this, we avoided commonly adopted oversimplifying assumptions, such as that the wetting phase is connected everywhere in the porous media through corners (Valvatne and Blunt, 2004), or that all corners in a pore element act as a bundle of channels with shared start and end points (Joekar-Niasar and Hassanizadeh, 2012), while also bypassing the great level of detail required to treat every corner in pore elements individually (Raeini et al., 2018).

As for our PNM conceptualization, instead of regarding networks as a representation of the porous space, with nodes and edges either occupied by the wetting or non-wetting fluid, we took networks as surrogates for the wetting-fluid occupation of the pore space itself. To do so, edges invaded by air are deleted from the network, making the size of the investigated system smaller at each simulation time step. This approach is suitable for transient flow – as in drainage or imbibition – when the fluid occupation of a pore is not expected to change more than once. The resultant evolving networks of liquid occupation could then be analyzed using some basic graph-theory concepts, e.g. degree of nodes, paths, and connected graph components, with related algorithms readily available in various well-established numerical libraries. Besides facilitating the model implementation, this particular PNM concept allows one to easily label fluid clusters and investigate their properties regardless of the network structure.

The straightforwardness of our modeling approach begets, however, limitations to the generalizability of the obtained results. Notably, the absence of viscous effects, the lack of a precise stability criterion for the capillary bridges, and the quasi-2D nature of the analyzed porous medium can interfere with our findings, as examined next.

For starters, despite the association of slow drainage 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 secondary drainage mechanism employs much less permeable channels than the primary one, and there should be, therefore, a preference for fluid to be mobilized from the main defending cluster rather than from isolated clusters. The oversight of this effect may justify the overestimation of

drainage via the secondary mechanism observed in Fig. 9, and could be addressed in a future study with a dynamic PNM. The overestimation of the contribution of the secondary drainage mechanism could also stem from the lack of a rigorous criterion for bridge snap-off, which should take into account the link of fluid arrangements to local capillary pressure and pore geometry. By equating the capillary pressure for snap-off to the throats' threshold for invasion, we are likely overstating the bridges' stability, which is corroborated by Fig. 15(a). For this reason, a criterion derived from a meticulous investigation of stable interfaces in constricted geometries, analogous to the one presented by Cox et al. (2023) regarding microfluidic devices, could enhance the quality of our results.

Finally, the very nature of our 2D PNM – based on experiments in quasi-2D porous media – imposes limitations on the applicability of the presented findings to flow in real 3D porous media. Even though capillary bridges and liquid rings are ubiquitous in unsaturated flows in 3D porous media (Scheel et al., 2008), identifying the conditions for the extra connectivity provided by them is considerably more difficult than in 2D geometries (Cejas et al., 2018). The description of such connections was presented by Kharaghani et al. (2021) for a body-center-cubic pack of monodispersed spheres, but expanding this concept to more complex lattices or irregular packings could be challenging. Moreover, fluid structures beyond rings and bridges could arise in drainage in 3D porous media, as recently demonstrated by Wu et al. (2023). In their work, liquid filaments were observed connecting liquid rings in throats invaded by the non-wetting phase at moderate capillary numbers ($Ca > 3 \times 10^{-4}$). Although these structures are expected to display very low permeability, they provide a channel for wetting-phase connectivity not taken into consideration in our model and may expose a way in which our results fail to represent real 3D flows. Nevertheless, we consider that the findings presented in Section 3 can be, at least qualitatively, generalized to more complex porous media geometries, and that simple idealized models, as described in Section 2, are instrumental in the understanding of complex physical phenomena.

## 5. Conclusions

In this study, a quasi-static pore-network model for slow drainage in granular porous media was presented. The model was based on the idealized bond-invasion-percolation method and incorporated a modification in the trapped-cluster identification algorithm, in order to reproduce the wetting-phase connectivity provided by chains of capillary bridges. The pore space was represented in the model as a regular honeycomb lattice 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 qualitatively 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 presented in Moura et al. (2019), or newly analyzed by the authors for this work. Despite the model's simple representation of the pore 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 be a valuable tool to interpret and complement experiments, in the quest to complete our understanding of drainage in granular porous media.

As a future development for the presented model for slow drainage, we suggest the implementation of a criterion for the capillary-bridge snap-off taking into consideration the shape of the bridges formed in the quasi-2D porous medium under investigation, as illustrated in
Fig. 1. Another relevant advancement in our current line of study could be achieved with the incorporation of viscous effects in the fluid displacement represented by the model. This modification could improve the quality of our findings related to the slow drainage regime, as well as allow the investigation of the relevance of primary and secondary drainage mechanisms under varying capillary numbers.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

We acknowledge the support of the University of Oslo, the Njord Center, and the Research Council of Norway through the PoreLab Center of Excellence (project number 262644), the FlowConn Researcher Project for Young Talent (project number 324555) and the M4: Mixing in Multiphase flows through Microporous Media Researcher Project for Young Talent (project number 325819).

## References

* Armstrong, Ryan T., Berg, Steffen, 2013. Interfacial velocities and capillary pressure gradients during haines jumps. Phys. Rev. E 88 (4), 043010.

* Auradou, Harold, Måløy, Knut Jørgen, Schmittbuhl, Jean, Hansen, Alex, Bideau, Daniel, 1999. Competition between correlated buoyancy and uncorrelated capillary effects during drainage. Phys. Rev. E 60 (6), 7224.

* Ayaz, Monem, Toussaint, Renaud, Schäfer, Gerhard, Måløy, Knut Jørgen, 2020. Gravitational and finite-size effects on pressure saturation curves during drainage. Water Resour. Res. 56 (10), e2019WR026279.

* Basirat, Farzad, Yang, Zhibing, Niemi, Auli, 2017. Pore-scale modeling of wettability effects on CO₂-brine displacement during geological storage. Adv. Water Resour. 109, 181–195.

* Birovljev, A., Furuberg, L., Feder, J., Jøssang, T., Måløy, K.J., Aharony, A., 1991. Gravity invasion percolation in two dimensions: Experiment and simulation. Phys. Rev. Lett. 67 (5), 584.

* Blunt, Martin J., 2001. Flow in porous media—pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 6 (3), 197–207.

* Blunt, Martin J., Bijeljic, Branko, Dong, Hu, Gharbi, Oussama, Iglauer, Stefan, Mostaghimi, Peyman, Paluszny, Adriana, Pentland, Christopher, 2013. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216.

* Brakke, Kenneth A., 1992. The surface evolver. Exp. Math. 1 (2), 141–165.

* Bultreys, Tom, De Boever, Wesley, Cnudde, Veerle, 2016. 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-Sci. Rev. 155, 93–128.

* Cejas, Cesare M., Hough, Lawrence A, Frétigny, Christian, Dreyfus, Rémi, 2018. Effect of geometry on the dewetting of granular chains by evaporation. Soft Matter 14 (34), 6994–7002.

* Chen, Chen, Duru, Paul, Joseph, Pierre, Geoffroy, Sandrine, Prat, Marc, 2017. Control of evaporation by geometry in capillary structures. From confined pillar arrays in a gap radial gradient to phyllotaxy-inspired geometry. Sci. Rep. 7 (1), 15110.

* Chen, Li, He, An, Zhao, Jianlin, Kang, Qinjun, Li, Zeng-Yao, Carmeliet, Jan, Shikazono, Naoki, Tao, Wen-Quan, 2022. Pore-scale modeling of complex transport phenomena in porous media. Prog. Energy Combust. Sci. 88, 100968.

* Chen, C., Joseph, P., Geoffroy, S., Prat, M., Duru, P., 2018. Evaporation with the formation of chains of liquid bridges. J. Fluid Mech. 837, 703–728. http://dx.doi.org/10.1017/jfm.2017.827.

* Cox, S.J., Davarpanah, A., Rossen, W.R., 2023. Interface shapes in microfluidic porous media: Conditions allowing steady, simultaneous two-phase flow. Transp. Porous Media 147 (1), 197–216.

* Fei, Linlin, Qin, Feifei, Zhao, Jianlin, Derome, Dominique, Carmeliet, Jan, 2022. Pore-scale study on convective drying of porous media. Langmuir 38 (19), 6023–6035.

* Flekkøy, Eirik G, Schmittbuhl, Jean, Løvholt, Finn, Oxaal, Unni, Måløy, Knut Jørgen, Aagaard, Per, 2002. Flow paths in wetting unsaturated flow: Experiments and simulations. Phys. Rev. E 65 (3), 036312.

* Fu, Ya-lu, Zhang, Biao, Zhu, Xun, Ye, Ding-ding, Sui, Pang-Chieh, Djilali, Ned, Liao, Qiang, 2022. Pore-scale modeling of mass transport in the air-breathing cathode of membraneless microfluidic fuel cells. Int. J. Heat Mass Transfer 188, 122590.

* Furuberg, Liv, Måløy, Knut Jørgen, Feder, Jens, 1996. Intermittent behavior in slow drainage. Phys. Rev. E 53 (1), 966.

* Gao, Jinfang, Xing, Huilin, Tian, Zhiwei, Pearce, Julie K., Sedek, Mohamed, Golding, Suzanne D., Rudolph, Victor, 2017. Reactive transport in porous media for CO₂ sequestration: Pore scale modeling using the lattice Boltzmann method. Comput. Geosci. 98, 9–20.

* Golparvar, Amir, Zhou, Yingfang, Wu, Kejian, Ma, Jingsheng, Yu, Zhixin, 2018. A comprehensive review of pore scale modeling methodologies for multiphase flow in porous media. Adv. Geo-Energy Res. 2 (4), 418–440.

* Guo, Lingyi, Chen, Li, Zhang, Ruiyuan, Peng, Ming, Tao, Wen-Quan, 2022. 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.

* Haines, William B., 1930. Studies in the physical properties of soil. V. The hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith. J. Agric. Sci. 20 (1), 97–116.

* Han, M., Youssef, S., Rosenberg, E., Fleury, M., Levitz, P., 2009. Deviation from Archie's law in partially saturated porous media: Wetting film versus disconnectedness of the conducting phase. Phys. Rev. E 79 (3), 031127.

* Herminghaus, S., 2005. Dynamics of wet granular matter. Adv. Phys. 54 (3), 221–261. http://dx.doi.org/10.1080/00018730500167855.

* Hoogland, Frouke, Lehmann, Peter, Mokso, Rajmund, Or, Dani, 2016a. Drainage mechanisms in porous media: From piston-like invasion to formation of corner flow networks. Water Resour. Res. 52 (11), 8413–8436.

* Hoogland, Frouke, Lehmann, Peter, Or, Dani, 2016b. Drainage dynamics controlled by corner flow: Application of the foam drainage equation. Water Resour. Res. 52 (11), 8402–8412.

* Joekar-Niasar, V., Hassanizadeh, S.M., 2012. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Crit. Rev. Environ. Sci. Technol. 42 (18), 1895–1976.

* Joekar-Niasar, V., Hassanizadeh, S.M., Leijnse, A., 2008. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transp. Porous Media 74, 201–219.

* Kallel, Wissem, Van Dijke, Marinus Izaak Jan, Sorbie, Kenneth Stuart, Wood, Rachel, 2017. Pore-scale modeling of wettability alteration during primary drainage. Water Resour. Res. 53 (3), 1891–1907.

* Kharaghani, Abdolreza, Mahmood, Hafiz Tariq, Wang, Yujing, Tsotsas, Evangelos, 2021. Three-dimensional visualization and modeling of capillary liquid rings observed during drying of dense particle packings. Int. J. Heat Mass Transfer 177, 121505.

* Løvoll, Grunde, Méheust, Yves, Måløy, Knut Jørgen, Aker, Eyvind, Schmittbuhl, Jean, 2005. Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study. Energy 30 (6), 861–872.

* Måløy, Knut Jørgen, Furuberg, Liv, Feder, Jens, Jøssang, Torstein, 1992. Dynamics of slow drainage in porous media. Phys. Rev. Lett. 68 (14), 2161.

* Måløy, Knut Jørgen, Moura, Marcel, Alex, Flekkøy, Eirik Grude, Toussaint, Renaud, 2021. Burst dynamics, upscaling and dissipation of slow drainage in porous media. Front. Phys. 718.

* Mansouri-Boroujeni, Mahdi, Soulaine, Cyprien, Azaroual, Mohamed, Roman, Sophie, 2023. How interfacial dynamics controls drainage pore-invasion patterns in porous media. Adv. Water Resour. 171, 104353.

* Masoudi, Mohammad, Fazeli, Hossein, Miri, Rohaldin, Hellevang, Helge, 2021. Pore scale modeling and evaluation of clogging behavior of salt crystal aggregates in CO₂-rich phase during carbon storage. Int. J. Greenh. Gas Control 111, 103475.

* Méheust, Yves, Løvoll, Grunde, Måløy, Knut Jørgen, Schmittbuhl, Jean, 2002. Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects. Phys. Rev. E 66 (5), 051603.

* Moebius, Franziska, Or, Dani, 2014. Pore scale dynamics underlying the motion of drainage fronts in porous media. Water Resour. Res. 50 (11), 8441–8457.

* Molaeimanesh, G.R., Akbari, M.H., 2014. A three-dimensional pore-scale model of the cathode electrode in polymer-electrolyte membrane fuel cell by lattice Boltzmann method. J. Power Sources 258, 89–97.

* Moura, M., Fiorentino, E.-A., Måløy, Knut Jørgen, Schäfer, Gerhard, Toussaint, Renaud, 2015. Impact of sample geometry on the measurement of pressure-saturation curves: Experiments and simulations. Water Resour. Res. 51 (11), 8900–8926.

* Moura, Marcel, Flekkøy, Eirik Grude, Måløy, Knut Jørgen, Schäfer, Gerhard, Toussaint, Renaud, 2019. Connectivity enhancement due to film flow in porous media. Phys. Rev. Fluids 4 (9), 094102.

* Mukherjee, Partha P., Kang, Qinjun, Wang, Chao-Yang, 2011. Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells—progress and perspective. Energy Environ. Sci. 4 (2), 346–369.

* Panda, Debashis, Bhaskaran, Supriya, Paliwal, Shubhani, Kharaghani, Abdolreza, Tsotsas, Evangelos, Surasani, Vikrant Kumar, 2022. Pore-scale physics of drying porous media revealed by lattice Boltzmann simulations. Drying Technol. 40 (6), 1114–1129.

* Payton, Ryan L., Sun, Yizhuan, Chiarella, Domenico, Kingdon, Andrew, 2022. Pore scale numerical modelling of geological carbon storage through mineral trapping using true pore geometries. Transp. Porous Media 141 (3), 667–693.
Prat, Marc, 2002. Recent advances in pore-scale models for drying of porous media. Chem. Eng. J. 86 (1-2), 153-164.

Primkulov, Bauyrzhan K., Pahlavan, Amir A., Fu, Xiaojing, Zhao, Benzhong, MacMinn, Christopher W., Juanes, Ruben, 2021. Wettability and Lenormand's diagram. J. Fluid Mech. 923, A34.

Primkulov, Bauyrzhan K., Talman, Stephen, Khaleghi, Keivan, Shokri, Alireza Rangriz, Chalaturnyk, Rick, Zhao, Benzhong, MacMinn, Christopher W., Juanes, Ruben, 2018. Quasistatic fluid-fluid displacement in porous media: Invasion-percolation through a wetting transition. Phys. Rev. Fluids 3 (10), 104001.

Primkulov, Bauyrzhan K., Zhao, Benzhong, MacMinn, Christopher W., Juanes, Ruben, 2022. Avalanches in strong imbibition. Commun. Phys. 5 (1), 52.

Raeini, Ali Q., Bijeljic, Branko, Blunt, Martin J., 2018. Generalized network modeling of capillary-dominated two-phase flow. Phys. Rev. E 97 (2), 023308.

Rayleigh, Lord, 1878. On the instability of jets. Proc. Lond. Math. Soc. 1 (1), 4-13.

Rieser, Jennifer M., Arratia, P.E., Yodh, A.G., Gollub, J.P., Durian, D.J., 2015. Tunable capillary-induced attraction between vertical cylinders. Langmuir 31 (8), 2421-2429. http://dx.doi.org/10.1021/la5046139, PMID: 25646573.

Romano, Marta, Chabert, Max, Cuenca, Amandine, Bodiguel, Hugues, 2011. Strong influence of geometrical heterogeneity on drainage in porous media. Phys. Rev. E 84 (6), 065302.

Ryazanov, Andrey V., Van Dijke, Marinus Izaak Jan, Sorbie, Kenneth Stuart, 2009. Two-phase pore-network modelling: Existence of oil layers during water invasion. Transp. Porous Media 80, 79-99.

Scheel, M., Seemann, R., Brinkmann, M., Michiel, M. Di, Sheppard, A., Herminghaus, S., 2008. Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime. J. Phys.: Condens. Matter 20 (49), 494236. http://dx.doi.org/10.1088/0953-8984/20/49/494236.

Su, Junwei, Wang, Le, Gu, Zhaolin, Zhang, Yunwei, Chen, Chungang, 2018. Advances in pore-scale simulation of oil reservoirs. Energies 11 (5), 1132.

Sun, Zhonghao, Santamarina, J. Carlos, 2019. Haines jumps: Pore scale mechanisms. Phys. Rev. E 100 (2), 023115.

Surasani, V.K., Metzger, T., Tsotsas, E., 2008. Influence of heating mode on drying behavior of capillary porous media: Pore scale modeling. Chem. Eng. Sci. 63 (21), 5218-5228.

Tahmasebi, Pejman, Sahimi, Muhammad, Kohanpur, Amir H., Valocchi, Albert, 2017. Pore-scale simulation of flow of CO₂ and brine in reconstructed and actual 3D rock cores. J. Pet. Sci. Eng. 155, 21-33.

Tuller, Markus, Or, Dani, 2001. Hydraulic conductivity of variably saturated porous media: Film and corner flow in angular pore space. Water Resour. Res. 37 (5), 1257-1276.

Vahid Dastjerdi, Samaneh, Karadimitriou, Nikolaos, Hassanizadeh, S. Majid, Steeb, Holger, 2022. Experimental evaluation of fluid connectivity in two-phase flow in porous media during drainage. Water Resour. Res. 58 (11), e2022WR033451.

Valvatne, Per H., Blunt, Martin J., 2004. Predictive pore-scale modeling of two-phase flow in mixed wet media. Water Resour. Res. 40 (7).

Vincent-Dospital, Tom, Moura, Marcel, Toussaint, Renaud, Måløy, Knut Jørgen, 2022. Stable and unstable capillary fingering in porous media with a gradient in grain size. Commun. Phys. 5 (1), 306.

Vorhauer, Nicole, Wang, Y.J., Kharaghani, A., Tsotsas, Evangelos, Prat, Marc, 2015. Drying with formation of capillary rings in a model porous medium. Transp. Porous Media 110, 197-223.

Wilkinson, David, 1984. Percolation model of immiscible displacement in the presence of buoyancy forces. Phys. Rev. A 30 (1), 520.

Wilkinson, David, Willemsen, Jorge F., 1983. Invasion percolation: A new form of percolation theory. J. Phys. A: Math. Gen. 16 (14), 3365.

Wu, Ting, Yang, Zhibing, Hu, Ran, Chen, Yi-Feng, 2023. Three-dimensional visualization reveals pore-scale mechanisms of colloid transport and retention in two-phase flow. Environ. Sci. Technol. 57 (5), 1997-2005.

Xu, Wei, Ok, Jeong Tae, Xiao, Feng, Neeves, Keith B., Yin, Xiaolong, 2014. Effect of pore geometry and interfacial tension on water-oil displacement efficiency in oil-wet microfluidic porous media analogs. Phys. Fluids 26 (9), 093102.

Zhai, Qian, Rahardjo, Harianto, Satyanaga, Alfrendo, Dai, Guoliang, 2020. Estimation of the soil-water characteristic curve from the grain size distribution of coarse-grained soils. Eng. Geol. 267, 105502.

Zhao, Xiucai, Blunt, Martin J., Yao, Jun, 2010. Pore-scale modeling: Effects of wettability on waterflood oil recovery. J. Pet. Sci. Eng. 71 (3-4), 169-178.

Zhao, Benzhong, MacMinn, Christopher W., Juanes, Ruben, 2016. Wettability control on multiphase flow in patterned microfluidics. Proc. Natl. Acad. Sci. 113 (37), 10251-10256.

Zhao, Benzhong, MacMinn, Christopher W., Primkulov, Bauyrzhan K, Chen, Yu, Valocchi, Albert J., Zhao, Jianlin, Kang, Qinjun, Bruning, Kelsey, McClure, James E., Miller, Cass T., et al., 2019. Comprehensive comparison of pore-scale models for multiphase flow in porous media. Proc. Natl. Acad. Sci. 116 (28), 13799-13806.

Zhao, Jianlin, Qin, Feifei, Kang, Qinjun, Derome, Dominique, Carmeliet, Jan, 2022. Pore-scale simulation of drying in porous media using a hybrid lattice Boltzmann: Pore network model. Drying Technol. 40 (4), 719-734.

Zhou, Dengen, Blunt, Martin, Orr, Jr., F.M., 1997. Hydrocarbon drainage along corners of noncircular capillaries. J. Colloid Interface Sci. 187 (1), 11-21.

Zhu, Guangpu, Yao, Jun, Li, Aifen, Sun, Hai, Zhang, Lei, 2017. Pore-scale investigation of carbon dioxide-enhanced oil recovery. Energy Fuels 31 (5), 5324-5332.

Zhu, Lijun, Zhang, Heng, Xiao, Liusheng, Bazylak, Aimey, Gao, Xin, Sui, Pang-Chieh, 2021. Pore-scale modeling of gas diffusion layers: Effects of compression on transport properties. J. Power Sources 496, 229822.