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

HAL
open science

# Gravity stabilized drainage in porous media with controlled disorder

Khobaib Khobaib, Paula Reis, Marcel Moura, Renaud Toussaint, Eirik Grude

Flekkøy, Knut Jørgen Måløy

► To cite this version:

Khobaib Khobaib, Paula Reis, Marcel Moura, Renaud Toussaint, Eirik Grude Flekkøy, et al.. Gravity stabilized drainage in porous media with controlled disorder. Physical Review Research, 2025, 7, 10.1103/PhysRevResearch.7.023040. insu-05082244

HAL Id: insu-05082244

https://insu.hal.science/insu-05082244v1

Submitted on 16 Aug 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Gravity stabilized drainage in porous media with controlled disorder

Khobaib Khobaib<sup>1,*</sup> Paula Reis<sup>1,†</sup> Marcel Moura<sup>1</sup> Renaud Toussaint<sup>1,2</sup>

Eirik Grude Flekkøy<sup>1,3</sup> and Knut Jørgen Måløy<sup>1,4,‡</sup>

<sup>1</sup>PoreLab, The Njord Centre, Department of Physics, University of Oslo, P. O. Box 1048 Blindern, 0316 Oslo, Norway

<sup>2</sup>Université de Strasbourg, CNRS, ENGEES, Institut Terre et Environnement de Strasbourg, UMR 7063 Strasbourg, France

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

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

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

(Received 20 September 2024; accepted 26 March 2025; published 11 April 2025)

This work presents an experimental, numerical, and theoretical investigation of slow drainage processes in porous media, focusing on systematic variations in pore-scale disorder. We examine a system comprising a monolayer of cylinders in a quasi-two-dimensional (2D) geometry, where the medium’s disorder is controlled via a disorder parameter, $\epsilon$. By systematically altering $\epsilon$, we randomly shift the positions of the cylinders and analyze the resulting effects on the gravity stabilized invasion front. Our results show that increasing $\epsilon$ has a significant impact on the stable width of the invasion front, with this width scaling with the disorder parameter according to an exponent $\beta = 0.57$, which closely matches the theoretical value predicted by percolation theory for 2D systems. Additionally, we derive a reciprocal relationship between the dimensionless fluctuation number and the disorder parameter $\epsilon$. Moreover, we find a correlation between the invasion front width and the size of trapped clusters in the system, and we explore how pore disorder affects the air saturation left behind the invasion front. Both our experimental and numerical findings exhibit exponent values that align with theoretical scaling predictions. This study advances the current understanding of multiphase flows in porous media by specifically addressing the influence of pore-scale disorder on the morphology of invasion fronts.

DOI: 10.1103/PhysRevResearch.7.023040

# I. INTRODUCTION

Liquid flow through porous media is highly relevant in various scientific and engineering domains. Porous media, which are found in both natural and industrial processes, cover a wide range of materials, including geological formations and manufactured goods such as filters and membranes. Understanding the dynamics of fluid flow within these complex networks of interconnected void spaces is essential for optimizing processes including fuel cells [1], solar technology [2], oil and gas engineering [3–8], subsurface carbon storage [9], water filtration into the soil [10], and chemical filtration [11]. For many decades, this multidisciplinary subject has drawn the attention of researchers from diverse scientific backgrounds such as physics, chemistry, biology, geoscience, hydrology, and engineering.

When a fluid displaces another immiscible fluid in a porous medium, the resulting displacement patterns exhibit a diverse range of shapes [12–18], from compact to complex

fractal structures [19,20]. These patterns primarily depend on the relevant forces, such as the viscous [12–14,21–25], capillary [14,22–30], and gravitational forces [25,30–36], the wetting properties [15,37–42], as well as the local geometry of the porous medium [43–45]. Among these, a crucial factor affecting fluid displacement patterns is the disorder of the porous medium. Extensive research has been conducted to investigate the impact of flow conditions, wettability, and gravity on fluid displacement in porous media. However, comparatively less attention has been given to studying and quantifying the effect of pore-scale disorder.

In drainage experiments at sufficiently small wetting angles $\theta$, the displacement patterns in porous materials are directly influenced by the movement of the menisci between the wetting and nonwetting fluids at the pore scale. Here $\theta$ is zero when the wetting fluid perfectly wets the solid matrix. In disordered media, pore-throat widths are not equal throughout the porous material. The nonwetting fluid can only penetrate a pore throat if the capillary pressure, $p$ (the pressure difference between the nonwetting and wetting fluids), exceeds the capillary threshold value, $p_t$, of that throat. This threshold pressure is directly related to the geometry of the pore throat. According to the Young-Laplace law [46], smaller pore throats correspond to higher capillary pressure thresholds, making it more difficult for the nonwetting fluid to invade them. The effect of pore-size disorder alone will result in an unstable invasion front [31,32]. An imposed gravitational or viscous field, on the other hand, can either stabilize or destabilize the front, introducing a characteristic length scale

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*Contact author: khobaibk@fys.uio.no</footnote>

<footnote>†Contact author: paula.reis@mn.uio.no</footnote>

<footnote>‡Contact author: maloy@fys.uio.no</footnote>
into the system with a fractal structure from the size of the pores to this length scale. In the case of a stabilizing field, for example when a low-density fluid enters a high-density fluid from the top in a gravitational field, this length scale dictates the width of the front and the size of the trapped clusters behind it [30,32]. In the case of a destabilizing field, such as viscous fingering, with a low-viscosity fluid invading a more viscous one, or gravitational fingering, with a low-density fluid invading a high-density fluid from the bottom, the length scale determines the characteristic width of the fingers [22,30]. The structure will thus typically be fractal up to this characteristic length scale, transitioning to a new fractal structure in the unstable scenario or to a flat front when the field stabilizes [22,32,47,48].

For large wetting angles, the interface traveling through nearby pores may overlap or touch [49–51], resulting in correlated motion between neighboring pores and diminishing the influence of the disorder on the instability [49–51]. This effect has been investigated systematically both in experiments [52,53] and simulations [49–52] when the viscous field destabilizes the invasion front, giving rise to smoother and thicker viscous fingers. At low injection rates where capillary forces dominate, compact structures are observed [52,53]. The idea of using percolation theory to study two-phase flow in porous media has been explored in the 1980s [54,55], and its theoretical development was significantly supported by the invasion percolation model by Wilkinson [56]. Since then, numerous experimental, numerical, and theoretical studies have been carried out on this topic [31,32,57–60]. The dimensionless fluctuation number $F = Ga/W_p$ [30,35,60], where $W_p$ is the width of the normalized capillary pressure threshold distribution, *a* is the characteristic pore-throat size, and *G* is an external pressure gradient field, gives the scaling of the width of the invasion front $\eta \propto F^{-\nu/(1+\nu)}$ [30–32,60]. Here $\nu = 4/3$ is the critical correlation length exponent in two-dimensional (2D) percolation. In our case $G = \Delta\rho g_0 \sin \alpha$, where $g_0$ is the acceleration of gravity and $\alpha$ is the tilt angle of the model with respect to the horizontal level. The fluctuation number *F* can be compared to other dimensionless numbers often used to classify flow patterns in porous materials, such as the Bond number $Bo$ which compares gravitational and capillary forces. However, unlike the fluctuation number *F*, the Bond number does not provide *a priori* information on the actual material disorder. Earlier investigations have focused on changing *F* or $Bo$ by changing the external field $G$, which is considerably easier experimentally than changing the disorder. Vincent-Dospital *et al.* [30] introduced a gradient in disorder and, consequently, in the average capillary pressure.

In this work, we focus on the case of stable drainage under the influence of a gravitational field, where a low-density fluid (air) displaces a high-density fluid (water-glycerol solution) from above. The system is characterized by a sufficiently small wetting angle to prevent the interface between neighboring pores from touching or overlapping. The main goal of this project is to systematically investigate the impact of disorder, i.e., of capillary pressure distribution width, while keeping all other parameters constant. We thus introduce experiments and simulations where we systematically change the disorder in a controlled way while keeping the grav-

<figure><img src="image_4.png" /><figcaption>FIG. 1. A sketch illustrating the pore-scale disorder in a 3D-printed model, where *r* is the radius of an arbitrary pore throat. The diameter of the cylinders is Ø and *a* denotes the separation between the cylinders in the initial triangular array (zero disorder). The dashed circles indicate the cylinder’s position at zero disorder conditions, while the solid black circles represent the cylinder positions when the disorder is introduced into the model.</figcaption></figure>

itational field and fluid parameters such as wetting angle, density, and surface tension fixed. The quantitative control of the disorder allows for an analytical prediction of the front width that compares well with our experimental and simulated results.

## II. MATERIALS AND METHODS

Our experimental and numerical methods are described in this section. In both cases, we studied quasi-2D porous networks formed by a monolayer of cylinders, a modification of the classic Hele-Shaw geometry [61] widely used in fluid mechanics studies. We first consider a setup in which the cylinders are arranged in a crystalline triangular assembly. We then introduce disorder into the network systematically by means of a parameter $\epsilon$ theoretically ranging from 0 to 1. The ordered crystalline case corresponds to $\epsilon = 0$ and $\epsilon = 1$ is the system with maximum disorder. A more detailed description of this parameter will be given in the coming section, but in essence the method involves moving each cylinder in a random direction by a random distance $\zeta$ uniformly selected in the interval $0 \le \zeta \le \zeta_{max}$, with $\zeta_{max} = \epsilon a/2$ where *a* is the initial separation between the cylinders in the crystalline array. In Fig. 1, we show a schematic representation of the ordered (dashed circles) and disordered (filled circles). For the experimental case, we limited this parameter to the range $0.1 \le \epsilon \le 0.9$ in order to avoid possible artifacts in the extreme cases. A truly crystalline case $\epsilon = 0$ cannot be created due to the inherent finite resolution of the 3D-printing technique employed (the surface fluctuations coming from the printing itself would dominate the dynamics in a case without disorder). The other extreme case with $\epsilon = 1$ is also avoided as the disorder here is such that neighboring cylinders could touch, thus changing the pore connectivity.
<figure><img src="image_5.png" /><figcaption>FIG. 2. A sketch illustrating the experimental setup: (a) A 3D-printed model is tightly fastened between an air-filled pressure cushion and a plexiglass plate. The top-sealed air cushion ensures a leakage-free environment. The model is placed on the illumination lamp; the entire setup is carefully tilted at an angle $\alpha = 1.5^{\circ}$. Capturing images from the top, the experimental design affords precise control over the gravitational field. Consequently, it enables a detailed examination of the capillary fluctuations at the pore level under a hydrostatic pressure gradient. (b) A top-view image where a nonwetting fluid (air) invades a wetting fluid (water-glycerol solution) in a porous model with a coordinate system $(x, y)$, and the dimensions of the model are 140 $\times$ 120 $\times$ 4.0 mm, where the cylinder's height is 2 mm and the plate's thickness under the cylinder is 2 mm. The width of the invasion front is denoted by $2\eta$, and the yellow dotted line indicates the average position of the invasion front. $\Lambda$ represents the maximum length of a trapped cluster.</figcaption></figure>

### A. Experimental method

A sketch of the experimental setup is shown in Fig. 2. The study involved conducting experiments on 3D-printed porous models including varying degrees of pore disorder. Each model was produced using Form 3B, a 3D printer from Form-labs, which employs stereolithography (SLA) technology. This process utilizes a laser technique to systematically cure clear liquid resin layer by layer, resulting in the creation of highly accurate and transparent 3D models, with spatial resolution of 50 µm. This technique enables the meticulous control of the porous network's geometry, offering specific tunings for the distribution of grain and pore-throat sizes within a printed model. The cylinders, with height $h = 2.0$ mm and diameter $\Ø = 1.0$ mm, are arranged in a monolayer on a rectangular plate having dimensions (140 $\times$ 120 mm). The initial pore-throat size of a model without disorder is $a = 1.0$ mm. The porosity of the models is $\phi \approx 73\%$. Following the printing process, each model undergoes a thorough cleaning in an ultrasonic bath containing isopropanol for a period of 40 min. Later, it is cured in a UV chamber for 50 min to enhance its structural integrity, reducing the vulnerability to any form of damage during the experiments. The isopropanol cleaning and UV curing are standard procedures in stereolithography printing. The final setup consists of a 3D-printed porous model, a plexiglass frame, and an air-filled pressure cushion. The air cushion seals the model from the top and provides uniform pressure across its surface. The printed porous model is attached with double-sided tape on 4-mm-thick plexiglass frames, which have a cylindrical inlet and outlet leading to broad channels spanning the entire model's width. This design aims to reduce boundary effects during liquid drainage. The cylindrical inlet and outlet allow the injection of wetting and nonwetting fluids throughout the experiment. The sides of each model are sealed with a gasket to avoid any leakage. Each

model is then placed on an LED lightbox setup and the whole setup is tilted above the horizontal at an angle $\alpha = 1.5^{\circ}$, thus allowing a small stabilizing gravitational component $g_0 \sin \alpha$ along the flow direction.

The models are initially filled with the wetting fluid composed of a 80% glycerol and 20% water solution mixture (the percentages are given by weight). In this mixture, Nigrosin dye (4g/l of water) is added to make the wetting phase dark, thus allowing for easier visualization and later segmentation via image analysis. The contact angle of the wetting phase with the cured resin and air cushion film is $60 \pm 10^{\circ}$ in the presence of the nonwetting phase (air). The dynamic viscosity of the water-glycerol mixture, density, and surface tension of the wetting phase with respect to air are $\mu = 44.45$ mPa s, $\rho = 1198$ kg m$^{-3}$ and $\gamma = 6.4 \times 10^{-2}$ N m$^{-1}$, respectively. We assume that the dynamic viscosity and density of the non-wetting phase (air) are negligible. The density difference $\Delta\rho$ between wetting and nonwetting fluid is then $\Delta\rho \approx \rho$. During the experiment, the wetting fluid is withdrawn at a constant flow rate $q = 1.00$ ml h$^{-1}$ using a syringe pump, see Fig. 2. Due to the slow drainage, viscous forces are negligible, resulting in a relatively small capillary number ($Ca = \mu V/\gamma) \approx 8 \times 10^{-7}$, where V represents the average flow velocity. The experiments are terminated once the nonwetting fluid reaches the outlet channel of the model.

#### 1. Designing of pore-scale disorder in 3D-printed porous media: A sketch

Figure 1 represents a sketch of the pore-scale disorder in a 3D-printed model, providing insights into how the disorder affects the arrangement of cylindrical structures within the porous material. In our system, we originally placed the cylinders on a triangular lattice and then introduced disorder into the model. The disorder causes the cylinders to move a
<figure><img src="image_6.png" /><figcaption>FIG. 3. The standard deviation of pore-throat diameter normalized with $2r$ plotted against the $W_r/(2r)$. The data points from left to right correspond to $\epsilon$ values ranging from 0.1 to 0.9.</figcaption></figure>

distance $\zeta$ in a random direction, where $\zeta = \epsilon ja/2$. Here $j$ is a random number chosen from a flat distribution between 0 and 1, and $a$ is the initial pore-throat size in the triangular array; see Fig. 1. In our 3D-printed models, we introduced disorder within the parameter range of $0.1 \le \epsilon \le 0.9$, using an initial pore-throat size of $a = 1.0$ mm. When $\epsilon > 0.9$, we noticed that the cylinders in the model start to touch each other, indicating a critical disorder point for printing. This is due to potential fluctuations in the cylinder's surfaces coming from the 3D-printing technique. This experimental limitation is, of course, not present in simulations, where we chose the disorder parameter in the interval $0.1 \le \epsilon \le 1.0$. In Fig. 1, the dashed circles show the original position of the cylinders under disorder-free conditions. In contrast, the solid black circles show the randomly shifted positions of the cylinders when the disorder is introduced into the model. Consider an arbitrary pore throat with a radius $r$ defined as the half distance between the neighboring cylinders that define the pore throat (see Fig. 1). Thus, the maximum range of the distribution $W_r$ of pore-throat radii $r$ is

$$W_r = [a(1 + \epsilon) - a(1 - \epsilon)]/2 = \epsilon a. \quad (1)$$

In Fig. 3, the normalized standard deviation of the pore-throat diameter is plotted against the width of the distribution $[\text{std}(2r)/2r$ vs $W_r/2r]$. The result shows a linear relationship with a slope of approximately 1.0, which validates that $W_r$ is proportional to standard deviation of the pore-throat diameter. Consequently, both $W_r = \epsilon a$ and $\text{std}(2r)$ can serve as characteristic measures for the width of the pore-throat distribution. Through systematic parameter variations of $\epsilon$, we explore the resulting effects on the overall drainage behavior in 3D-printed porous model.

## 2. Analyzing capillary pressure thresholds and pore structures in 3D-printed porous models

We assess the distribution of pore-throat sizes by analyzing initial images corresponding to each $\epsilon$ value obtained from experiments, when the porous network is still fully saturated

with a water-glycerol solution [see Figs. 4(a)–4(c)]. Initially, we identify the position and average radius of each cylinder in the network, marked by yellow dots. Subsequently, to determine distances between consecutive cylinder centers, we employ the Delaunay triangulation algorithm [62], shown by red lines in Figs. 4(d)–4(f). Next, we produce the Voronoi lattice shown in green lines. The vertices of the Voronoi lattice correspond to the pore centers, while the intersections between the Delaunay triangulation lines and the Voronoi lattice indicate pore-throat centers. By subtracting cylinder radii from the distances between their centers, we directly determine the distribution of pore-throat sizes. In Figs. 4(g)–4(i) we show the probability density function (PDF) associated with the pore-throat size distributions. We see that the distribution of pore-throat sizes becomes wider with $\epsilon$, as expected from Eq. (1). The actual distribution of capillary pressure thresholds ($p_t$) was determined using the Young-Laplace equation as follows:

$$p_t(r) = \gamma \cos \theta \left( \frac{1}{r} + \frac{2}{h} \right), \quad (2)$$

where $\theta$ is the contact angle between the air-liquid interface and the solid surface, $r$ is half the distance between adjacent cylinders in the 3D-printed model, and $h$ is the height of the cylinders.

We consider the cumulative distribution function (CDF), denoted as $f(p)$, which gives the probability that the random variable $p_t$ is less than or equal to a specific value $p$. This quantity is represented by the integral of the PDF, denoted as $N(p_t)$ (normalized capillary pressure threshold distribution, which integrates to 1, over the range of capillary pressures), from negative infinity up to $p$, i.e.,

$$f(p) = \int_{-\infty}^{p} N(p_t) dp_t. \quad (3)$$

Figures 4(j)–4(l) shows the cumulative distribution function associated with the capillary pressure threshold $p_t$ [calculated from Eq. (2)].

## B. Numerical method

As a complement to the experiments, the effects of pore-size disorder in drainage were also investigated numerically, using a 2D bond-invasion-percolation model with trapping [54–56]. Such models conceptualize the pore space as networks of pores (sites) connected by pore throats (bonds) and represent slow drainage as a sequence of discrete steps, in each of which one bond is invaded.

In the case of drainage under the influence of gravity, each bond of the network is associated with an invasion threshold $p_t$, related to the pore throats' sizes, and a local capillary pressure $p$, related to the pore throats' position [32,35,63,64]. The difference between these two values defines the bond's likelihood of invasion. In each step of our simulations, the bonds belonging to the interface of the main defending cluster—the wetting-phase cluster connected to the network's outlet—are ranked based on $p - p_t$, and the one with the highest value is invaded.

The following steps were implemented to generate networks similar to the 3D-printed porous media used in the
<figure><img src="image_7.png" /><figcaption>FIG. 4. [(a)–(c)] Top view of a small section of 3D-printed porous model filled with a water-glycerol solution for disorder levels of $\epsilon = 0.1, 0.6$, and $0.9$. [(d)–(f)] Image analysis is employed in experimental investigations to determine the distribution of pore-throat sizes for three different $\epsilon$. [(g)–(i)] Probability density functions (PDFs) of normalized pore-throat sizes for the pore geometries with the respective disorder levels. [(j)–(l)] Cumulative distribution functions (CDFs) of the capillary pressure $p$. The porosity ($\phi$) and average pore-throat size ($\overline{2r}$) remain consistent across all levels of disorder in the 3D-printed models.</figcaption></figure>

experiments: Instead of using the actual 3D-printed model, based on a similar geometry, a triangular lattice was employed as a framework for the cylinder positions without disorder. Then, these positions were shifted according to the disorder parameter $\epsilon$, as described in Sec. II A1 and seen in Fig. 1. Using the final position of the cylinders as a set of points, a Voronoi diagram was created, in which the cells' faces correspond to the pore throats, and the cell's vertices correspond to the pores. In this way, the resulting pore networks describe irregular hexagonal lattices. Finally, the values of $p_t$ are calculated using Eq. (2), and the local capillary pressure, including the contribution from gravity, is given by:

$$p(x) = p_0 + \Delta \rho g_0 x \sin \alpha, \quad (4)$$

where $p_0$ is a reference capillary pressure at the model's outlet, at $x = 0$, $x$ is a linear coordinate pointing upwards tangent to the cell and $\alpha$ is the tilt angle of the model with respect to the horizontal plane, as seen in Fig. 2.

The employed invasion percolation (IP) model [54–56] was implemented in the same way described in Ref. [64], except for the representation of film flow effects [65–67].

Given the geometrical aspects and wettability of the 3D-printed models, we assumed that flow through films behind the main invasion front would be negligible.

The numerical results presented in Sec. IV were obtained using 50 different random realizations of the pore networks for each level of disorder $\epsilon$. In this way, we could assert that the observations regarding front width, cluster size, and saturation are consistently related to the porous media disorder, not a particular porous medium model. Moreover, the simulated porous medium domain was equivalent to that used in the 3D-printed models (14 × 12 cm) and the same fluid and wettability properties presented in Sec. II A were used in Eqs. (2) and (4). Consequently, numerical and experimental results could be directly compared, as shown in Fig. 5. The images represent the porous medium fluid occupation at two-thirds of the total drainage time, with the main invasion front highlighted in red. Cases corresponding to $\epsilon = [0.1, 0.5, 0.9]$ are displayed to illustrate how the effect of disorder on the front with, cluster sizes and saturation are similarly captured by both the experimental and numerical methods adopted in this study.
<figure><img src="image_8.png" /><figcaption>FIG. 5. [(a)–(c)] The experimental images for disorder levels of $\epsilon = 0.1, 0.5$, and $0.9$. Similarly, (d)–(f) show the simulated images for the same disorder levels. In these images, the black color represents the wetting phase (water-glycerol solution), while the light color indicates the nonwetting phase (air). The cylinders are depicted in a light color but can be identified by their distinct circumference. The red line indicates the drainage front and the black arrow represents the flow direction. See also the corresponding movies S1 and S2 in the Supplementary Material [68].</figcaption></figure>

### III. THEORY: SCALING BEHAVIOR OF THE INVASION FRONT

The experiment begins with the invasion of the nonwetting fluid by drainage of the wetting fluid which has a higher density and viscosity. By carefully withdrawing the wetting fluid at a controlled, slow rate, we ensure that capillary forces govern the drainage process entirely, and in the meantime, gravity stabilizes the fluid front. Each pore throat has a capillary pressure threshold $p_t$ that must be overcome in order to invade that pore throat. The relationship between the occupation probability in percolation theory denoted as $f$, and the capillary pressure, denoted as $p$, is then given by the cumulative distribution function $f(p)$, Eq. (3) [25,69,70].

The difference between the occupation probability $f$ and the critical occupation probability $f_c$ can be written as

$$f - f_c = \int_{p_c}^{p} N(p_t) dp_t, \qquad (5)$$

where $p_c$ is the critical percolation pressure [69,71].

Now, performing a Taylor expansion of $N(p_t)$ around $p_c$ in Eq. (5) and considering only the lowest order term in $p - p_c$, we obtain

$$f - f_c = N(p_c)(p - p_c). \qquad (6)$$

There are no trapped clusters of nonwetting fluid in the displaced fluid (wetting fluid); therefore, it is realistic to assume

that the viscous pressure drop in the wetting fluid will depend linearly on the length scale. Let us consider the capillary pressure at a height $x$. Since the gravitational field also exhibits a linear dependence on the length scale $x$ (see Eq. (4)), we can write

$$f - f_c = N(p_c) \left( \Delta \rho g_0 a \sin \alpha - \frac{q \mu a}{\kappa A} \right) \frac{(x - x_0)}{a}, \qquad (7)$$

where $x_0$ represents the height along the front where the capillary pressure reaches the percolation threshold, $\kappa$ represents the permeability of the model and $A$ is the cross-sectional area of the porous medium. Equation (7) can be written as

$$f - f_c = F \frac{(x - x_0)}{a}, \qquad (8)$$

where $F$ is the dimensionless number called fluctuation number which governs the invasion process [25,60]

$$F = N(p_c) \left( \Delta \rho g_0 a \sin \alpha - \frac{q \mu a}{\kappa A} \right), \qquad (9)$$

the quantity $1/N(p_c)$ represents the typical width of the normalized capillary pressure threshold distribution. The dimensionless number $F$ depends on the capillary pressure and gravitation field as well as the viscous field.

Now, we will make use of an assumption introduced by Sapoval [25,32,56,72,73] for the correlation length in percolation, denoted as $\xi$, which scales with the width of the front denoted as $\eta$ (standard deviation of the front along the
position x), $\xi \propto (x - x_0) = \eta$. This assumption is based on the observation that the largest trapped cluster of wetting fluid is controlled by the invasion front width. According to percolation theory, the correlation length $\xi$ is given by [69]

$$ \frac{\xi}{a} \propto |f - f_c|^{-\nu}, \qquad (10) $$

with $\nu$ denoting the critical exponent and equal to $\frac{4}{3}$.

By using Sapoval's argument [73] and using Eqs. (10) and (8) we get [25,32,56,72,73]

$$ \frac{\eta}{a} \propto F^{\frac{\nu}{1+\nu}}, \qquad (11) $$

where the critical exponent $\nu$ is dependent on the dimensionality of the system. For 2D $\nu = 4/3$ and for 3D $\nu \approx 0.88$ [69]. Now we calculate the variation of the capillary pressure threshold with the corresponding pore disorder from Eq. (2). Using a first-order Taylor expansion of $p_t(a/2 + \delta r)$ in $\delta r = r - a/2$ around the initial pore radius $a/2$, the change in capillary pressure $\delta p_t = p_t(a/2) - p_t(r)$ for a small variation in pore-throat size $\delta$ is

$$ \delta p_t = \gamma \cos \theta \frac{4\delta r}{a^2}. \qquad (12) $$

Due to the honeycomb lattice structure in our model [Fig. 4(d)–4(f)], the critical bond percolation threshold is $f(p = p_c) = 0.6527$ [71]. The critical pore size $r_c$ corresponds to the critical pressure threshold $p_c$, determined from CDF ($f(p)$), see Figs. 4(j)–4(l). The normalized distribution at a critical pore length $r = r_c$ denoted as $N(r_c)$ is inversely related to the characteristic width of the distribution $W_r$. Thus, we express the relationship as

$$ N(r_c) \approx \frac{1}{W_r}. \qquad (13) $$

Similarly, the normalized distribution at the critical percolation pressure $p_c$ denoted as $N(p_c)$ is also inversely related to the characteristic width of the pressure distribution $W_p$. Thus, we can express this relationship as

$$ N(p_c) \approx \frac{1}{W_p}. \qquad (14) $$

Then identifying $W_r$ with $\delta r$ and $W_p$ with $\delta p_t$ and using Eqs. (12) and (1) we get

$$ W_p = W_r \frac{4\gamma \cos \theta}{a^2} = \frac{4\epsilon\gamma \cos \theta}{a}. \qquad (15) $$

In our experiment, the wetting fluid is withdrawn at a sluggish flow rate, resulting in a negligible viscous pressure ($\approx 10^{-5} \text{N m}^{-2}$). Therefore, the F number is only determined by the normalized capillary pressure threshold distribution at $p_c$ and the hydrostatic pressure gradient. Equation (9) can be simplified as

$$ F = N(p_c) \Delta \rho g_0 a \sin \alpha, \qquad (16) $$

using Eqs. (11), (14), (15), and (16), we can write

$$ \frac{\eta}{a} \propto \left( \frac{a^2 \Delta \rho g_0 \sin \alpha}{\epsilon \gamma \cos \theta} \right)^{-\beta} \propto \epsilon^{\beta}, \qquad (17) $$

where the exponent $\beta = \frac{\nu}{1+\nu}$ is equal to 0.57.

<figure><img src="image_9.png" /><figcaption>FIG. 6. Front width $\eta$ plotted against time for different levels of pore disorder parameters $\epsilon$. The black dotted vertical lines denote the interval used for evaluating the mean width.</figcaption></figure>

## IV. RESULTS AND DISCUSSION

### A. Effect of pore disorder on the width of the invasion front

In this section, we experimentally and numerically investigate the impact of pore disorder on the width of the invasion front. For this, the front width $\eta$ is calculated as the standard deviation of the x coordinate of the front pixels from the experimental images (see Fig. 2), or the front network nodes, in the numerical case.

As an illustration of the experimentally measured front width values, in Fig. 6, $\eta$ is plotted against the time for various levels of disorder $\epsilon$. For each $\epsilon$, we observe an initial growth in the front width, which gradually transitions into a fluctuating but stabilized state. This stabilization occurred roughly 150 min into the experiments (this corresponds to $\approx 6\%$ of the pore volume of the samples). In the plots, we can also observe a direct correlation between the magnitude of pore disorder and the increase of the front width. To obtain a single representative front width value for each experiment, we calculated the mean front width from measurements collected between 200 and 800 min. In this region, the front width is relatively stable, and there is no noticeable influence from the boundaries of the model.

A similar method was used to calculate the mean front width for each simulated case. After obtaining the standard deviation of the front nodes' positions at each simulation time step, values corresponding to the beginning and end of drainage were discarded, and the mean of the remaining values was calculated. These mean values of front width, denoted by $\bar{\eta}$, are normalized by mean pore-throat size $2r$ and used for further analysis in Fig. 7.

Figure 7(a), shows the relationship between the degree of pore disorder $\epsilon$ and the normalized mean width of the invasion front ($\bar{\eta}/2r$). To ensure the robustness and reproducibility of our findings, we conducted several experiments for each $\epsilon$, shown by black solid spheres, primarily using newly printed models, but also on some previously used models. As the results were indistinguishable in both scenarios, we could assume that any aging effect due to the absorption
<figure><img src="image_10.png" /><figcaption>FIG. 7. The analysis focused on the normalized mean invasion front width ($\overline{\eta}/2r$) as a function of $\epsilon$ in (a) and as a function of $F$ number in (b). The black solid sphere symbol represents experimental data. Multiple experiments are performed for each value of $\epsilon$ to ensure reproducibility. The simulations consist of 50 different networks for each $\epsilon$, as indicated by the blue sphere symbol. The average across all simulations is represented by the red hexagon symbol. The black solid line signifies the theoretically predicted slope of 0.57 and -0.57.</figcaption></figure>

of water during the experiment or from the atmosphere was negligible.

In addition to experimental data, results from numerical simulations are presented in the figure. The blue spheres symbolize 50 different network simulations for each $\epsilon$, and the red hexagon represents the average across all simulations. Both experimental and numerical results show that the width of the invasion front gradually increased with the increase of pore disorder. As a guide to the eye, a black solid line with a slope of 0.57 is shown. This slope corresponds to the theoretically predicted scaling from Eq. (17).

In Fig. 7(b), we illustrate the relationship between the normalized width of the invasion front and the dimensionless $F$ number, calculated from Eq. (16), where $N(p_c)$ was estimated as follows:

$$N(p_c) = \left. \frac{df(p)}{dp} \right|_{p=p_c} \qquad (18)$$

In our models, the pores are connected by throats in a honeycomb lattice, see Figs. 4(d)–4(f), meaning that each lattice site connects to three neighboring sites. Therefore, the critical bond percolation threshold is equal to $f(p = p_c) = 0.6527$ [71], and the corresponding $p_c$ values can be obtained for each $\epsilon$ from the $f(p)$ curves illustrated in Figs. 4(j)–4(l).

The scaling observed in the Fig. 7(b) confirms the theoretical prediction from Eq. (11). Using Eqs. (11) and (17), we see the reciprocal relationship between fluctuation number and pore disorder

$$F \propto \epsilon^{-1}. \qquad (19)$$

By comparing the experimental data with the simulation results depicted in Fig. 7(a) and 7(b), it becomes evident that the 2D numerical simulations consistently reproduce the observations made using the 3D-printed porous models. In Fig. 7(a), some inconsistency can be observed between the simulated and experimental data, especially at small $\epsilon$ values. Similarly, in Fig. 7(b), the fluctuation numbers differ when comparing experiments to simulations. In this case, the results

show that the magnitude of the $F$ numbers measured from 3D-printed models exceeds what we found in the simulations, for the same $\epsilon$ values. This means that larger pore-throat-size distributions are found in the 3D-printed models in contrast to the simulated networks, indicating the presence of manufacturing artifacts, including induced roughness and nonuniform model thickness. These differences, however, do not significantly impact our analysis of the effects of disorder on the drainage front width, as theoretical predictions can be supported by both numerically and experimentally obtained data, within a small margin of error.

## B. Understanding the influence of pore disorder on trapped cluster length

Previous studies [52,74] have explored the impact of the pore-size disorder on immiscible viscous fluids displacement. They suggest that reducing disorder in the porous medium enhances displacement efficiency and reduces the size of the trapped clusters.

Here the length of the trapped clusters left behind the invasion front is also investigated. In Figs. 8(a) and 8(b), the normalized maximum cluster length in the direction of fluid flow, $(\Lambda/2r)$, is plotted as a function of $\epsilon$ and the $F$ number. The plots reveal a consistent pattern, indicating a correlation between maximum trapped cluster size and both $\epsilon$ and the $F$ number, and in line with the trends reported from previous studies [52,74]. More specifically, the results indicate a slope of 0.57 for the plot of $(\Lambda/2r)$ vs $\epsilon$, and -0.57 when $(\Lambda/2r)$ is plotted against $F$ number, confirming the reciprocity between these values pointed in Eq. (19). Besides that, these scaling coefficients are the same as the ones presented in Fig. 7, for the relationship between front width, $\epsilon$, and $F$ number. This consistency suggests that the width of the invasion front is indeed primarily influenced by the size of trapped clusters, in line with what is expected from Sapoval’s argument [73] and explained in Sec. III.
<figure><img src="image_11.png" /><figcaption>FIG. 8. Plotting the normalized maximum cluster length $\Lambda/2r$ as a function of $\epsilon$ in (a), and $F$ number in (b), which reveals distinctive trends. Similarly to Fig. 7, symbols denote individual experiments and simulations. The presence of a black solid line signifies the predicted slopes, providing a visual reference for comparative analysis.</figcaption></figure>

Similarly to the results presented in Sec. IV A, satisfactory qualitative alignment between experimental findings and simulated data was observed in Fig. 8. The discrepancies between the two data sets, as mentioned in the previous section, can be attributed to artifacts within the printed model.

### C. Influence of pore disorder on air saturation left behind the invasion front

In Secs. IV A and IV B, we established a direct correlation between pore disorder and the widening of the invasion front and the increase of trapped cluster sizes. Now, we investigate the influence of pore disorder on the saturation of the nonwetting fluid behind the invasion front at the end of the drainage.

In previous studies, pressure saturation relationships were examined for a broad range of capillary numbers through drainage experiments conducted on quasi-two-dimensional porous models [48,71,75]. Similarly, Ayaz et al. [63] conducted both experimental and numerical analyses to explore the saturation of a nonwetting fluid during drainage under different degrees of gravitational forces. In that study, a scaling relation between nonwetting fluid saturation and the Bond number was proposed, as

$$S_{nw}^{F} \propto \text{Bo}^{\frac{-\nu(D-d)}{1+\nu}}, \quad (20)$$

where the Bond number, $\text{Bo} = \Delta\rho g_0 a^2 / \gamma$, represents the ratio of the gravitational field to capillary pressure on a single-pore length scale, $D$ is the fractal dimension of the invasion structure with a value of 1.82 for 2D systems [56], and $d$ is the spatial dimension 2.

Therefore, in the quasi-2D models used in our analysis, the proposed exponent $\varphi = \frac{-\nu(D-d)}{1+\nu}$ is equal to 0.10. In a more general description, however, when the disorder in the capillary pressure thresholds is taken into consideration, one must consider the fluctuation number rather than the Bond number. This gives [25]

$$S_{nw}^{F} \propto F^{\frac{-\nu(D-d)}{1+\nu}}. \quad (21)$$

According to the theoretical prediction, Eq. (21) becomes

$$S_{nw}^{F} \propto F^{0.10}, \quad (22)$$

and using Eqs. (19) and (22), we obtain

$$S_{nw}^{F} \propto \epsilon^{-0.10}, \quad (23)$$

To verify these theoretical predictions, Figs. 9(a) and 9(b) illustrate the nonwetting phase saturation obtained experimentally and numerically as a function of the level of disorder $\epsilon$ and the $F$ number. In these plots, guide-to-the-eye lines corresponding to the theoretical exponents are also shown.

In both figures, we can observe that our experimental and simulated data satisfactorily align with the theoretically predicted scaling, thereby reinforcing the robustness of our findings and substantiating the validity of the theoretical framework. Nonetheless, it is relevant to mention that the scaling presented in Eq. (21) is extremely sensitive to the observed fractal dimension of the nonwetting percolating cluster. Based on slow drainage experiments using quasi-2D etched micromodels, Lenormand and Zarcone [26] reported $D$ ranging from 1.80 to 1.84, which would lead to $\varphi$ values between 0.091 and 0.114. These variations could be linked to the small deviations between theoretical and measured values presented in Fig. 9.

### D. Estimation of fluctuation number from more general pore structures

Even though the experiments and simulations in this work used a (disordered) honeycomb pore structure with constant pore connectivity, the theory presented is in principle valid to more general pore structures and the specific details of the pore geometry used here do not enter the derivation of the scaling relationship in Eq. (11). The fluctuation number calculated from Eq. (16) is valid in the case of any slow drainage process and can also be generalized for the case of faster flows where viscous effects are also present, as seen in Eq. (9) [25,60]. Although the calculation of $N(p_c)$ in principle depends on knowledge of the pore connectivity
<figure><img src="image_12.png" /><figcaption>FIG. 9. The final saturation ($S_{\text{nw}}^F$) of nonwetting fluid is plotted as a function of $\epsilon$ in (a) and F number in (b). Similarly to Fig. 7, symbols denote individual experiments and simulations. The presence of a black solid line signifies the theoretically predicted slope.</figcaption></figure>

via the bond percolation threshold, the approximation to this quantity given in Eq. (14) is expected to produce reasonable values as long as the distribution of pore-throat sizes is not too far from a uniform distribution. This result is potentially very useful in scenarios where the specific connectivity of a sample is not known but the pore-throat size distribution is obtained via other indirect methods such as mercury intrusion porosimetry [76,77] or gas adsorption techniques [78,79]. If the methodology used gives access to the width $W_p$ of the capillary pressure thresholds distribution, then Eq. (14) can be applied directly to estimate $N(p_c)$. If only the with $W_r$ of the pore-throat size distribution is known, then the width of the pressure threshold distribution can be determined from the Young-Laplace equation as $W_p = W_r 4\gamma \cos\theta / a^2$ prior to using Eq. (15). After computing $N(p_c)$ from Eq. (14), Eqs. (9) and (11) can then be applied to compute the fluctuation number $F$ and the scaling behavior of the invasion front width. This finding enables a straightforward estimation of how pore disorder influences the invasion front width in systems where only the average pore-throat size $a$ and the width of the pore-throat size distribution $W_r$ are known.

V. CONCLUSIONS

In this paper, we present a theoretical, experimental, and numerical analysis of the influence of pore-network disorder on invasion patterns during slow drainage in a porous medium. By controlling the disorder via a single parameter,

we explored how the invasion front width, maximum trapped cluster length, and nonwetting phase saturation are influenced by the underlying disorder. Our experimental and simulation results show that the invasion structure, including the residual fluid left behind the fluid front, is primarily controlled by the dimensionless fluctuation number $F$. This quantity has been explored in the past [25,30,35,60] in connection with variations in viscous and gravitational forces and systematic gradients in pore geometry. Our findings indicate that the invasion front width is largely determined by the size of the trapped clusters. In the 2D case, the experimental and simulation results are consistent with the theoretically predicted scaling behavior.

This analysis highlights the significant impact of pore-scale disorder on slow-draining processes in porous media. By establishing scaling laws and analyzing pore-scale disorder, we provide a framework for predicting and controlling multiphase flow patterns in a wide range of applications. Future research can extend these findings by investigating more complex geometries, such as 3D systems, mixed wetting conditions, and heterogeneous porous media.

ACKNOWLEDGMENTS

We acknowledge Tom Vincent-Dospital, Kevin Pierce, Per Arne Rikvold, Federico Lanza, and Fredrik Kvalheim Eriksen for useful discussions. We also thank the support from the University of Oslo and from the Research Council of Norway, through Projects No. 262644 and No. 324555.

[1] R. Anderson, L. Zhang, Y. Ding, M. Blanco, X. Bi, and D. P. Wilkinson, A critical review of two-phase flow in gas flow channels of proton exchange membrane fuel cells, *J. Power Sources* **195**, 4531 (2010).

[2] S. H. Hwang, C. Kim, H. Song, S. Son, and J. Jang, Designed architecture of multiscale porous TiO$_2$ nanofibers for dye-sensitized solar cells photoanode, *ACS Appl. Mater. Interfaces* **4**, 5287 (2012).

[3] L. Lake, *Enhanced Oil Recovery* (Prentice Hall, Englewood Cliffs, NJ, 1989).

[4] F. M. Orr and J. J. Taber, Use of carbon dioxide in enhanced oil recovery, *Science* **224**, 563 (1984).

[5] K. V. Veld and O. R. Phillips, The economics of enhanced oil recovery: Estimating incremental oil supply and Co$_2$ demand in the powder river basin, *Energy J. 31*, 31 (2010).
[6] M. Shabani Afrapoli, S. Alipour, and O. Torsæter, Fundamental study of pore scale mechanisms in microbial improved oil recovery processes, *Transp. Porous Media* **90**, 949 (2011).

[7] J. Yan, X. Luo, W. Wang, R. Toussaint, J. Schmittbuhl, G. Vasseur, F. Chen, A. Yu, and L. Zhang, An experimental study of secondary oil migration in a three-dimensional tilted porous medium, *AAPG Bull.* **96**, 773 (2012).

[8] G. Vasseur, X. Luo, J. Yan, D. Loggia, R. Toussaint, and J. Schmittbuhl, Flow regime associated with vertical secondary migration, *Mar. Pet. Geol.* **45**, 150 (2013).

[9] R. L. Payton, M. Fellgett, B. L. Clark, D. Chiarella, A. Kingdon, and S. Hier-Majumder, Pore-scale assessment of subsurface carbon storage potential: implications for the UK Geoenergy Observatories project, *Pet. Geosci.* **27**, petgeo2020-092 (2021).

[10] J. Lipiec, J. Kuś, A. Słowińska-Jurkiewicz, and A. Nosalewicz, Soil porosity and water infiltration as influenced by tillage methods, *Soil Tillage Res.* **89**, 210 (2006).

[11] N. Tufenkji and M. Elimelech, Correlation equation for predicting single-collector efficiency in physicochemical filtration in saturated porous media, *Environ. Sci. Technol.* **38**, 529 (2004).

[12] J.-D. Chen and D. Wilkinson, Pore-scale viscous fingering in porous media, *Phys. Rev. Lett.* **55**, 1892 (1985).

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

[14] R. Lenormand, E. Touboul, and C. Zaarcone, Numerical models and experiments on immiscible displacement in porous media, *J. Fluid Mech.* **189**, 165 (1988).

[15] B. Zhao, C. W. Mac Minn, and R. Juanes, Wettability control on multiphase flow in patterned microfluidics, *Proc. Natl. Acad. Sci. USA* **113**, 10251 (2016).

[16] R. Lenormand, Flow through porous media: Limits of fractal pattern, *Proc. R. Soc. Lond. A* **423**, 159 (1989).

[17] B. Sandnes, H. A. Knudsen, K. J. Måløy, and E. G. Flekkøy, Labyrinth patterns in confined granular-fluid systems, *Phys. Rev. Lett.* **99**, 038001 (2007).

[18] B. Sandnes, E. G. Flekkøy, K. J. Måløy, and H. See, Patterns and flow in frictional fluid dynamics, *Nat. Commun.* **2**, 288 (2011).

[19] B. B. Mandelbrot, *The Fractal Geometry of Nature* (W. H. Freeman, San Francisco, CA, 1982).

[20] J. Feder, *Fractals* (Plenum, New York, 1988).

[21] D. A. Weitz, J. P. Stokes, R. C. Ball, and A. P. Kushnick, Dynamic capillary pressure in porous media: Origin of the viscous-fingering length scale, *Phys. Rev. Lett.* **59**, 2967 (1987).

[22] G. Løvoll, Y. Mèheust, R. Toussaint, J. Schmittbuhl, and K. J. Måløy, Growth activity during fingering in a porous Hele Shaw cell, *Phys. Rev. E* **70**, 026301 (2004).

[23] K. T. Tallakstad, H. A. Knudsen, T. Ramstad, G. Løvoll, K. J. Måløy, R. Toussaint, and E. G. Flekkøy, Steady-state two-phase flow in porous media: Statistics and transport properties, *Phys. Rev. Lett.* **102**, 074502 (2009).

[24] K. T. Tallakstad, G. Løvoll, H. A. Knudsen, T. Ramstad, E. G. Flekkøy, and K. J. Måløy, Steady-state, simultaneous two-phase flow in porous media: An experimental study, *Phys. Rev. E* **80**, 036308 (2009).

[25] K. J. Måløy, M. Moura, A. Hansen, E. G. Flekkøy, and R. Toussaint, Burst dynamics, upscaling and dissipation of slow drainage in porous media, *Front. Phys.* **9**, 796019 (2021).

[26] R. Lenormand and C. Zarcone, Invasion percolation in an etched network: Measurement of a fractal dimension, *Phys. Rev. Lett.* **54**, 2226 (1985).

[27] L. Furuberg, J. Feder, A. Aharony, and T. Jøssang, Dynamics of invasion percolation, *Phys. Rev. Lett.* **61**, 2117 (1988).

[28] M. Moura, K. J. Måløy, and R. Toussaint, Critical behaviour in porous media flow, *Europhys. Lett.* **118**, 14004 (2017).

[29] B. Xiao, Q. Huang, H. Chen, X. Chen, and G. Long, A fractal model for capillary flow through a single tortuous capillary with roughened surfaces in fibrous porous media, *Fractals* **29**, 2150017 (2021).

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

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

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

[33] V. Frette, J. Feder, T. Jøssang, and P. Meakin, Buoyancy fluid migration in porous media, *Phys. Rev. Lett.* **68**, 3164 (1992).

[34] G. Wagner, A. Birovljev, P. Meakin, J. Feder, and T. Jøssang, Fragmentation and migration of invasion percolation clusters: Experiments and simulations, *Phys. Rev. E* **55**, 7015 (1997).

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

[36] M. Muharrik, T. Suekane, and A. Patmonoaji, Effect of buoyancy on fingering growth activity in immiscible two-phase flow displacement, *J. Fluid Sci. Technol.* **13**, JFST0006 (2018).

[37] C. Cottin, H. Bodiguel, and A. Colin, Influence of wetting conditions on drainage in porous media: A microfluidic study, *Phys. Rev. E* **84**, 026311 (2011).

[38] R. Holtzman and E. Segre, Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore filling, *Phys. Rev. Lett.* **115**, 164501 (2015).

[39] R. Holtzman, M. Dentz, R. Planet, and J. Ortin, The origin of hysteresis and memory of two-phase flow in disordered media, *Commun. Phys.* **3**, 222 (2020).

[40] P. de Gennes, F. Brochard-Wyart, and D. Quéré, *Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves* (Springer, New York, 2004).

[41] B. K. Primkulov, S. Talman, K. Khaleghi, A. Rangriz Shokri, R. Chalaturnyk, B. Zhao, C. W. MacMinn, and R. Juanes, Quasistatic fluid-fluid displacement in porous media: Invasion-percolation through a wetting transition, *Phys. Rev. Fluids* **3**, 104001 (2018).

[42] B. K. Primkulov, A. A. Pahlavan, X. Fu, B. Zhao, C. W. MacMinn, and R. Juanes, Signatures of fluid-fluid displacement in porous media: Wettability, patterns and pressures, *J. Fluid Mech.* **875**, R4 (2019).

[43] H. S. Rabbani, D. Or, Y. Liu, C.-Y. Lai, N. B. Lu, S. S. Datta, and H. A. Stone, Suppressing viscous fingering in structures porous media, *Proc. Natl. Acad. Sci. USA* **115**, 4833 (2018).

[44] N. B. Lu, C. A. Browne, D. B. Amchin, J. K. Nunes, and S. S. Datta, Controlling capillary fingering using pore size gradients in disordered media, *Phys. Rev. Fluids* **4**, 084303 (2019).
[45] K. Singh, H. Scholl, M. Brinkmann, M. D. Michiel, M. Scheel, S. Herminghaus, and R. Seemann, The role of local instabilities in fluid invasion into permeable media, *Sci. Rep.* **7**, 444 (2017).

[46] J. Bear, *Dynamics of Fluids in Porous Media* (Elsevier, New York, 1972).

[47] R. Toussaint, G. Løvoll, Y. Méheust, K. J. Måløy, and J. Schmittbuhl, Influence of pore-scale disorder on viscous fingering during drainage, *Europhys. Lett.* **71**, 583 (2005).

[48] G. Løvoll, M. Jankov, K. Måløy, R. Toussaint, J. Schmittbuhl, G. Schäfer, and Y. Méheust, Influence of viscous fingering on dynamic saturation-pressure curves in porous media, *Transp. Porous Media* **86**, 305 (2011).

[49] M. Cieplak and M. O. Robbins, Dynamical transition in quasistatic fluid invasion in porous media, *Phys. Rev. Lett.* **60**, 2042 (1988).

[50] M. Cieplak and M. O. Robbins, Influence of contact angle on quasistatic fluid invasion of porous media, *Phys. Rev. B* **41**, 11508 (1990).

[51] R. Hu, T. Lan, G.-J. Wei, and Y.-F. Chen, Phase diagram of quasi-static immiscible displacement in disordered porous media, *J. Fluid Mech.* **875**, 448 (2019).

[52] R. Holtzman, Effects of pore-scale disorder on fluid displacement in partially-wettable porous media, *Sci. Rep.* **6**, 36221 (2016).

[53] D.-S. Wu, R. Hu, T. Lan, and Y.-F. Chen, Role of pore-scale disorder in fluid displacement: Experiments and theoretical model, *Water Resour. Res.* **57**, 16 (2021).

[54] R. G. Larson, L. E. Scriven, and H. T. Davis, Percolation theory of two phase flow in porous media, *Chem. Eng. Sci.* **36**, 57 (1981).

[55] R. Chandler, J. Koplik, K. Lerman, and D. Willemsen, Capillary displacement and percolation in porous media, *J. Fluid Mech.* **119**, 249 (1982).

[56] D. Wilkinson and J. F. Willemsen, Invasion percolation: A new form of percolation theory, *J. Phys. A: Math. Gen.* **16**, 3365 (1983).

[57] O. I. Frette, K. J. Måløy, and J. Schmittbuhl, Immiscible displacement of viscosity-matched fluids in two-dimensional porous media, *Phys. Rev. E* **55**, 2969 (1997).

[58] E. Aker, K. J. Måløy, and A. Hansen, Viscous stabilization of 2D drainage displacements with trapping, *Phys. Rev. Lett.* **84**, 4589 (2000).

[59] E. Aker, K. J. Måløy, and A. Hansen, Dynamics of stable viscous displacement in porous media, *Phys. Rev. E* **61**, 2936 (2000).

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

[61] H. S. Hele-Shaw, The flow of water, *Nature* (London) **58**, 34 (1898).

[62] D. T. Lee and B. J. Schachter, Two algorithms for constructing a Delaunay triangulation, *Int. J. Comput. Inform. Sci.* **9**, 219 (1980).

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

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

[65] P. Concus and R. Finn, On the behavior of a capillary surface in a wedge, *Proc. Natl. Acad. Sci. USA* **63**, 292 (1969).

[66] P. Concus and R. Finn, On capillary free surfaces in the absence of gravity, *Acta Math.* **132**, 177 (1974).

[67] Q.-Z. Huang, S.-Y. Hsu, H. Chen, C.-W. Tsao, K. Lamorski, and C. Sławinúski, Old and new wetting liquids separation in grain-based pattern micromodel during wetting cycles, *Int. Agrophys.* **38**, 447 (2024).

[68] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.7.023040 for Experimental Movie (Movie S1) and Simulation Movie (Movie S2) for disorder levels 0.1, 0.5, and 0.9.

[69] D. Stauffer, *Introduction to Percolation Theory* (Taylor & Francis, London, 1994).

[70] S. Roux and E. Guyon, Temporal development of invasion percolation, *J. Phys. A* **22**, 3693 (1989).

[71] M. Moura, E. Florentino, K. Måløy, G. Schäfer, and R. Toussaint, Impact of sample geometry on the measurement of pressure-saturation curves: Experiments and simulations, *Water Resour. Res.* **51**, 8900 (2015).

[72] J.-F. Gouyet, M. Rosso, and B. Sapoval, Fractal structure of diffusion and invasion fronts in three-dimensional lattices through the gradient percolation approach, *Phys. Rev. B* **37**, 1832 (1988).

[73] E. Sapoval, M. Rosso, and J. Gouyet, The fractal nature of a diffusion front and the relation to percolation, *J. Phys Lett.* **46**, 149 (1985).

[74] Z. Wang, K. Chauhan, J.-M. Pereira, and Y. Gan, Disorder characterization of porous media and its effect on fluid displacement, *Phys. Rev. Fluids* **4**, 034305 (2019).

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

[76] E. W. Washburn, Note on a method of determining the distribution of pore sizes in a porous material, *Proc. Natl. Acad. Sci. USA* **7**, 115 (1921).

[77] R. Lal and M. K. Shukla, *Principles of Soil Physics* (CRC Press, London, 2019).

[78] E. P. Barrett, L. G. Joyner, and P. P. Halenda, The determination of pore volume and area distributions in porous substances. I. computations from nitrogen isotherms, *J. Am. Chem. Soc.* **73**, 373 (1951).

[79] K. S. W. Sing, Reporting physisorption data for gas/solid systems with special reference to the determination of surface area and porosity (provisional), *Pure Appl. Chem.* **54**, 2201 (1982).