# Interaction between corner and bulk flows during drainage in granular porous media

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

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

► To cite this version:

Paula Reis, Gaute Linga, Marcel Moura, Per Arne Rikvold, Renaud Toussaint, et al. Interaction between corner and bulk flows during drainage in granular porous media. Advances in Water Resources, 2025, 198, pp.104914. 10.1016/j.advwatres.2025.104914. hal-05402343

HAL Id: hal-05402343

https://univoak.hal.science/hal-05402343v1

Submitted on 7 Dec 2025

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 Attribution - ShareAlike 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# Interaction between corner and bulk flows during drainage in granular porous media

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

## ARTICLE INFO

### Keywords:

Drainage

Granular porous media

Corner flow

Pore-scale modeling

Capillary bridge

## ABSTRACT

Drainage of a liquid by a gas in porous media can be broken down into two main mechanisms: a primary piston-like displacement of the interfaces through the bulk of pore bodies and throats, and a secondary slow flow through corners and films in the wake of the invasion front. In granular porous media, this secondary drainage mechanism unfolds in connected pathways of pendular structures, such as capillary bridges and liquid rings, formed between liquid clusters. To represent both mechanisms, we proposed a dynamic dual-network model for drainage, considering that a gas displaces a wetting liquid from quasi-2D granular porous media. For this model, dedicated analyses of the capillary bridge shapes and hydraulic conductivity were conducted so that the secondary drainage mechanism could be properly quantified at finite speeds. With the model, an investigation of the wetting-phase connectivity and flow during drainage was carried out, covering a broad range of flow conditions. Results indicate that the span of liquid-connected structures in the unsaturated region, as well as their ability to contribute to flow, varies significantly with Capillary and Bond numbers.

## 1. Introduction

The displacement of a liquid by a gaseous phase in granular porous media is a central component of diverse natural and engineered sub-surface processes, such as CO₂ geological sequestration (Juanes et al., 2010; Aryana and Kovscek, 2012; Datta and Weitz, 2013; Zacharoudiou et al., 2018), enhanced oil recovery (Alvarado and Manrique, 2010; Afzali et al., 2018), and water retention in soils (Or and Assouline, 2013; Hoogland et al., 2016a,b, 2017). Understanding and properly designing these processes, therefore, depend on how sensibly the pore-scale dynamics of liquid mobilization by gas can be characterized and upscaled to the geological length scale of interest (Toussaint et al., 2012; Moebius and Or, 2014). In this work, we focus on the pore-scale description of such flows under the assumption that the liquid phase preferentially wets the grains constituting the porous medium and that the viscosity and density of the gas can be neglected. In particular, we investigate how the wetting-phase connectivity and flow through corners during drainage can be effectively represented in pore-scale models, leading to more realistic predictions of transport properties and liquid residual saturations under a broad range of flow conditions.

Drainage displacements are commonly sorted into different categories according to the balance between capillary, viscous, and gravitational forces taking place during the flow (Lenormand, 1990). In the absence of significant viscous or gravitational forces, the motion of gas as it invades a liquid-filled random porous medium delineates a pattern known as capillary fingering (Lenormand and Zarcone, 1989; Måløy et al., 1992). Under this regime, the largest pore throats at the fluid's interface are drained one at a time, since they exhibit the lowest values of capillary pressure thresholds for invasion, $P_t$, given by the Young-Laplace Eq. (1). Following this criterion only, the interface advancement traps clusters of the wetting phase and forms fractal fingers that can be well described by invasion percolation theory (Wilkinson and Willemse, 1983; Lenormand and Zarcone, 1989).

$$P_t = \gamma \left( \frac{1}{r_1} + \frac{1}{r_2} \right) \qquad (1)$$

where $\gamma$ is the interfacial tension between the phases and $r_1$ and $r_2$ are the principal radii of curvature of the interface moving through a pore throat.

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

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

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

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

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

<footnote>* Corresponding author.
E-mail address: paula.reis@mn.uio.no (P. Reis).</footnote>
Departing from capillary-dominated to faster flows, the pattern outlined by the gas motion is modified by the onset and evolution of viscous instabilities, leading to viscous fingering (Homsy, 1987). In this case, the progression of the interface is not governed solely by the size of the available pore throats but also by the pressure drop across them (Måløy et al., 1985). This altered criterion for pore invasion is a consequence of the viscous pressure gradient developed in the liquid, which descends along the direction of the flow and favors drainage at the most advanced points of the front. As a result, the efficiency of liquid mobilization is low, and the gas tends to flow through thin fingers aligned with the flow direction (Løvoll et al., 2004, 2011).

In contrast to fingered flows, the invasion front between gas and liquid during drainage may be stable when gravitational effects come into play, with the denser liquid phase lying initially under the gaseous phase (Løvoll et al., 2005). As proposed by Méheust et al. (2002), the crossover for viscous instability in the presence of gravity occurs when the generalized bond number, $Bo^* = Bo - Ca$, falls below zero. Thus, when the bond number – given by $Bo = \Delta\rho g a^2 / \gamma$, where $\Delta\rho$ is the density difference between the phases, $g$ is the gravitational component in the direction of the flow, and $a$ is the typical pore size – exceeds the capillary number – given by $Ca = v\mu a^2 / \gamma k$, where $v$ and $\mu$ are the Darcy velocity and the viscosity of the wetting phase, and $k$ is the porous medium permeability – the hydrostatic pressure gradient in the liquid overcomes the destabilizing viscous pressure gradient, suppressing the unbounded growth of fingers and favoring the invasion of pores far from the tip of the front.

The preceding description of drainage flow regimes – capillary fingering, viscous fingering, and stable front displacement – focuses on the evolution of the main invasion front as it percolates through a porous medium. This front defines the border between a saturated region, labeled here as the main defending cluster, and an unsaturated region where liquid can be found in clusters, corners, and films (Hoogland et al., 2016a; Moura et al., 2019). While the role of corners in wetting-phase transport has been extensively investigated in imbibition (Zhao et al., 2016; Golmohammadi et al., 2021; Primkulov et al., 2021, 2022; Liu et al., 2022; Zhao et al., 2022; Pavuluri et al., 2023; Suo, 2024), few studies have been dedicated to adequately representing drainage via corners in the wake of the main invasion front progression. Next, we present a brief review of studies concerning the role of corner flow in drainage. Then, we establish how the present work fits in with the current literature on this topic.

In a study concerning water retention in soils, Hoogland et al. (2016a) proposed a distinction between two main drainage mechanisms: a rapid piston-like displacement at the front, and a slow gravity-driven corner flow behind the front. In that work, the authors used invasion-percolation theory and experiments on small columns filled with sand and glass beads to predict the saturation at which the network of water-filled pores is disrupted and drainage is dominated by corner flow. In addition, they suggested a star-shaped pore geometry to estimate water retention and transport in natural pores. In an accompanying study (Hoogland et al., 2016b), the foam drainage equation (FDE) (Or and Assouline, 2013) was employed to predict macroscopic drainage dynamics in the slow corner-flow zone behind the front. The alternative use of the FDE was based on the idea that the flow through channels formed between bubbles in foams, known as Plateau borders, could represent the flow in crevices between soil particles. As a result, the wetting phase is assumed to be fully connected by liquid-filled corners behind the front.

The hypothesis of overall wetting-phase connectivity in unsaturated granular porous media can, however, lead to overestimated water removal rates and inaccurate predictions of liquid residual saturations (Bryant and Johnson, 2003). For this reason, conceptualizing flow in a network of angular pores and throats with interconnected liquid-filled corners may be inadequate. In fact, when the flow through the bulk of pores and throats is disrupted in granular porous media, pathways of corner flow may remain, where liquid rings formed at

grain–grain contact points and capillary bridges join disconnected elements of the network of pores and throats. These corner-flow pathways evolve with time and do not coincide with the original network of pores and throats (Moura et al., 2019). In both quasi-2D (Aursjø et al., 2014; Chen et al., 2018; Moura et al., 2019) and 3D (Scheel et al., 2008; Kharaghani et al., 2021; Dong et al., 2022) grain packings, it has been verified in visualization experiments that liquid located in rings around contact points between smooth solid surfaces and capillary bridges can form large locally connected structures behind the front. Still, these do not span the entire unsaturated zone, allowing the drainage of only a fraction of liquid clusters in this region.

Such description of wetting-phase connectivity has been represented in a few pore-scale models for drying (Vorhauer et al., 2015; Chen et al., 2018; Kharaghani et al., 2021), contaminant transport (Chen and Guo, 2023), and drainage (Bryant and Johnson, 2003; Reis et al., 2023), in granular media. Vorhauer et al. (2015) proposed a pore-network model (PNM) of drying in quasi-2D regular porous media in which liquid clusters could be formed not only by interconnected pores and throats in the network but also by adjacent liquid rings. This approach was later expanded to 3D porous media consisting of simple cubic lattices of spheres by Kharaghani et al. (2021), leading to better drying rate estimates than PNM with clusters of pores and throats only. Similarly, Chen et al. (2018) proposed a numerical model to reproduce drying in a circular quasi-2D porous medium composed of cylinders arranged as neighboring spirals between parallel plates. In this particular geometry, no clusters of liquid-filled pores are formed and long chains of capillary bridges connect the central part of the circle to its open outer rim, enhancing the drying rate. Based on this simplified corner-flow geometry, their model was the first to simulate viscous flows through chains of capillary bridges, by employing viscous resistance values obtained with direct numerical simulation (DNS) of Stokes flow on single bridge shapes. This, combined with realistic bridge breakup pressure estimates, enabled their model to reasonably predict the maximum capillary bridge chain length for different imposed evaporation rates.

The effects of pendular structures on PFAS (per- and polyfluoroalkyl substances) transport in granular media were also recently investigated by Chen and Guo (2023), using a PNM. In their model, water occupies the pore network both as thin films covering the grains, and as rings at grain–grain contact points. Using a threshold capillary pressure for disconnecting the pendular rings at drier conditions, their results suggest that the presence of such structures is a primary factor controlling PFAS transport in unsaturated soil. As for models representing drainage, Bryant and Johnson (2003) developed a quasi-static PNM based on a 3D dense packing of spheres in which pendular rings could act as connecting points for liquid in neighboring pores/throats. Although their model only included short-range connectivity rules that could not account for the whole unsaturated region, estimates of residual liquid saturation matched experimental results significantly better than models assuming unrestricted liquid connectivity. With a similar goal, Reis et al. (2023) proposed a modified invasion-percolation (IP) model for slow drainage under the effects of capillary and gravitational forces. Their model could identify long-range extra connectivity provided by liquid rings and bridges, qualitatively reproducing experiments presented by Moura et al. (2019) in quasi-2D porous media consisting of spherical beads trapped between parallel planes.

In the present work, the wetting-phase connectivity concept presented in Reis et al. (2023) was incorporated in a dynamic PNM for drainage, in which a dual-lattice network approach was used to represent both networks of connected pores and throats, and of connected rings and capillary bridges. The interest in moving from an IP model to a dynamic PNM is twofold: investigating the effect of moderate/high Ca on corner flow, and more adequately representing corner flow during slow drainage. Using a quasi-static model, we could not differentiate the time scales required to drain a pore belonging to the front from
a pore belonging to a cluster connected by corner flow. To further improve this point, viscous resistance in bridges was calculated using DNS of Stokes flow through various single bridge shapes, similarly to the work of Chen et al. (2018), and implemented in the network of corner flow. Besides the new model, we present a comprehensive evaluation of how the balance between capillary, viscous, and gravitational forces influence the wetting-phase connectivity and flow in the unsaturated region trailing the main drainage front in granular porous media.

## 2. Pore-scale model

The pore-scale model proposed in this work to represent drainage with corner flow belongs to the category of Dynamic Pore-Network Models (Joekar-Niasar and Hassanizadeh, 2012). In these models, the pore space is conceptualized as networks of geometrically-simplified elements, in which fluid configurations are presumed according to flow conditions and wetting properties. Based on these simplifications, multiphase flow through the networks is normally solved by imposing mass conservation on nodes (representing pores) and assuming Poiseuille flow on the edges (representing throats). Capillary effects are incorporated into the calculations using the Young–Laplace equation (see Eq. (1)) on network elements containing more than one phase. Gravitational effects can also be included, by computing the individual contribution to the hydrostatic pressure of the fluids occupying each network element. While Dynamic PNMs are not able to resolve geometric details related to the fluids' interfaces and pore shapes, as achieved with DNS models, the possibility of simulating bigger domains can offer insights into flow phenomena stretched over several thousands of pores (10³ – 10⁴), as explored in this study.

Within the dynamic PNM category, the one presented here was specifically designed to reproduce drainage in a quasi-2D porous medium consisting of a Hele-Shaw cell filled with a monolayer of glass beads. This design of porous medium has been employed in several experimental studies of drainage (Måløy et al., 1985, 1992; Méheust et al., 2002; Løvoll et al., 2004, 2005, 2011; Toussaint et al., 2012; Moura et al., 2015, 2019; Ayaz et al., 2020) and represented in a quasi-static PNM in Reis et al. (2023). Such structure is particularly beneficial to our study, as it allows us to clearly distinguish between flow phenomena in saturated and unsaturated regions while retaining fundamental 3D features of corner flow in granular materials. Behind the main invasion front, liquid rings remain between each bead and the Hele-Shaw cell planes, linking otherwise disconnected liquid-filled pores and throats. Notably, draining such a porous medium evidences how the pathways of corner flow do not coincide with the main network of pores and throats (Moura et al., 2019), a key aspect explored in this work and detailed in Section 2.1.

Another particular aspect of the proposed model comes from the choice of represented wetting and non-wetting fluids, namely a glycerol-water mixture and air. Given their significant difference in viscosity and density, the pressure gradient within the non-wetting phase is assumed to be negligible. Further details of how this assumption affects the flow calculation in each network element, as well as the model's governing equations, are discussed in Sections 2.2 and 2.3. As a final noteworthy feature of the presented PNM, estimates of the viscous resistance through capillary bridges were obtained with DNS of Stokes flow through a wide range of bridge shapes, as discussed in Section 2.2.2.

### 2.1. Dual-lattice network representing pores, throats, and corner-flow pathways

The representation of Hele-Shaw cells filled with glass beads as networks of pores and throats was carried out here in the same way as proposed by Reis et al. (2023). The disposition of the beads between the planes was assumed to form a regular triangular lattice (TL), leading to a hexagonal lattice (HL) of pores connected by throats. During drainage

simulations, starting from a saturated pore space, the original pore network with HL topology evolved according to its fluid occupation. Nodes representing pores in this lattice could be labeled as filled by gas or liquid, while edges representing throats were deleted from the original grid when completely filled by gas. In this way, the remaining connected parts of the hexagonal lattice network correspond directly to the liquid clusters in the porous medium. More details on this network approach to represent fluids in pores and throats, as well as its relation to graph theory concepts, can be found in Reis et al. (2023).

Using the hexagonal-lattice network as a scaffold, the flow through the bulk of pores and throats was represented in the model. As proposed in Hoogland et al. (2016a) and Moura et al. (2019), the fast piston-like flow of gas as it invades the bulk of liquid-filled pores and throats is responsible for mobilizing most of the drained volume, and is therefore termed as the primary drainage mechanism. In opposition, the slow corner-flow dominated drainage in the unsaturated region (Hoogland et al., 2016a) is termed the secondary drainage mechanism (Moura et al., 2019). Besides the difference in time scales required for drainage via both mechanisms, it has been verified experimentally that they also do not happen in coinciding paths in granular media. For this reason, we chose to represent corner flow through an evolving secondary network built on a triangular lattice – the dual lattice of the hexagonal pore network.

The schematic representation of a Hele-Shaw cell filled with beads undergoing drainage in Fig. 1 helps illustrate the dual-lattice network approach. In the sketched image, depicting a top view of the porous medium, large gray circles represent the beads and blue areas represent the liquid occupation of the porous medium, at the frontier between saturated and unsaturated zones. Along with the sketch, the corresponding dual-lattice network representation of the wetting-phase connectivity is shown. From Figs. 1(a) to 1(c), we see the non-wetting phase invading two new pores and throats, and the resulting alterations in the network representing pores connected by throats. At this point, the hexagonal-lattice (HL) network still corresponds to the region connected by liquid. From Figs. 1(d) to 1(g), further advancement of the gas into the porous medium leaves behind capillary bridges between adjacent beads, which can establish permeable enough paths for the liquid to flow (Chen et al., 2017, 2018; Moura et al., 2019). These flow paths are perpendicular to the ones in the hexagonal lattice, and could not be represented in the previously established network of pores connected by throats. Therefore, to account for the flow through the capillary bridges, we created new nodes and edges in the dual lattice of the original hexagonal pore network. In the new lattice, edges (double lines in Fig. 1) represent the capillary bridges, while the nodes (green diamond shapes in Fig. 1) are the bridges' endpoints, coinciding with the glass beads' positions. As dual lattices do not share nodes or edges, connecting elements were also created, so that the liquid flow could be represented in a single network. In this way, the new links, illustrated in Fig. 1 as dashed lines, represent the flow between liquid clusters and the capillary bridges. As a final remark, the flow through capillary bridges may also be interrupted, if the capillary pressure for the breakup of a bridge is reached, as illustrated in Fig. 1(h). In this case, an isolated liquid ring can be formed between the bead and each Hele-Shaw cell plane, as represented by the blue-filled circle on top of the bead. Additional details about the breakup process are given in Section 2.2.3.

Fig. 2 represents a more complex liquid configuration on a larger domain, which includes the main defending cluster, a cluster connected to the front by a chain of two capillary bridges, a small disconnected cluster, a disconnected bridge, and several isolated liquid rings. In this figure, an experimental photo Moura et al. (2019) from a similar liquid configuration can also be seen in the bottom-left insert.

### 2.2. Flow through different network elements

The dual-lattice network described in Section 2.1 represents the wetting-phase occupation in the modified Hele-Shaw cell porous
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Illustration of the dual-network representation of flow paths, as the quasi-2D granular porous medium is drained. Network elements are defined in the bottom right insert. Small circle and diamond shapes represent the network nodes in the hexagonal and triangular lattices, respectively, and correspond to the endpoints of pore throats and capillary bridges. Simple, double, and dashed lines represent the network edges, representing pore throats, capillary bridges, and the connections between them, in that order. From (a) to (c), gas fully invades two consecutive pore throats, leading to the deletion of their representative edges from the hexagonal-lattice network. From (d) to (g), the movement of the gas-liquid interface leads to the formation of capillary bridges, which are then represented with new nodes and edges on the dual-lattice of the network of pores and throats. From (g) to (h) the transformation of the network due to the breakup of a bridge is illustrated.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Schematic representation of saturated and unsaturated zones during drainage in a Hele-Shaw cell filled with glass beads and the corresponding flow networks. Large gray circles represent glass beads, and blue areas correspond to the liquid occupation of the pore space. As for the network representation: on the hexagonal lattice, small blue-colored circles represent liquid-filled pores, while white ones are filled with gas. Edges representing throats only exist if filled with liquid. As seen in the bottom-left insert, clusters connected by capillary bridges are disconnected in this network. On the triangular lattice, edges representing capillary bridges establish pathways for corner flow.</figcaption></figure>

medium, schematically illustrated in **3**. During drainage, liquid can only flow through the portions of this network connected to the liquid outlet, with the rest being permanently trapped. During this process, the pressure in the air-invaded portions of the porous medium is assumed to remain constant, due to this phase's relatively low density and viscosity.

A converging-diverging shape was adopted in this work for the pore throats (see Section 2.2.1) (Aker et al., 1998). In this way, we do not specify a particular shape or flow resistance to the pores, and nodes in the networks act as mere throat connecting points. Therefore, to calculate the liquid flow out of the porous medium we only need
<figure><img src="image_6.png" /><figcaption>Fig. 3. Schematic representation of a tilted rectangular Hele-Shaw cell filled with a monolayer of spherical glass beads. X and Y are the global coordinates inside the cell, and β is the inclination angle of the cell with respect to the horizontal plane. The air inlet is positioned at the top of the cell, while the liquid outlet is positioned at the bottom of the cell. In this way, gravitational pressure fields stabilize the flow.</figcaption></figure>

to solve the flow through edges, which can represent four different scenarios: (i) liquid-filled throats located inside clusters, (ii) throats containing air and liquid separated by a meniscus at the clusters' borders with the non-wetting phase, (iii) capillary bridges, and (iv) the junctions between capillary bridges and clusters. Altogether, these liquid structures compose a single network that connects to a liquid outlet, where the flow is locally governed by Darcy's law.

In Section 2.2.1, we describe the flow calculation through pore throats, represented by edges in the hexagonal lattice. Then, in Sections 2.2.2 and 2.2.3, we detail the procedure used to calculate the flow through capillary bridges, represented by edges in the triangular lattice, and the criterion for capillary bridge breakup. Finally, in Section 2.2.4, the flow between capillary bridges and clusters is detailed.

### 2.2.1. Flow through pore throats

In the proposed model, pore throats formed between adjacent glass beads in the Hele-Shaw cells were idealized as converging-diverging cylinders with parabolic profiles given by Eq. (2).

$$r(x') = a + b \left(x' - \frac{l}{2}\right)^2 \qquad (2a)$$

$$a = r_t \qquad (2b)$$

$$b = r_b \left(\frac{2}{l}\right)^2 \qquad (2c)$$

where l is the length of the throat, $r_t$ is the radius at the throat's constriction, corresponding to half the distance between the surfaces of adjacent beads, and $r_b$ is the radius of the beads. Values of l, $r_b$, and the distribution of $r_t$ were based on the experimental work of Moura et al. (2019) and are the same as used in Reis et al. (2023).

With the profile given by Eq. (2), the flow through pore throats is calculated assuming Poiseuille flow and varies according to the fluid occupation status. At liquid-filled throats, the contributions of viscous and gravitational effects to the pressure are considered. At throats containing both phases, as illustrated in Fig. 4, capillary pressure effects are also taken into account. Next, we present the procedure to calculate the flow through throats containing a meniscus separating two phases. From that, single-phase flow can be obtained as a particular case of the presented equations, in which the meniscus position $x'_m = 0$ and the capillary pressure $\Delta P_{cap} = 0$.

Using Fig. 4 as a guide, we assume that the pressure drop across the pore throat, $\Delta P = P_i - P_{i+1}$, can be decomposed into three components, $\Delta P = \Delta P_{visc} + \Delta P_{grav} + \Delta P_{cap}$, in which $\Delta P_{visc}$ represents the pressure variation due to viscous dissipation during flow, $\Delta P_{grav}$ represents the hydrostatic pressure and $\Delta P_{cap}$ represents the capillary pressure. Each of these components of the pressure drop is calculated according to Eqs. (3):

$$\Delta P_{visc} = \frac{8q}{\pi} \left( \int_{0}^{x'_{m}} \frac{\mu_{g}}{r(x')^{4}} dx' + \int_{x'_{m}}^{l} \frac{\mu_{l}}{r(x')^{4}} dx' \right) \qquad (3a)$$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Schematic representation of a converging-diverging pore throat containing a meniscus between the gas and liquid phases.</figcaption></figure>

$$\Delta P_{grav} = \rho_l g (h_{i+1} - h_i) \left(\frac{l - x'_{m}}{l}\right) \qquad (3b)$$

$$\Delta P_{cap} = \gamma \left( \frac{\cos \theta_b}{r(x'_m)} + \frac{\cos \theta_p}{r_b} \right) \qquad (3c)$$

where q is the flow rate through the throat, $x'_m$ is the meniscus position, $\mu_g$ and $\mu_l$ are the viscosities of gas and liquid, respectively, $\rho_l$ is the density of liquid, g is the gravitational acceleration, $h_{i+1}$ and $h_i$ are elevations at the throat's endpoints, calculated as $(Y_{i+1} - Y_i) \sin \beta$, with (X, Y) being the global coordinates of the Hele-Shaw cell and $\beta$ its angle of tilt to the horizontal plane, as shown in Fig. 3, $\theta_b$ is the contact angle between the air-liquid interface and the beads' surface, and $\theta_p$ is the contact angle between the air-liquid interface and the Hele-Shaw cell planes.

Based on the decomposition of $\Delta P$ into viscous, gravitational, and capillary components presented in Eq. (3), we define the hydraulic conductance of each throat as $\Gamma = q/\Delta P_{visc}$, and calculate the throat's flow rate as $q = \Gamma(\Delta P - \Delta P_{grav} - \Delta P_{cap})$. This equation for flow through individual throats is later incorporated into the model's governing equations for flow through the networks, as described in Section 2.3.

### 2.2.2. Flow through capillary bridges

To model flow through the network of capillary bridges, we follow a similar approach to Chen et al. (2018) and use DNS to compute the conductance of individual liquid bridges. By inspection of the experimental images, we find that capillary bridges often connect two adjacent beads, as seen in the bottom left insert in Fig. 2, establishing liquid flow paths distinct from those between pore bodies and throats. We thus seek an effective expression for the conductance $\Gamma$ of a link in the dual network representing two beads connected by a capillary bridge, as a function of the bead-bead separation $\ell$ and the capillary pressure $p_c$ in the bridge (defined here as the pressure difference between the outside and the inside of the bridge, hence typically $p_c > 0$).

We assume quasi-static conditions, such that the shape of a bridge is solely given by its equilibrium shape, i.e. the shape that minimizes the free energy of the air-liquid interface (constrained by liquid volume conservation). As such, any flow through the bridge does not alter its shape. This approximation is good in the range of low capillary numbers considered here, which coincides with the range where corner flow can play an important role. We note that the liquid volume of the bridge must be a function of the capillary pressure and the distance $\ell$ between the beads, hence $\Gamma = \Gamma(p_c, \ell)$. By dimensional analysis (see Appendix A), we find

$$\Gamma = \frac{d^4}{\mu \ell} k_0 \left( \frac{p_c d}{2\gamma} \frac{\ell}{d} \right), \qquad (4)$$
<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) Liquid bridge geometry. The quarter of the domain drawn with solid lines is used for CFD calculations by symmetry considerations. The solid beads have centers at (x, y, z) = (0, 0, 0) and (ℓ, 0, 0) and are not shown directly here. (b) Visualization of the simulated flow field in the computational domain showing streamlines with velocity magnitude |u| and pressure p.</figcaption></figure>

where $k_0(\tilde{\kappa}, \tilde{\ell})$ is a dimensionless function that only depends on the shape of the bridge, encoded in the dimensionless curvature $\tilde{\kappa} = \kappa d$ and bead separation $\tilde{\ell} = \ell/d$, both normalized by the bead diameter $d$. The capillary pressure $p_c$ is related to the mean curvature $\kappa$ by the Young–Laplace law, $p_c = 2\kappa\gamma$.

We consider $z$ as the out-of-plane coordinate and a periodic array of beads along $x$ with centers at $y = 0, z = d/2$. Thus the connected bridge and net flow is along the $x$ direction. The liquid bridge geometry, as well as the coordinate system, is shown schematically in Fig. 5(a). In the $y$ direction, the surface will be symmetric about the plane $y = 0$. Moreover, in the $x$ direction, the surface is symmetric about the centers of the beads and the midplanes between two adjacent beads. These symmetry considerations limit our domain to a quarter of the original domain (cf. the part of Fig. 5(a) drawn with black solid lines), allowing for increased numerical resolution in the DNS. Note that because of a small but finite gravity, the domain is not symmetric in the $z$ direction. This symmetry breaking generally has a small effect on the conductance, but significant effects on the connectivity and breakup of bridges, cf. Section 2.2.3.

In order to obtain liquid bridge geometries in which to calculate the conductance by DNS, we used the software Surface Evolver (Brakke, 1992) to find the equilibrium air–liquid interface by surface energy minimization in the described solid geometry. We considered a wide range of different scaled curvatures $\tilde{\kappa}$ and bead–bead distances $\tilde{\ell}$ around the values encountered in the experiments by Moura et al. (2019). The symmetry planes were implemented as neutrally wetting boundaries (90° contact angle, which is mathematically equivalent). We then smoothed the curved liquid–air surface using isotropic explicit remeshing (Hoppe et al., 1993) via PyMeshLab (Muntoni and Cignoni, 2021). This smoothing step was key to achieving well-defined surface normals required for imposing stress-free boundary conditions (see Eqs. (6) below). We verified that the resulting meshes were not significantly altered and still physically valid, e.g. compatible with the imposed contact angles, as reported in Appendix B. The interior of the resulting closed surface meshes were tetrahedralized using TetGen (Hang, 2015).

Assuming creeping flow conditions, the fluid velocity $\mathbf{u}$ and pressure field $p$ inside the liquid bridge are governed by the Stokes equations:
$$ \mu \nabla^2 \mathbf{u} = \nabla p, \quad \nabla \cdot \mathbf{u} = 0, \qquad (5) $$

where $\mu$ is the dynamic liquid viscosity. We decompose the pressure into a linear and a fluctuating part, $p = Fx + p'$, where $F = \Delta p/\ell$ is the imposed macroscopic pressure gradient driving a net flow through the bridge, and $p'$ represents local pressure variations. As for the surface, the velocity field solution is symmetric about $y = 0$. Because of the linearity of the Stokes Eqs. (5) the resulting flow field is antisymmetric about the planes $x = 0$ and $x = \ell/2$. We further assume no-slip conditions at the solid boundaries. This amounts to the following boundary conditions:

$$ u_y = \frac{\partial u_x}{\partial y} = \frac{\partial u_z}{\partial y} = 0 \quad \text{at } y = 0, \qquad (6a) $$

$$ \mathbf{u} = \mathbf{0} \quad \text{at } z = \pm d/2, \text{ and } \sqrt{x^2 + y^2 + z^2} = d/2, \qquad (6b) $$

$$ \frac{\partial u_x}{\partial x} = u_y = u_z = p' = 0 \quad \text{at } x = 0, \ell/2. \qquad (6c) $$

On the air–water interface, we impose a free-slip condition:

$$ \mathbf{u} \cdot \hat{\mathbf{n}} = 0, \qquad (6d) $$

$$ (\mathbf{Du} \cdot \hat{\mathbf{n}}) \times \hat{\mathbf{n}} = \mathbf{0}, \qquad (6e) $$

where $\mathbf{Du} = (\nabla \mathbf{u} + \nabla \mathbf{u}^T)/2$ is the symmetrized strain rate, and $\hat{\mathbf{n}}$ is the normal vector to the air–water interface. Note that Eq. (6d) reflects the quasi-static approximation, i.e. the interface shape is fixed regardless of flow. Eq. (6e) encodes that no shear stress acts along the interface, which follows from neglecting the viscosity in the air phase.

We solve numerically Eqs. (5) and (6) to obtain the shape factor $k_0(\tilde{\kappa}, \tilde{\ell})$ as the total flux $\Phi$ in the case of $\mu = F = d = l$ (see Appendix A). The numerical implementation relying on the finite element method is described in Appendix C in more detail. Note that to obtain $\Phi$ (and thus $k_0$) we must multiply the numerically obtained flux by 2 because of our symmetry-reduced computational domain. Finally, knowing $k_0$, Eq. (4) can be used to retrieve the physical conductance.

An example of numerically obtained velocity and pressure fields is shown in Fig. 5. It is worth noting that the fastest flow occurs in the corners between the beads and the confining plates, i.e. near the top and bottom, while the flow in the near-contact region between two adjacent beads is negligible. This observation is also reflected in the pressure gradients, which are more prominent close to the computational inlet plane, i.e. $x = 0$.

Fig. 6 shows the results for the shape factor $k_0$ obtained by DNS for different bead–bead separations $\ell = d + \delta$ and mean curvatures
<figure><img src="image_9.png" /><figcaption>Fig. 6. Shape factor $k_0$ obtained by DNS for different bead-bead separations $\ell = d + \delta$ and mean curvatures $\kappa$.</figcaption></figure>

$\kappa$. Notably, we find that $k_0$ is independent of $\ell$, evidenced by the numerical data collapse onto a single curve, $k_0(\kappa)$. We find that the function

$$k_0(\bar{\kappa}) = \left( \frac{A}{\bar{\kappa} + B} \right)^3, \quad (7)$$

describes the data remarkably well, with $A = 0.0869$ and $B = 0.417$ determined by a nonlinear fit of Eq. (7) to the data. A physical justification for the functional form of Eq. (7), as well as for the independence of $\ell$, is given in Appendix D. We stress that the values of $A, B$ are specific to the chosen values of the contact angles $\theta_b = 20^\circ, \theta_p = 60^\circ$, but in accordance with the derivation in Appendix D we expect the approximate form of Eq. (7) to be fairly general.

### 2.2.3. Breakup of capillary-bridges

The liquid flow paths established by capillary bridges, as detailed in Section 2.2.2, may be disrupted in case the capillary pressure for the breakup of a bridge is reached during the flow. As the capillary pressure increases, the surface of a bridge becomes highly curved, making its structure progressively more slender. At a critical pressure value, the bridge becomes unstable, leading to the breakup. The effect of capillary pressure on the bridges' shapes is illustrated in Fig. 7. In the two images on the left side (a), we see a lateral and a top view of a capillary bridge modeled with the software Surface Evolver. This bridge is formed between two adjacent spheres, which rest between two parallel planes, representing a small sample of a Hele-Shaw cell filled with glass beads. In these images, the capillary pressure corresponds to 231 Pa, the spheres have a diameter of 1 mm, and are separated by a gap of 0.3 mm. In 7(b), we have the same sphere configuration, but the bridge is at a capillary pressure of 297 Pa. It is noticeable how increasing the capillary pressure made the bridge thinner, which could eventually lead to its rupture, as shown in 7(c). In this third set of images, all the liquid contained at the former capillary bridge is redistributed to the four liquid rings formed between the beads and the Hele-Shaw planes, ceasing the possibility of flow. This scenario helps illustrate how we can have liquid at the corners of granular media while not necessarily establishing an overall wetting-phase continuity.

This phenomenon was observed in experiments conducted by Moura et al. (2019) and could lead to important drainage interruptions in the

wake of the invasion front. It was verified numerically (Reis et al., 2023) that even the breakup of single capillary bridges could disconnect and trap significant portions of liquid in quasi-2D unsaturated granular media. To incorporate these effects in our model, we adopted the same criterion for capillary bridge breakup as proposed in Reis et al. (2023). In this way, for the sake of simplicity, the capillary pressure at the breakup of a bridge in a given throat is approximated by the maximum capillary pressure at the piston-like invasion of the same throat (Eq. (3c), with $x'_m = 1/2$), given the similar curvatures experienced by the gas-liquid interfaces in both scenarios.

### 2.2.4. Flow between capillary bridges and liquid clusters

In Sections 2.2.1 and 2.2.2, the conductances of pore throats – represented as network elements of the hexagonal lattice – and capillary bridges – represented as network elements of the triangular lattice – were defined. In order to have a single network containing elements from both lattices, we also need to define the flow through connecting links between capillary bridges and pore bodies, as shown in Figs. 1 and 2 by the dashed edges.

To estimate the conductance of such links, we leveraged the resemblance of their flow paths to that of a capillary bridge. As suggested in Fig. 5(b) and further discussed in Appendix D, most of the viscous pressure drop over a capillary bridge occurs at the corners around the solid-solid contact points between the beads and the planes. A similar flow path through corners is followed by the liquid as it goes from a capillary bridge to a pore body. Given that the length of these corners in a connecting link is half of that of a capillary bridge, we approximated these links' conductances as twice the conductance described in Eq. (4) and (7) for a bridge.

### 2.3. Governing equations

With the flow through pore throats, capillary bridges, and their connecting links defined in Section 2.2, we can write the equations for the drainage of the networks described in Section 2.1. To do so, first, we define the incidence matrix $C$, relating the network's nodes and edges,
<figure><img src="image_10.png" /><figcaption>Fig. 7. Examples of capillary bridges modeled with the software Surface Evolver. Here, the spherical beads placed between planes, represented by the pink mesh, have a diameter of 1 mm and are separated by a distance of 0.3 mm. The blue regions of the capillary bridge mesh represent the contact surfaces with the Hele-Shaw cell planes, while the black mesh represents the interface with air. We can notice that the contact angle between the liquid and the spheres, $\theta_b = 20^{\circ}$, differs from the contact angle with the cell planes, $\theta_p = 60^{\circ}$. The lateral and top views of the bridge shown in (a) correspond to a capillary pressure of 231 Pa. In (b), the same spherical bead configuration is shown, but the capillary bridge shape corresponds to a capillary pressure of 297 Pa. Higher values of capillary pressure would lead to breakup of this bridge, as shown in (c). In this case, the capillary pressure is reduced to 123 Pa, as the total liquid of the bridge is redistributed into the 4 liquid rings between the planes and beads. Solid and fluid properties described in Section 2 were used in these simulations.</figcaption></figure>

as indicated in Eq. (8).

$$c_{ij} = \begin{cases} 1 & \text{if edge } j \text{ enters node } i \\ -1 & \text{if edge } j \text{ leaves node } i \\ 0 & \text{otherwise} \end{cases} \qquad (8)$$

Using C, the system of equations for flow in the networks is presented in Eqs. (9). Two sets of variables are defined in this system: the pressure $P_i$ in every node $i$ in the network, and the position of the meniscus $x'_{m,j}$, for every throat containing two phases. To relate these variables, volume balance Eq. (9a) is enforced in every node, and an equation for the motion of the menisci Eq. (9b) is written as a function of the throat's flow rate and cross-sectional area. The choice of edge orientation in C does not affect the network flow calculations for throats containing only liquid. Throats containing a meniscus are oriented leaving the gas-filled node and entering the liquid-filled node so that the sign of the capillary pressure is consistent with the flow direction. Eq. (9c) describes how the flow through each throat or capillary bridge is written using the matrix C and concepts detailed in Section 2.2.

$$\sum_{j=1}^{n_v} c_{ij} q_j = 0 \qquad (9a)$$

$$\frac{dx'_{m,j}}{dt} = \frac{q_j}{\pi r (x'_{m,j})^2} \qquad (9b)$$

$$q_j = -\Gamma_j \left\{ \sum_{k=1}^{n_v} c_{kj} \left[ P_k + h_k \rho_l g \left( \frac{l - x'_{m,j} \delta_m}{l} \right) \right] + \delta_m \gamma \left( \frac{\cos \theta_b}{r(x'_{m,j})} + \frac{\cos \theta_p}{r_b} \right) \right\} \qquad (9c)$$

where $q_j$ is the volumetric flow rate through a network edge $j$, $n_v$ is the total number of edges in the network, $n_h$ is the total number of nodes in the network, $\delta_m$ is used to differentiate the calculation of flow through throats containing a meniscus, $\delta_m = 1$, from the flow through liquid-filled throats, capillary bridges, and connecting links, $\delta_m = 0$.

These equations are numerically solved in the proposed model by, first, approximating Eq. (9b) with the Implicit Euler method, then applying the Newton-Raphson method to the resulting system of nonlinear equations, for each time step $dt$. During the drainage simulation, time steps are chosen so that the maximum displacement among the menisci does not exceed approximately 10% of the throat's length.

In order to have networks representing only the liquid occupation of the porous medium, as introduced in Section 2.1, whenever a meniscus reaches the full length of a throat, the network shifts as indicated in

Fig. 8. In this case, the edge representing the throat is deleted, the occupation of its end node is changed from liquid to gas, and new menisci are created at the edges incident on that node. Once the pore-throat network is changed, its structure is checked for the creation of new capillary bridges. The procedure to identify the formation of bridges used here is the same as presented in Reis et al. (2023). If a new bridge is formed, new nodes and edges in the triangular lattice network are created to account for the corner flow (see Fig. 1). At any time step, these bridges may also be deleted, in case the average pressure among its end nodes exceeds the bridge's breakup pressure, defined in Section 2.2.3. The solution of Eqs. (9) combined with the network's updates is successively carried out until the non-wetting phase percolates through the pore space.

## 2.4. Simulation parameters

The networks used to generate the results presented in Section 3 represent Hele-Shaw cells with dimensions 12 cm × 14 cm (width vs. height), filled with a monolayer of glass beads with a diameter equivalent to 1 mm. Presuming a regular arrangement of glass beads (see Section 2.1), this leads to HL networks with 13203 nodes and, initially, 19723 edges. Values of $r_t$ are assigned randomly to the edges, following the pore-size distribution reported in Moura et al. (2019). All edges have the same length of $l = 1$ mm.

As for the simulated drainage conditions, the rectangular porous medium cell is initially saturated with liquid, with the top face open to the air inflow. The bottom face of the cell is connected to a liquid outlet, to which a flow rate is prescribed, while the lateral faces are impermeable. The porous medium cells can also have a tilt angle $\beta$ with respect to the horizontal plane so that the flow is affected by gravitational forces (see Fig. 3). Under these conditions, drainage is carried out until the air percolates the cell.

Therefore, we assume that nodes at the top of the networks are initially filled with air, while all others contain liquid. At the bottom of the network, all nodes are connected by high-conductance edges to a single outlet node, to which a fixed drainage rate is assigned. In the presented results, the following set of flow rates was tested $q_{net} = [10^{-3}, 3.5 \cdot 10^{-2}, 10^{-2}, 10^{-1}, 1, 10]$mL min$^{-1}$. As for the network tilt angle, the following values were used: $\beta = [0, 10, 20, 30, 45, 60, 75, 90]^{\circ}$. In this manner, a large range of capillary and Bond numbers were explored in this study. For each set of tested drainage parameters, 40 different networks were used in the simulations, so that the results reflect the flow conditions, not a particular porous medium realization with randomly distributed pore-throat sizes.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** (a) A throat contains an advancing meniscus located at its midsection. (b) The advancing meniscus reaches the throat's full length $x'_m = l$, prompting a shift in the pore-throat network. (c) In this case, the edge representing the throat is deleted, the occupation of the node changes from liquid to gas, and new menisci are created on the edges incident on the newly gas-filled node, with initial position $x'_m = 0$.</figcaption></figure>

To compare the results with experimental data presented in Moura et al. (2019), the wetting phase represented in the model is a glycerol-water mixture, with viscosity and density equivalent to $\mu_l = 5 \cdot 10^{-2}$ Pa s and $\rho_l = 1.2 \cdot 10^3$ kg m$^{-3}$, respectively. The non-wetting phase is air, within which we considered the effects of viscous and gravitational forces on the pressure field to be negligible. For the calculation of the conductance of throats containing both liquid and gas, nonetheless, a viscosity of $\mu_g = 1.8 \cdot 10^{-5}$ Pa s was used to represent the air, so that numerical problems were avoided due to $g_j \to \infty$ when the meniscus reached the throat's length. Besides, different contact angles are assumed between the fluids' interface and the glass beads, $\theta_b = 20^\circ$, and between the fluids' interface and the Hele-Shaw cell planes $\theta_p = 60^\circ$, which were covered with a contact paper sheet for experiments presented in Moura et al. (2019). These contact angles are static and do not vary spatially, throughout the simulated porous medium.

### 3. Results

As illustrated in Fig. 9, the broad range of drainage flow conditions adopted in this study led to distinct flow regimes, including capillary fingering, stable displacement, and viscous fingering. In the axis of Fig. 9, Bond numbers were calculated as $Bo = \rho_l g \sin \beta a^2 / \gamma$, where *a* is the typical pore size, of 1 mm. Capillary numbers were calculated as $Ca = q_{net} \mu_l a^2 / A \gamma k_{net}$, where *A* is the networks transversal area, of 1.2 cm$^2$, and $k_{net}$ is the average network permeability of $3.8 \cdot 10^{-5}$ cm$^2$. Note that this definition of Ca estimates the viscous pressure drop in a porous medium using its permeability and Darcy's law. It differs from $Ca^* = \mu_l v / \gamma$ (where *v* is the fluid velocity), commonly used in pore-scale studies of multiphase flow, in which the viscous pressure drop in a capillary tube is considered.

In each image, white represents the gas occupation of a pore network at breakthrough, while black represents both the liquid-filled regions and the solid matrix. As proposed by Méheust et al. (2002), we can notice that a stable drainage front was only achieved when $Bo^* = Bo - Ca > 0$. On the top row of flow cases, corresponding to $Bo = 0$, we notice the transition from the capillary fingering type of instability to viscous fingering, as Ca increases from $2.8 \cdot 10^{-5}$ to $2.8 \cdot 10^{-1}$. Unstable flows are also present in the last column, when $Ca = 2.8 \cdot 10^{-1}$, as the gravitational pressure gradient is not sufficient to bound the growth of viscous fingers. These results suggest that the model proposed in Section 2 can appropriately represent the interplay between capillary, viscous, and gravitational forces during drainage in porous media at finite speeds.

In the next sections, detailed images of the cases delineated in red, which correspond to a subset of the stable displacement cases, are shown separately for clarity. These cases correspond to tilt angles of $\beta = [20, 30, 60]^\circ$ and imposed liquid withdrawal rates of $q_{net} = [0.01, 0.1, 1] \text{mL min}^{-1}$.

#### 3.1. Total drainage through capillary-bridge-connected paths

Fig. 10 represents the regions of a pore network that were invaded by gas under different flow conditions, separated by drainage mechanism. Green areas correspond to the primary drainage mechanism, during which pores along the main invasion front are invaded in a fast piston-like displacement. Yellow areas correspond to the secondary drainage mechanism, indicating pores that had belonged to clusters but were drained using connected pathways formed by capillary bridges and liquid rings. All nine cases in this figure represent stable gas displacements, and, as expected, decreasing the Bond number and/or increasing the capillary number led to larger bypassed liquid clusters and reduced primary drainage efficiency. From these clusters, we can verify that a more significant fraction could be later drained via corner flow in slower displacements, when compared to fast drainage, for the same Bond numbers.

The additional time required for clusters to be drained to the front through corners can be inferred from Fig. 11. These images illustrate the same cases shown in Fig. 10, but now the colors, from dark blue to dark red, represent ascending gas invasion times. The horizontal color stratification seen in Fig. 11 indicates the stable movement of the gas front from the top of the pore networks, where the gas inlet is located, to their bottom, from which liquid is withdrawn. Small isolated patches with distinctively later invasion times than their neighbors indicate that, under some flow conditions, clusters can be drained through corners far behind the invasion front. Examples of that are the red patches at the center of Fig. 11(a) and the bright spots at the top of Fig. 11(i). Conversely, secondary drainage mechanism events may also take place predominantly in the vicinity of the main invasion front, as indicated in Fig. 11(g). In this case, a significant part of the liquid removal is attributed to the secondary drainage mechanism, as seen in Fig. 10(g), but no noticeable delay is associated with these events.

The total amount of liquid mobilized from the pore networks through the secondary drainage mechanism, represented as a fraction of the pore volume (PV), is presented in Fig. 12. In this graph, continuous lines represent numerical results, and the dashed line represents experiments performed in Moura et al. (2019). The line colors from dark blue to yellow indicate the capillary number of the displacement, in ascending order. The gray line is associated with numerical results obtained with the modified invasion-percolation (IP) model proposed by Reis et al. (2023), in which the effects of viscous forces are not taken into consideration. It can be verified here that, for the same Bond number, increasing the drainage flow rate can reduce significantly the corner flow, as observed before in Fig. 10. This effect is related to the different time scales associated with the primary and secondary drainage mechanisms. An increase in Ca produces a faster movement of the drainage invasion front, but the gravity-driven drainage from the clusters behind the front remains slow. In this way, clusters connected to the front by paths of capillary bridges are allowed less time to drain before eventually becoming disconnected from the fast-moving front. For low and moderate Bo, even the results obtained with a
<figure><img src="image_12.png" /><figcaption>Fig. 9. Gas occupation of a pore network at breakthrough, for various investigated flow conditions. As observed, stable displacement, capillary, and viscous fingering regimes were covered in our analyses.</figcaption></figure>

very low imposed drainage rate, $Ca \approx 3 \cdot 10^{-5}$, differ considerably from the IP results, highlighting the importance of representing viscous effects in the new PNM proposed in this study. For high Bond numbers, however, connected clusters tend to be drained closer to the front, as seen in Fig. 11(g), reducing the effect of $Ca$ in the secondary drainage mechanism. Consequently, as the influence of gravitational forces increases, both the invasion-percolation and dynamic models lead to similar predictions of corner-flow drainage, given that $Ca \le 3 \cdot 10^{-3}$.

It is also noticeable in Fig. 12 that, for low capillary numbers, the effect of gravitational forces on the amount of liquid drained through corners is non-monotonic. This tendency was also verified in the experiments conducted by Moura et al. (2019), as seen in the graph's dashed and continuous medium blue lines. We attribute this effect to two main reasons. During the stable displacement regime, increasing the Bond number leads to flatter invasion fronts and a smaller fraction of liquid remaining in clusters in the unsaturated zone (Birovljev et al., 1991; Måløy et al., 2021). Then, as previously pointed out in Reis et al. (2023), although gravitational effects drive the drainage through corners in this region, at high $Bo$ there is less liquid to be drained, reducing the importance of the secondary drainage mechanism. On top

of that, we noticed a strong correlation between gravitational effects and the probability of breakup of capillary bridges, as presented in Fig. 13. For all evaluated drainage flow rates, there is a significant increase in the fraction of capillary bridges formed behind the invasion front that reaches the breakup capillary pressure, as $Bo$ increases. With this, the connectivity between liquid clusters and the drainage front is more likely to be interrupted, reducing corner flow. More details on how viscous and capillary forces affect the connectivity between clusters and invasion front are presented in the next section.

## 3.2. Active zone for corner flow

As introduced by Moura et al. (2019), the active zone for corner flow denotes the region trailing the invasion front where corner-flow drainage events are more likely to happen. Here, we measure the width of the area containing clusters connected to the invasion front by paths of capillary bridges as a surrogate for the active zone for corner flow. During drainage, this area undergoes large fluctuations, as the liquid long-range connectivity often relies on single capillary bridges that snap off, and on clusters that are drained (Moura et al., 2019; Reis et al., 2023). Still, for the stable drainage regime, these fluctuations tend to
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Representation of the gas occupation of the pore networks at breakthrough, under various flow conditions. In green, we have pores and throats occupied under the primary drainage mechanism. In yellow, we have the regions of the network only drained due to the liquid connectivity established by capillary bridges and liquid rings.</figcaption></figure>

stabilize around a mean, leading to a well-defined finite region where secondary drainage may occur.

Fig. 14 shows the extent of the liquid-connected zone above the front under different drainage conditions. Unlike Figs. 10 and 11, here we show in color the pores and throats filled with liquid, while the black background represents the gas-filled and solid portions of the porous media. In the images, the main defending cluster is represented in light blue, the pores belonging to the invasion front are highlighted in yellow, and the clusters connected to the front by capillary bridges are shown in pink. Clusters that are disconnected from the front, and therefore trapped, are shown in green. For all nine flow conditions represented in Fig. 14, the depicted snapshots correspond to a time $t = 0.75t_{bl}$, where $t_{bl}$ is the breakthrough time. The comparison between the cases represented in Fig. 14 suggests that variations in the relevance of both gravitational and viscous forces during the flow significantly affect the liquid connectivity in the region trailing the invasion front. In fact, their effects on the active zone extent seem much more pronounced than those on the front width, which are well described in the literature (Birovljev et al., 1991; Auradou et al., 1999; Méheust

et al., 2002; Løvoll et al., 2005; Ayaz et al., 2020; Mâløy et al., 2021; Toussaint et al., 2012; Vincent-Dospital et al., 2022).

As previously reported in both experimental (Moura et al., 2019) and numerical (Reis et al., 2023) slow drainage studies, increasing the Bond number leads to a compaction of the active zone, with connected clusters becoming progressively limited to the vicinity of the drainage front. For the cases with $Ca \approx 3 \cdot 10^{-4}$, this reduction of the active zone with Bo explains why corner-flow drainage events could take place way past the front passage with low Bo, as seen in Fig. 11(a), but not with higher Bo, as in Fig. 11(g). For the cases with $Ca \approx 3 \cdot 10^{-3}$, we also notice the formation of limited regions where clusters are connected to the front by capillary bridges, which decrease with rising gravitational forces. Interestingly, as seen in Fig. 12, the amount of liquid drained through corners with $Ca \approx 3 \cdot 10^{-3}$ increases monotonically with Bo. These results suggest that extensive liquid connectivity in the unsaturated region is not linked to increased secondary drainage. In reality, the flow at low Bo and high Ca accentuates the difference between the characteristic drainage time scales associated with the
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Representation of the time of the gas invasion of each drained pore and throat at breakthrough. From dark blue to dark red, the colors indicate early to late times of invasion. Black regions in the images represent liquid clusters that were not drained.</figcaption></figure>

movement of the front and the corner-dominated flow behind it. Therefore, the conditions leading to the formation of extensive networks of clusters connected by capillary bridges are the same leading to reduced drainage efficiencies.

For the drainage cases with $Ca \approx 3 \cdot 10^{-2}$, shown in Figs. 14(c), (f), and (i), we observe that the liquid-connected zone above the front spans the entire simulated porous-medium domain. These results, therefore, cannot indicate whether the connected region would be reduced as Bo increases, as seen for slower flows. The analysis of these images, however, clarifies the difference in the behavior between the curves representing $Ca < 3 \cdot 10^{-2}$ and $Ca \ge 3 \cdot 10^{-2}$ in Fig. 12. For $Ca < 3 \cdot 10^{-2}$, the liquid amount drained through corners either decreases or stabilizes, as Bo goes from intermediate to high values. As briefly discussed in Section 3.1, although gravity is the main driving force for flow in the unsaturated region, its positive effects on drainage can be counterbalanced by the reduced liquid connectivity span seen in Fig. 14. On the other hand, for $Ca \ge 3 \cdot 10^{-2}$, Fig. 12 points to a linear increase in corner-flow drainage with Bo. With a seemingly constant width of the connected liquid region for all evaluated Bond

numbers, this counterbalancing effect vanishes, and the gravitational influence becomes predominantly positive on the secondary drainage mechanism.

As a way to quantify the active zone width for all investigated stable drainage cases, the mean distance between the pores in the active zone to the pores in the invasion front is presented in Fig. 15(a). For each pore in the connected region, this distance is measured as the minimum Euclidean distance to any point on the front. Then, the measured values for all connected pores and all times $t \ge 0.75t_{bh}$ are averaged. For $Ca < 3 \cdot 10^{-4}$, we can observe that the average distance from the clusters connected by capillary bridges to the front is very small, encompassing a region approximately 10 pores wide. For $Ca \approx 3 \cdot 10^{-3}$, as previously suggested in Figs. 14(b), (e), and (h), the effects of gravitational forces on the active zone are prominent, with the average distance from connected pores to the front decreasing threefold within the range of Bo in which the displacement is stable. Finally, for $Ca \approx 3 \cdot 10^{-2}$, we verify that the measured distances are nearly constant with Bo and consistent with an active zone that spans the entire length of the pore networks used in this study.
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Total corner flow drainage obtained at breakthrough, measured as pore volumes (PV) vs. Bo. Here, all flow conditions tested in this work, quasi-static results from Reis et al. (2023) and experiments from Moura et al. (2019) are represented.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Fraction of capillary bridges that suffered breakup during drainage under different flow conditions.</figcaption></figure>

Additionally, the average distance between pores drained via corner flow and the front was measured and presented in Fig. 15(b). This information, along with the results in Fig. 15(a), can indicate whether drained pores from connected clusters tend to be evenly distributed within the connected region, or concentrate at a certain distance from the front. At low Bond numbers, drainage is more likely to occur in a subset of the connected-region pores closer to the invasion front, especially at higher Ca. At $Bo \approx 6 \cdot 10^{-2}$, for instance, the average distance from secondary drainage events to the front is the same for $Ca \approx 3 \cdot 10^{-2}$ and $Ca \approx 3 \cdot 10^{-4}$, even though the total connected region is much narrower for the slower flow.

Congruently, as Bo increases, the probability of drainage events in pores from the connected region farther from the front becomes higher. This effect is clear in the $Ca \approx 3 \cdot 10^{-2}$ curve in Fig. 15(b). For this case, the width of the liquid-connected region above the front is constant (given the pore-network size limitation), and we can observe that the average distance from secondary drainage events to the front is extended with gravitational forces. For $Ca < 3 \cdot 10^{-2}$, the distance from drainage events to the front decreases with Bo, as the whole connected region also narrows.

## 4. Discussion

The results obtained with the dynamic pore-network model proposed in this study suggest that variations in viscous and gravitational forces acting during drainage substantially impact liquid connectivity and flow in the unsaturated zone. Similar to the effects observed on the front width under the stable displacement regime, the width of the region where clusters can still be drained by corner flow increases with the capillary number and decreases with the Bond number. When it comes to drainage efficiency, increasing Ca reduces not only the amount of liquid mobilized by the primary drainage mechanism – as the displacement gradually shifts from stable to viscous fingering – but also hinders the secondary drainage mechanism. Even though the liquid content and connectivity in the unsaturated zone increase with Ca, the different time scales related to the movement of the front and the drainage through liquid rings and capillary bridges impede substantial liquid mobilization from the clusters. The influence of Bo on the total drainage from the unsaturated zone to the front can also be split into contrasting effects. As Bo increases, the gravity-driven corner flow from the clusters to the front is favored. Still, the formation of smaller liquid
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Representation of the main defending cluster (light blue), invasion front (yellow), clusters connected to the invasion front by capillary bridges and liquid rings (pink), and clusters that lost connectivity to the front (green). From (a) to (i), different flow conditions are represented, at a time $t = 0.75t_{br}$, where $t_{br}$ is the breakthrough time.</figcaption></figure>

clusters and more frequent breakup of capillary bridges reduce the permeability of the corner-flow network, potentially reducing drainage.

Besides the novelty of the results, the model proposed in this study is in itself relevant, as it successfully follows an unconventional approach to representing corner flow in granular porous media. The identification of liquid-connectivity networks decoupled from the pore-throat networks (Bryant and Johnson, 2003; Vorhauer et al., 2015; Chen et al., 2017, 2018; Moura et al., 2019; Kharaghani et al., 2021; Reis et al., 2023) can provide a more physically-realistic image of the phase distribution among the grains during drainage. On top of identifying liquid-connected paths formed by liquid rings and capillary bridges, the adopted dual-network approach, along with the study of bridge shapes and conductances, provides a suitable way to quantify corner-flow drainage. In contrast, previous works in the literature either tracked the liquid-connected structures in the unsaturated region in quasi-static models (Vorhauer et al., 2015; Kharaghani et al., 2021; Reis et al., 2023), or calculated the viscous flow through predictable linear chains of capillary bridges (Chen et al., 2017, 2018).

As for this study's potential applications, the connectivity and flow of liquid in granular media during drainage can be linked, for instance, to water retention in soils. Understanding such processes is relevant for environmental and agricultural applications, as it controls the plant water availability and biological functioning of soils (Hoogland et al., 2016a). Therefore, further studies with the proposed model could focus on obtaining macro water-flow properties in soil for large continuous models. Additionally, understanding water flow and distribution in granular media can have important implications for the transport and mixing of solutes in soils. Although the results presented in Section 3 focused on flow, identifying the span and form of liquid-connected structures in unsaturated soil can be instrumental to understanding how far and how fast nutrients or contaminants can travel within this type of porous medium.

Nevertheless, some simplifying assumptions of the presented model and analyses should be taken into consideration prior to their application. A first important simplification is the adoption of a 2D hexagonal lattice to represent the network of pore bodies and throats. This geometry retains some important aspects of 3D granular porous
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** (a) Average Euclidean distance from pores belonging to clusters connected to the front by capillary bridges. (b) Average Euclidean distance from pores effectively drained due to corner flow to the front.</figcaption></figure>

media, e.g. the coordination number (Dong and Blunt, 2009), and we expect that the type of connectivity established by bridges and liquid rings defined in this work to be valid also in 3D (Bryant and Johnson, 2003; Scheel et al., 2008). Still, the results presented in Section 3 could vary in actual 3D porous media, as the liquid phase in the unsaturated region tends to be more connected than in the 2D case (Wilkinson, 1984). The development of a PNM retaining the key concepts presented in this work, but representing 3D or irregular granular media, would likely face many challenges, in particular with keeping track of the formation of capillary bridges and addressing the greater variety of shapes they could exhibit. Besides, the resulting model could be too computationally expensive for practical use. A potentially feasible next step in this direction is representing the porous medium as a simple cubic pack of spheres, as done in the quasi-static PNM for drying proposed by Kharaghani et al. (2021). In this case, the formation of capillary bridges is still predictable and their shapes would resemble the family of liquid structures identified in the experimental work of Scheel et al. (2008).

The effects of wettability on the prevalence of drainage through capillary bridges and corners were also not investigated in this work. The presented results reflect the flow through a water-wet Hele-Shaw cell ($\theta_p = 60°$) filled with strongly water-wet spherical beads ($\theta_b = 20°$) where no dynamic or spatial variations in contact angle were considered. Different combinations of $\theta_p$ and $\theta_b$ would lead to changes in the bridges' shapes, affecting their stability and transport properties. On top of that, mixed-wetting conditions, commonly found in natural subsurface porous media (Garfi et al., 2022; Irannezhad et al., 2023b,a), could

lead to interesting modifications in the presented results. As presented by Irannezhad et al. (2023b,a), displacement flows in mixed-wetting media produce dramatically varying macroscopic displacement patterns, which could alter the observed flow dynamics between clusters and drainage front through liquid rings and bridges.

Another source of uncertainty in the results could be the simplified shape used to calculate the conductance of the model's pore throats. While the conductances of the capillary bridges were based on DNS of Stokes flow through their shapes, the throats were approximated by converging-diverging circular tubes. This choice allowed us to have analytical solutions for their conductances and capillary pressures, but do not represent the throat's shapes to a high level of accuracy. On top of that, the conductances of the capillary bridges also incorporate some simplifications: we considered that neither the flow through their structures nor the presence of neighboring bridges would lead to distortions in their shapes. In any case, such simplifications are typically adopted in PNM and should not significantly impact the level of predictability expected for this modeling approach.

As a final remark, we considered in this work that liquid in the unsaturated region could only be connected and flow through the main network or pore bodies and throats – as inside liquid clusters – or by the network of capillary bridges and liquid rings. Different studies in the literature point out the additional connectivity provided by crevices on rough grains (Bryant and Johnson, 2003) and adsorbed films covering the entire grain surface (Zhang and Guo, 2024). Yet, such structures are expected to have thicknesses – hence permeabilities – much lower than the liquid flow paths incorporated in our model. As our focus here was to quantify the extra drainage allowed by corner flow within the time scale of the drainage front advancement, we assumed that including the flow through crevices and adsorbed films would not lead to significant deviations in the presented results.

## 5. Conclusion

In this work, we proposed a fully-implicit dynamic pore-network model for drainage in granular media, when the wetting phase is a liquid and the non-wetting phase is an inviscid and weightless gas. In the model, dual-lattice networks for the flow through pores connected by throats and the flow through capillary bridges and liquid rings were adopted. With this approach, we could efficiently reproduce and distinguish drainage events associated with the primary and secondary drainage mechanisms (Hoogland et al., 2016a; Moura et al., 2019; Reis et al., 2023). To appropriately quantify the liquid mobilization through the secondary drainage mechanism, a dedicated analysis of the capillary bridges' shapes and capacity to accommodate flow was conducted and used as input for the pore-network model.

Using the resulting model, we investigated drainage in quasi-2D granular porous media under a broad range of capillary and Bond numbers. In particular, the widths of the regions where clusters are still connected to the receding liquid front, as well as the amount of liquid drained from them, were calculated and examined in Section 3. The results suggested that liquid connectivity and flow in the unsaturated region strongly depend on drainage conditions, challenging the assumption that the wetting phase is always fully connected by corner flow.

This study's findings and new modeling approach can be relevant to natural and engineered processes in the subsurface, such as the water cycle in soils and gas storage in saturated geological formations. Therefore, future developments in the presented line of work should include expanding the model to accommodate irregular and 3D networks, aiming at a more fitting representation of real granular media. The evaluation of wettability effects on the observed results, including intermediate and mixed-wetting systems, would also be a fundamental step towards a better representation of drainage in subsurface porous media. In addition, conducting a more rigorous analysis of the flow and interface displacement through the networks of pores and throats found
in granular media – as done here for the capillary bridges using DNS – could improve the accuracy of the results. It could also be interesting to incorporate in the model the possibility of flow through crevices and adsorbed films on the grains' surface, to investigate if these structures would significantly impact drainage under particular flow conditions. In a similar direction, the effects of evaporation in the unsaturated region could also be considered, especially in the limit of slow flows, when liquid removal by drainage and drying potentially have similar characteristic time scales.

### CRediT authorship contribution statement

**Paula Reis:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Gaute Linga:** Writing – review & editing, Writing – original draft, Software, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Marcel Moura:** Writing – review & editing, Visualization, Funding acquisition, Conceptualization. **Per Arne Rikvold:** Writing – review & editing. **Renaud Toussaint:** Writing – review & editing, Investigation, Conceptualization. **Eirik Grude Flekkøy:** Writing – review & editing, Supervision, Funding acquisition, Conceptualization. **Knut Jørgen Måløy:** Writing – review & editing, Supervision, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### 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.

### Acknowledgments

We acknowledge the support of the University of Oslo, Norway, the Njord Centre, Norway, and the Research Council of Norway through the PoreLab Center of Excellence (project number 262644), the Researcher Projects for Young Talents FlowConn, Norway (project number 324555) and M4, Norway (project number 325819), and the CNRS/University of Oslo IRP D-FFRACT, Norway.

## Appendix A. Scaling the Stokes equations for flow through liquid bridges

We introduce the scaled variables

$$ \tilde{\mathbf{x}} = \mathbf{x}/d, \quad \tilde{\mathbf{u}} = \mathbf{u}/U, \quad \text{and} \quad \tilde{p} = p'/P, \qquad (A.1) $$

where we recall that $p' = p - Fx$ and $F = \Delta p/\ell$ is the imposed macroscopic pressure gradient. Plugging these definitions into the Stokes Eqs. (5), we obtain

$$ \frac{\mu U}{Fd^2} \tilde{\nabla}^2 \tilde{\mathbf{u}} = \frac{P}{Fd} \tilde{\nabla} \tilde{p} - \tilde{\mathbf{x}}, \quad \tilde{\nabla} \cdot \tilde{\mathbf{u}} = 0, \qquad (A.2) $$

where $\tilde{\nabla}$ is the gradient operator in the scaled spatial coordinate $\tilde{\mathbf{x}}$. Choosing the scalings

$$ U = \frac{Fd^2}{\mu}, \qquad (A.3) $$

$$ P = Fd, \qquad (A.4) $$

means that the scaled velocity field $\tilde{\mathbf{u}}$ resulting from solving the scaled Stokes equations,

$$ \tilde{\nabla}^2 \tilde{\mathbf{u}} = \tilde{\nabla} \tilde{p} - \tilde{\mathbf{x}}, \quad \tilde{\nabla} \cdot \tilde{\mathbf{u}} = 0, \qquad (A.5) $$

can be related to the physical field $\mathbf{u}$ by (A.3):

$$ \mathbf{u}(\mathbf{x}) = \frac{Fd^2}{\mu} \tilde{\mathbf{u}} \left( \frac{\mathbf{x}}{d} \right). \qquad (A.6) $$

To quantify the conductance we compute the cross-sectional flux (for arbitrary $x$):

$$ \Phi = \int_{\Omega_x} \hat{\mathbf{x}} \cdot \mathbf{u}(\mathbf{x}) d^2 \mathbf{x} = \frac{Fd^4}{\mu} \int_{\bar{\Omega}_x} \hat{\mathbf{x}} \cdot \bar{\mathbf{u}}(\bar{\mathbf{x}}) d^2 \bar{\mathbf{x}} = \frac{Fd^4}{\mu} k_0. \qquad (A.7) $$

Here, we have identified the dimensionless shape factor $k_0$. Note that if the integral is taken over half the cross section, like the computational domain above, a factor 2 must be included here. The conductance $\Gamma$ used for pore-network modeling is implicitly defined by:

$$ \Phi = \Gamma \Delta P \qquad (A.8) $$

where $\Delta p = F\ell$ is the total bead-bead pressure difference. This yields

$$ \Gamma = \frac{\Phi}{F\ell} = \frac{k_0 d^4}{\mu\ell}. \qquad (A.9) $$

Note that since the shape factor $k_0$ only depends on the shape of the bridge, it may only depend on the length scales involved: the bead diameter $d$, the bead-bead separation $\ell$, and the mean curvature $\kappa$. (Note that we here assume that we can neglect any hysteresis effects that could make $k_0$ multi-valued or history dependent.) From these dimensional quantities we can construct exactly two independent dimensionless quantities, such that the Buckingham $\pi$ theorem (Buckingham, 1914) implies:

$$ k_0 = k_0 \left( \kappa d, \frac{\ell}{d} \right). \qquad (A.10) $$

## Appendix B. Validation of bridge shapes

We verify the shapes obtained by Surface Evolver by fitting a circle to the shape to the air-liquid interface that crosses the inlet plane (i.e. the plane that crosses the bead center; see Fig. D.17 below). The effective contact angles can be then be read off by measuring the intersection of the fitted circle and the top plane or the bead. As shown in the left panel of Fig. B.16, the measured angles are in good agreement with the imposed angle. As a second test, we can check the exact relation (D.7) against the data. As shown in the right panel of Fig. B.16, the agreement with the numerical data is excellent. Furthermore, Eq. (D.8) is shown to yield a very good approximation for the data range considered, establishing confidence in the modeling approximations used in Appendix D.

## Appendix C. Finite element implementation of flow through liquid bridges

We want to solve the dimensionless Stokes Eqs. (5) (with $\mu = F = 1$) with the boundary conditions (6) using the finite element method. The main challenge is to implement the free-slip condition on the air-water interface $\Gamma$, which we opt to enforce weakly using Nitsche's method (Nitsche, 1971; Freund and Stenberg, 1995; Urquiza et al., 2014; Gjerde and Scott, 2022). The corresponding variational problem can be expressed as follows. Find $(\mathbf{u}, p) \in \mathcal{V} \times \mathcal{Q}$ such that for all $(\mathbf{v}, q) \in \mathcal{V} \times \mathcal{Q}$:

$$ \int_{\Omega} [2\mathbf{Du} : \mathbf{Dv} - p\nabla \cdot \mathbf{v} - q\nabla \cdot \mathbf{u}] dV - \int_{\Gamma} [((2\mathbf{Du} - p\mathbf{I}) \cdot \hat{\mathbf{n}}) \cdot \hat{\mathbf{n}})(\mathbf{v} \cdot \hat{\mathbf{n}}) + (((2\mathbf{Dv} - q\mathbf{I}) \cdot \hat{\mathbf{n}}) \cdot \hat{\mathbf{n}})(\mathbf{u} \cdot \hat{\mathbf{n}})] dS + \beta \int_{\Gamma} h^{-1}(\mathbf{u} \cdot \hat{\mathbf{n}})(\mathbf{v} \cdot \hat{\mathbf{n}})dS = \int_{\Omega} v_x dV. \qquad (C.1) $$

Here, the velocity subspace $\mathcal{V}$ encodes the Dirichlet conditions $u_y = 0$ at $y = 0$; $\mathbf{u} = \mathbf{0}$ at $z = \pm 1/2$ and $\sqrt{x^2 + y^2 + z^2} = 1/2$; and $u_y = u_z = 0$ at $x = 0, \ell/2$. The pressure subspace $\mathcal{Q}$ includes the Dirichlet conditions $p = 0$ at $x = 0$ and $x = \ell/2$. For the discretized (into tetrahedral cells) domain $\Omega$, we used the Taylor-Hood combination $P_2(\Omega)^3$ for $\mathcal{V}$ and $P_1(\Omega)$ for $\mathcal{Q}$. For $h$, we used the globally minimal mesh size $h_{\min}$ measured as the smallest distance between any two vertices of a cell. As
<figure><img src="image_19.png" /><figcaption>**Fig. B.16.** Validation of bridge shapes. The symbols correspond to the same simulations as plotted in Fig. 6, and were obtained by the fitting procedure described in Appendix B. (a) Cosine of contact angles at the bead and at the plates versus mean curvature. The solid (dashed) line represents the imposed contact angle at the plates (beads). (b) Distance $s_c$ from the z axis to the center of the fitted circle as a function of the radius of principal curvature. The solid line corresponds to Eq. (D.7) while the dashed line corresponds to the leading-order approximation (D.8).</figcaption></figure>

the penalization parameter we heuristically found $\beta = 100$ to yield good results, but we verified that the results were not sensitive to this choice. The variational form (C.1) with the boundary conditions above was discretized and solved numerically using the FEiCS framework (Logg et al., 2012). The resulting linear system was solved using a built-in direct LU solver. The implementation is available at the Git repository https://github.com/gautelinga/liquid bridge conductance.

## Appendix D. Theoretical modeling of bridge conductance

Here we motivate the functional form for the shape factor $k_0$ used for modeling liquid bridge conductance $\Gamma$.

We consider a bead of radius $R = d/2$ centered at the origin. At an in-plane distance $s$ from the bead center, under (or over) the bead, the thickness $h(s)$ of the liquid corner is given by

$$h = R - \sqrt{R^2 - s^2} \approx \frac{s^2}{d}, \qquad (D.1)$$

where the last expression follows from Taylor expansion to the leading order in $s/R$. This approximation is quite good even for moderate $s/R$, e.g. for a typical maximum value of $s$ where the flow is most constrained, $s = 0.6R$, the approximation gives $h \approx 0.18R$ compared to the exact $0.2R$.

The lubrication approximation for the flux per area in the corner yields

$$q = -\frac{h^3}{12\mu} \nabla p \approx -\frac{s^6}{12\mu d^3} \nabla p. \qquad (D.2)$$

We now assume that the flow is azimuthal in the region under the bead where the flow resistance is largest. Eq. (D.2) then becomes

$$q_{\theta} \approx -\frac{s^5}{12\mu d^3} \partial_{\theta} p. \qquad (D.3)$$

We linearize the pressure in the azimuthal angle $\theta$ such that

$$\partial_{\theta} p \approx -\frac{\Delta p}{2\varphi}, \qquad (D.4)$$

where $\Delta p$ is the total pressure drop (bead-bead) and $\varphi(\ell)$ is an angle over which the pressure drop occurs, which depends weakly on the bead-bead separation $\ell$. We model the latter angle by the linear relation

$$\varphi \sim \frac{\ell}{d}, \qquad (D.5)$$

meaning that the further the beads are apart, the larger will be the angle over which the pressure drop occurs. The total flux $\Phi$ is given by

<figure><img src="image_20.png" /><figcaption>**Fig. D.17.** Top half of the inlet plane geometry, with the wetting corner indicated in light blue color.</figcaption></figure>

radially integrating Eq. (D.4) to some outer radius $s$ (at an arbitrary $\theta$). We skip all numerical prefactors for convenience:

$$\Phi \sim \int_{0}^{s} q_{\theta}(s')ds' \sim \frac{\Delta p}{\mu d^{2}\ell} \int_{0}^{s} s's ds' \sim \frac{\Delta p \cdot s^{6}}{\mu d^{2}\ell}. \qquad (D.6)$$

The remaining task is to relate the outer radius $s$ to the mean curvature $\kappa$. As an intermediate quantity, we will use the largest principal curvature, $1/r$, of the air-liquid interface at the inlet plane, which is sketched in Fig. D.17. Assuming the interface to have a partially circular shape, $r$ is the radius of the corresponding circle. We denote $s_c$ as the in-plane distance from an axis $s$ passing through the bead center, to the center of the circle (see Fig. D.17).

By considering the inlet geometry shown in Fig. D.17, we find an exact relation between $s_c$ and $r$:

$$s_c = \sqrt{(\cos \theta_b + \cos \theta_p)dr + r^2 \sin^2 \theta_p}, \qquad (D.7)$$

where $\theta_b, \theta_p$ are the contact angles between the air-liquid interface and the bead and the plates, respectively. For small radii $r$, we may approximate

$$s \approx s_c \approx \sqrt{(\cos \theta_b + \cos \theta_p)dr} \sim \sqrt{dr}. \qquad (D.8)$$

Inserting the latter expression into Eq. (D.6), we obtain

$$\Phi \sim \frac{\Delta p dr^3}{\mu \ell}. \qquad (D.9)$$

The last step is to relate $r$ to the mean curvature $\kappa$. Rather than modeling the smallest principal curvature $1/R_\parallel$, where $R_\parallel \sim s$ and
<figure><img src="image_21.png" /><figcaption>Fig. D.18. Mean curvature of the air-liquid interface versus the principal curvature at the inlet plane. A linear fit is superimposed at the numerical data. The symbols correspond to the same simulations as plotted in Fig. 6.</figcaption></figure>

using the formula

$$ \kappa = \frac{1}{2} \left( \frac{1}{r} - \frac{1}{R_{\parallel}} \right), \qquad (D.10) $$

we found that our numerical data was remarkably well fitted by the simple linear expression

$$ \kappa = \frac{a}{r} - \frac{\bar{\kappa}_0}{d}, \qquad (D.11) $$

with the dimensionless coefficients $a = 0.342$ and $\bar{\kappa}_0 = 0.312$ (which we expect to depend only on $\theta_b, \theta_p$). The fit is shown in Fig. D.18. Note that (D.10) and (D.11) together suggest that $1/R_{\parallel}$ and $1/r$ are linearly related. Inserting (D.11) into (D.9) yields

$$ \Phi \sim \frac{\Delta pd^4}{\mu c^4} (\bar{\kappa} + \bar{\kappa}_0)^{-3}, \qquad (D.12) $$

where $\bar{\kappa} = \kappa d$ is the nondimensionalized mean curvature. From this expression we can extract the conductance $\Gamma = \Phi/\Delta p$ and thereby, upon comparison with (4), the shape factor $k_0$:

$$ k_0 \sim (\bar{\kappa} + \bar{\kappa}_0)^{-3}. \qquad (D.13) $$

This motivates the expression

$$ k_0(\bar{\kappa}) = \left( \frac{A}{\bar{\kappa} + B} \right)^3, \qquad (D.14) $$

used to fit, with A, B as fitting parameters, our numerical data for liquid bridge conductances.

## Data availability

Data will be made available on request.

## References

Afzali, S., Rezaei, N., Zendehboudi, S., 2018. A comprehensive review on enhanced oil recovery by water alternating gas (WAG) injection. Fuel 227, 218–246.

Aker, E., Jørgen MÅløy, K., Hansen, A., Batrouni, G.G., 1998. A two-dimensional network simulator for two-phase flow in porous media. Transp. Porous Media 32, 163–186.

Alvarado, V., Manrique, E., 2010. Enhanced oil recovery: an update review. Energies 3 (9), 1529–1575.

Aryana, S.A., Kovscek, A.R., 2012. Experiments and analysis of drainage displacement processes relevant to carbon dioxide injection. Phys. Rev. E 86 (6), 066310.

Auradou, H., Måløy, K.J., Schmittbuhl, J., Hansen, A., Bideau, D., 1999. Competition between correlated buoyancy and uncorrelated capillary effects during drainage. Phys. Rev. E 60 (6), 7224.

Aursjø, O., Erpelding, M., Tallakstad, K.T., Flekkøy, E.G., Hansen, A., Måløy, K.J., 2014. Film flow dominated simultaneous flow of two viscous incompressible fluids through a porous medium. Front. Phys. 2, 63.

Ayaz, M., Toussaint, R., Schäfer, G., Måløy, K.J., 2020. Gravitational and finite-size effects on pressure saturation curves during drainage. Water Resour. Res. 56 (10), e2019WR026279.

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

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

Bryant, S., Johnson, A., 2003. Wetting phase connectivity and irreducible saturation in simple granular media. J. Colloid Interface Sci. 263 (2), 572–579.

Buckingham, E., 1914. On physically similar systems; illustrations of the use of dimensional equations. Phys. Rev. 4 (4), 345.

Chen, C., Duru, P., Joseph, P., Geoffroy, S., Prat, M., 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, S., Guo, B., 2023. Pore-scale modeling of PFAS transport in water-unsaturated porous media: Air-water interfacial adsorption and mass-transfer processes in thin water films. Water Resour. Res. 59 (8), e2023WR034664.

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.

Datta, S.S., Weitz, D.A., 2013. Drainage in a model stratified porous medium. Europhys. Lett. 101 (1), 14002.

Dong, H., Blunt, M.J., 2009. Pore-network extraction from micro-computerized-tomography images. Phys. Rev. E—Stat. Nonlinear, Soft Matter Phys. 80 (3), 036307.

Dong, Y., Wang, L., Wei, C., 2022. Characterization of liquid-vapor interfaces in pores during evaporation. Water Resour. Res. 58 (9), e2021WR031908.

Freund, J., Stenberg, R., 1995. On weakly imposed boundary conditions for second order problems. In: Proceedings of the Ninth Int. Conf. Finite Elements in Fluids. Venice, pp. 327–336.

Garfi, G., John, C.M., Rücker, M., Lin, Q., Spurin, C., Berg, S., Krevor, S., 2022. Determination of the spatial distribution of wetting in the pore networks of rocks. J. Colloid Interface Sci. 613, 786–795.

Gjerde, I., Scott, L., 2022. Nitsche's method for Navier-Stokes equations with slip boundary conditions. Math. Comp. 91 (334), 597–622.

Golmohammadi, S., Ding, Y., Kuechler, M., Reuter, D., Schlueter, S., Amro, M., Geistlinger, H., 2021. Impact of wettability and gravity on fluid displacement and trapping in representative 2D micromodels of porous media (2D sand analogs). Water Resour. Res. 57 (10), e2021WR029908.

Hang, S., 2015. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Software 41 (2), 11.

Homsy, G.M., 1987. Viscous fingering in porous media. Annu. Rev. Fluid Mech. 19 (1), 271–311.

Hoogland, F., Lehmann, P., Mokso, R., Or, D., 2016a. Drainage mechanisms in porous media: From piston-like invasion to formation of corner flow networks. Water Resour. Res. 52 (11), 8413–8436.

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

Hoogland, F., Lehmann, P., Or, D., 2017. Capillary flows across layers and textural interfaces-Pathways and colloid transport considerations in unsaturated layered porous media. J. Colloid Interface Sci. 504, 294–304.

Hoppe, H., DeRose, T., Duchamp, T., McDonald, J., Stuetzle, W., 1993. Mesh optimization. In: Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques. pp. 19–26.

Irannezhad, A., Primkulov, B.K., Juanes, R., Zhao, B., 2023a. Characteristics of fluid-fluid displacement in model mixed-wet porous media: patterns, pressures and scalings. J. Fluid Mech. 967, A27.

Irannezhad, A., Primkulov, B.K., Juanes, R., Zhao, B., 2023b. Fluid-fluid displacement in mixed-wet porous media. Phys. Rev. Fluids 8 (1), L012301.

Joekar-Niasar, V., Hassanizadeh, S., 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.

Juanes, R., MacMinn, C.W., Szulczewski, M.L., 2010. The footprint of the CO₂ plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. Transp. Porous Media 82, 19–30.

Kharaghani, A., Mahmood, H.T., Wang, Y., Tsotsas, E., 2021. Three-dimensional visualization and modeling of capillary liquid rings observed during drying of dense particle packings. Int. J. Heat Mass Transfer 177, 121505.

Lenormand, R., 1990. Liquids in porous media. J. Phys.: Condens. Matter. 2 (S), SA79.

Lenormand, R., Zarcone, C., 1989. Capillary fingering: percolation and fractal dimension. Transp. Porous Media 4, 599–612.

Liu, Y., Berg, S., Ju, Y., Wei, W., Kou, J., Cai, J., 2022. Systematic investigation of corner flow impact in forced imbibition. Water Resour. Res. 58 (10), e2022WR032402.

Logg, A., Mardal, K.-A., Wells, G., 2012. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS book, vol. 84, Springer Science & Business Media.

Løvoll, G., Jankov, M., Måløy, K., Toussaint, R., Schmittbuhl, J., Schäfer, G., Méheust, Y., 2011. Influence of viscous fingering on dynamic saturation-pressure curves in porous media. Transp. Porous Media 86, 305–324.
Løvoll, G., Méheust, Y., Måløy, K.J., Aker, E., Schmittbuhl, J., 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.

Løvoll, G., Méheust, Y., Toussaint, R., Schmittbuhl, J., Måløy, K.J., 2004. Growth activity during fingering in a porous Hele-Shaw cell. Phys. Rev. E 70 (2), 026301.

Måløy, K.J., Feder, J., Jøssang, T., 1985. Viscous fingering fractals in porous media. Phys. Rev. Lett. 55 (24), 2688.

Måløy, K.J., Furuberg, L., Feder, J., Jøssang, T., 1992. Dynamics of slow drainage in porous media. Phys. Rev. Lett. 68 (14), 2161.

Måløy, K.J., Moura, M., Hansen, A., Flekkøy, E.G., Toussaint, R., 2021. Burst dynamics, upscaling and dissipation of slow drainage in porous media. Front. Phys. 718.

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

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

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

Moura, M., Flekkøy, E.G., Måløy, K.J., Schäfer, G., Toussaint, R., 2019. Connectivity enhancement due to film flow in porous media. Phys. Rev. Fluids 4 (9), 094102.

Muntoni, A., Cignoni, P., 2021. PyMeshLab. Zenodo, http://dx.doi.org/10.5281/zenodo.4438750.

Nitsche, J., 1971. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. In: Abhandlungen aus dem mathematischen Seminar der Universität Hamburg, vol. 36, Springer, pp. 9–15.

Or, D., Assouline, S., 2013. The foam drainage equation for unsaturated flow in porous media. Water Resour. Res. 49 (10), 6258–6265.

Pavuluri, S., Holtzman, R., Kazeem, L., Mohammed, M., Seers, T.D., Rabbani, H.S., 2023. Interplay of viscosity and wettability controls fluid displacement in porous media. Phys. Rev. Fluids 8 (9), 094002.

Primkulov, B.K., Pahlavan, A.A., Fu, X., Zhao, B., MacMinn, C.W., Juanes, R., 2021. Wettability and Lenormand's diagram. J. Fluid Mech. 923, A34.

Primkulov, B.K., Zhao, B., MacMinn, C.W., Juanes, R., 2022. Avalanches in strong imbibition. Commun. Phys. 5 (1), 52.

Reis, P., Moura, M., Linga, G., Rikvold, P.A., Toussaint, R., Flekkøy, E.G., Måløy, K.J., 2023. A simplified pore-scale model for slow drainage including film-flow effects. Adv. Water Resour. 182, 104580.

Scheel, M., Seemann, R., Brinkmann, M., Di Michiel, M., 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.

Suo, S., 2024. Spontaneous imbibition of a liquid film wetting a wall-mounted cylinder corner. Soft Matter 20 (3), 578–587.

Toussaint, R., Måløy, K.J., Méheust, Y., Løvoll, G., Jankov, M., Schäfer, G., Schmittbuhl, J., 2012. Two-phase flow: Structure, upscaling, and consequences for macroscopic transport properties. Vadose Zone J. 11 (3), vzj2011–0123.

Urquiza, J.M., Garon, A., Farinas, M.-I., 2014. Weak imposition of the slip boundary condition on curved boundaries for Stokes flow. J. Comput. Phys. 256, 748–767.

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

Vorhauer, N., Wang, Y., Kharaghani, A., Tsotsas, E., Prat, M., 2015. Drying with formation of capillary rings in a model porous medium. Transp. Porous Media 110, 197–223.

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

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

Zacharoudiou, I., Boek, E.S., Crawshaw, J., 2018. The impact of drainage displacement patterns and Haines jumps on CO₂ storage efficiency. Sci. Rep. 8 (1), 15561.

Zhang, W., Guo, B., 2024. Anomalous adsorption of PFAS at the thin-water-film air-water interface in water-unsaturated porous media. Water Resour. Res. 60 (3), e2023WR035775.

Zhao, B., MacMinn, C.W., Juanes, R., 2016. Wettability control on multiphase flow in patterned microfluidics. Proc. Natl. Acad. Sci. 113 (37), 10251–10256.

Zhao, J., Qin, F., Kang, Q., Qin, C., Derome, D., Carmeliet, J., 2022. A dynamic pore network model for imbibition simulation considering corner film flow. Water Resour. Res. 58 (7), e2022WR032332.