# A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media

Zhongzheng Wang, Jean-Michel Pereira, Emilie Sauret, Saman Aryana, Zhang Shi, Yixiang Gan

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

## To cite this version:

Zhongzheng Wang, Jean-Michel Pereira, Emilie Sauret, Saman Aryana, Zhang Shi, et al.. A pore-resolved interface tracking algorithm for simulating multiphase flow in arbitrarily structured porous media. Advances in Water Resources, 2022, 162, pp.104152. 10.1016/j.advwatres.2022.104152. hal-03905451

HAL Id: hal-03905451

https://enpc.hal.science/hal-03905451

Submitted on 18 Dec 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
A pore-resolved interface tracking algorithm for simulating
multiphase flow in arbitrarily structured porous media

Zhongzheng Wang,<sup>1,2,3</sup> Jean-Michel Pereira,<sup>3</sup> Emilie Sauret,<sup>1</sup>

Saman A. Aryana,<sup>4</sup> Zhang Shi,<sup>2,5</sup> and Yixiang Gan<sup>2,6,*</sup>

<sup>1</sup>*School of Mechanical, Medical and Process Engineering, Faculty of Engineering, Queensland University of Technology, QLD 4001, Australia*

<sup>2</sup>*School of Civil Engineering, The University of Sydney, NSW 2006, Australia.*

<sup>3</sup>*Navier, Ecole des Ponts, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France.*

<sup>4</sup>*Department of Chemical Engineering, University of Wyoming, Laramie, WY 82071, USA*

<sup>5</sup>*School of Mechanical and Mining Engineering, The University of Queensland, QLD 4072, Australia*

<sup>6</sup>*The University of Sydney Nano Institute (Sydney Nano), The University of Sydney, NSW 2006, Australia*
Abstract

We present an interface tracking algorithm that incorporates essential pore-scale invasion mechanisms for simulating multiphase flow in porous media in capillary dominated regime. The algorithm allows the capture of the actual shape of porous medium with complex geometrical features, including both convex and concave particles or irregular channels. It also incorporates the sharp edge pinning effect such that pinning of menisci at corners can be naturally taken into account. We firstly conduct simulations under imbibition and drainage conditions in a single junction micro-model and in disordered porous media for validation. It is shown that the invasion mechanisms at pore scale are accurately captured, and the macroscopic invasion morphology in porous media are reproduced. A series of simulations for fluid displacement processes are carried out across a wide range of intrinsic contact angles in a two-dimensional representation of a Berea sandstone. The transition between stable displacement and capillary fingering is characterized using displacement efficiency, fractal dimension, and fluid-fluid interfacial length. This work deepens the understanding of fluid-fluid displacement processes in porous media. Furthermore, the presented interface tracking algorithm offers a highly efficient tool for further exploration of effects of wettability, geometry, and topology on multiphase flow in realistic porous materials.

Keywords: Multiphase flow, interface tracking algorithm, pore-network model, wettability, porous media

* Corresponding author. Email: yixiang.gan@sydney.edu.au
# I. INTRODUCTION

The phenomenon where one fluid displaces another in porous materials is ubiquitous in many natural and industrial processes, such as water infiltration into soils [1], hydrogen fuel cells [2], enhanced oil recovery [3, 4], and carbon geological sequestration [5, 6]. Understanding the mechanisms of multiphase flow in porous media is of vital importance for prediction and control of fluid transport in various engineering problems. It has been established that the pattern formation of fluid-fluid displacement depends on both fluid properties and flow conditions [7–12], which can be described by the Capillary number, reflecting the relative importance of viscous forces to capillary forces, and viscosity ratio. Depending on the viscosity of fluids, interfacial tension, and injection velocity, different invasion morphologies including capillary fingering, viscous fingering, and stable displacement have been observed numerically and experimentally, both in 2D microfluidics (e.g., Hele-Shaw cells) and 3D granular materials [8, 13–18]. At the same time, the wetting condition of the solid skeleton, characterized by the contact angle formed at the fluid-fluid-solid triple line, can also greatly impact the pattern formation during the fluid-fluid displacement processes [13, 14, 19–30]. Specifically, when the solid shifts from non-wetting (drainage, contact angle $θ > 90°$) to wetting (imbibition, $θ < 90°$) to the invading fluid, the invasion front is smoothed, corresponding to the transition from viscous/capillary fingering towards compact displacement.

Pore-network models have been extensively applied to study multiphase flow in porous media due to their favorable computational cost compared with conventional computational methods, i.e., Navier Stokes solvers [11, 22, 27, 29–34]. Pore-network model are rigorous tools for reproducing macroscopic invasion morphologies, and consequently revealing potential connections between pore-scale mechanisms and macroscopic observations. In the seminal work by Cieplak and Robbins [22, 33], an interface-tracking algorithm considering different pore-filling events was shown to capture the transition from capillary fingering to compact displacement. The *overlap* event, during which two menisci merge and form a new meniscus, is found to be the key mechanism that smooths the invasion front and suppresses trapping of defending phase, leading to an overall more stable displacement. The interface-tracking algorithm originated from the work by Cieplak and Robbins [22, 33] has been adopted, or extended, to study the impact of wettability, viscosity, disorder, and corner flows during fluid-fluid displacement processes [27, 29, 34, 35]. For example, by adding the dynamic
effect, from the pore-network simulation results across a wide range of capillary numbers, the likelihood of the overlap event is found to strongly correlate with total fluid-fluid interfacial area and finger width [27].

However, one limitation of the aforementioned pore-network models is that they are only applicable to porous media composed of perfectly spherical grains, whereas in most natural systems such as sand packs and surrogate models such as microfluidics, solid surface with irregular shape or sharp edges are often encountered. Due to these non-smooth surfaces, pinning of menisci during displacement process may be observed [36, 37], which leads to effective contact angles deviating from the intrinsic one, altering the capillary pressure resistance at local pore/throat. The sharp edge pinning phenomenon has also been studied by Wu et al [38] (termed “capillary valve effect”). It was demonstrated that the pinning of meniscus needs to be considered in order to obtain accurate details of two-phase flows, especially during imbibition processes. Recently, an extended pore-network model (EPONM) with consideration of particle shape was proposed to investigate immiscible two-phase flows in the capillary dominated regime [30]. It is shown that the presence of sharp corners, and the subsequent meniscus pinning, leads to an effective contact angle that deviates from the intrinsic one. Nevertheless, only particles represented by simple polygons arranged on a regular triangular lattice are considered [30].

In this work, we further extend the interface tracking algorithm for multiphase flow, making it applicable to porous media of arbitrary pore structures and grain shapes such as reconstructed pore geometries from 2D scans of rock samples. We first present the algorithm along with an explanation of the physical mechanisms that are captured by the method. This is followed by simulations of two-phase flow in single junction micro-models and disordered porous media for validation purpose. Then, to demonstrate the capability of the algorithm, we carry out systematic simulations in representative porous media that are constructed based on a Berea sandstone across a wide range of contact angles. We analyze the effects of wettability on multiphase flow in this realistic porous medium by focusing on a few important metrics that are often used to characterize the displacement patterns, including sweep efficiency, fractal dimension, and fluid-fluid interfacial area.
<figure><img src="image_3.png" /><figcaption>FIG. 1. **Mechanisms of menisci advancement.** (A) Meniscus with intrinsic contact angle $\theta_0$ gets pinned at the sharp corner with the corner angle $\alpha$. The meniscus adjusts shape without contact line movement until reaching an upper limit at $\theta_{\text{unpin}}$, then it continues to advance. (B) A *burst* event (left) and *touch* event (right). The blue-dashed line marks the meniscus shape associated with the maximum capillary pressure resistance. Menisci number increases by one after *burst* or *touch*. (C) An *overlap* event, during which two menisci (light blue) merge into one (dark blue). (D) A local snapshot of (C). The bottom meniscus reaches its critical unpin contact angle marked by red-dashed lines, moving upward (from light blue to dark blue) to the next computational node (black dots) after which an *overlap* event is detected and executed.</figcaption></figure>

## II. METHODS

For a given two-dimensional porous medium represented in a grayscale image, the grain boundary surfaces are detected and discretized into computational nodes (i.e., black dots in Fig. 1(D)) according to a prescribed mesh size. This is followed by the initialization of a set of arcs with radius $r$ representing the menisci (fluid-fluid interfaces). The initial value of $r$ is chosen according to the local pore geometry and the effective contact angle $\theta$, where the effective contact angle is the angle measured within the invading phase at the triple contact point, which is equal to the intrinsic contact angle $\theta_0$ at perfectly smooth surface, and could be larger than $\theta_0$ at sharp corners [Fig. 1(A)] due to the sharp edge pinning effect as we will discuss later. For every meniscus, different potential advancement events are considered. Here, we consider (i) a *burst* event, corresponding to the capillary pressure associated with meniscus shape indicated by the blue-dashed lines in Fig. 1(B); (ii) a *touch* event, corresponding to the capillary pressure associated with meniscus shape indicated by the blue-solid lines in Fig. 1(B); and (iii) an *unpin* event (for both left and right triple points), corresponding to the capillary pressure when the effective contact angle reaches the
<figure><img src="image_4.png" /><figcaption>FIG. 2. **Schematic showing the determination of advancement events.** (A) The orange and purple circles are associated with the menisci shapes for *unpin* events at the right (with the invading direction as the frame of reference, i.e., black arrow) and left triple-contact points, respectively. (B) The blue-solid and green-solid curves are associated with the menisci shape for a *touch* and a *burst* event, respectively. The red rectangle is the search box, with all computation nodes in red color being considered as potential targets for the *touch* event, and the point marked by purple square is determined as the target *touch* point, which makes the largest angle with two triple points. In the case where the middle grain is absent, the meniscus shape for a *touch* event is marked by the blue-dashed line. (C) To determine whether a region is trapped after the advancement of a meniscus, a searching path following the clockwise direction is initialized from the left triple point. The search path follows dry nodes, jumps to the left triple point when reaching a right triple point, and stops when reaching the outlet. A region is labelled trapped if the path returns to the starting point. Menisci 1 and 2 show a trapped and an active meniscus, respectively.</figcaption></figure>

unpin angle $\theta_{\text{unpin}} = \theta_0 + (180^\circ - \alpha)$ [Fig. 1(A)]. This is to capture the sharp edge pinning phenomenon where the equilibrium state of effective contact angle $\theta$ cannot exceed $\theta_{\text{unpin}}$ according to a purely geometrical extension of Young-Dupre equation [36, 37].

Fig. 2(A-B) shows in more detail the determination of meniscus advancement events. In Fig. 2(A), the circle associated with the *unpin* event at the right triple point (with the invading direction as the frame of reference, i.e., black arrow) is shown in orange, and the circle associated with the left triple-contact point are shown in purple. In this scenario, as the pressure within the invading fluid builds up and the meniscus bends forward, the *unpin* event at the right triple point (in orange) occurs first. In Fig. 2(B), the blue-solid and green-
solid curves represent the menisci shapes for a *touch* and a *burst* event, respectively. The red box denotes the search box with red dots inside be considered as potential targets for the *touch* event. The purple square is determined to be the target *touch* point that makes the largest angle with two triple points. Note that the size of the search box should be greater than the largest pore size of the porous medium. For the *burst* event, the meniscus (green-solid curve) is a semicircle that corresponds to the smallest circle (largest pressure) for the given triple points. In this case, the *touch* event occurs before the *burst* event. The critical advancement pressure $P_c$ for different events is calculated according to $P_c = \gamma/r$, with $\gamma$ the interfacial tension and $r$ the meniscus radius. The final $P_c$ is chosen according to the event that occurs first (not necessarily the one with greatest/smallest $P_c$ values). For the meniscus in Fig. 2(A-B), the final $P_c$ is the one associated with the *unpin* event at the right triple point. Note that, if the middle grain in Fig. 2(B) was absent, the corresponding shape for a *touch* event would be represented by the blue-dashed curve. This *touch* event clearly takes place after the *burst* event, i.e., a post-burst event, and the $P_c$ associated with this event is thus overwritten by the *burst* pressure (the meniscus has to reach burst state before finally touching the grain).

After updating the critical capillary pressure at each computational step, the meniscus with minimum $P_c$ advances. The determination of the critical capillary pressure $P_c$ for a given meniscus considering different types of instability events is based on the instability event that occurs first based on the local pore geometry and contact angle, regardless of the value of $P_c$. Globally, however, the meniscus with the smallest resistance (based on $P_c$ value) advances. After each advancement, (iv) trapping of the defending phase is examined [Fig. 2(C)] and (v) an *overlap* event is checked and executed [Fig. 1(C)]. For a meniscus after an advancement, the trapping of the defending phase is checked by initializing a search path from the left triple point. The path follows the dry nodes [gray dots in Fig. 1(C)] in the clockwise direction, jumps to the left triple point once it reaches a right triple point, and stops when it reaches the outlet. A region is labelled trapped if the path returns to the starting point. Menisci 1 and 2 in Fig. 2(C) show a trapped and an active meniscus, respectively. Once a trapping event occurs, all menisci associated with the trapped region are deactivated to prevent further movement.

Fig. 3(A) shows the flow chart of the algorithm. A sub-region of the grayscale image of a representation of a Berea sandstone [15], and the corresponding snapshot of local menisci
configuration during fluid displacement simulation, are shown in Fig. 3(B) and Fig. 3(C), respectively. Specifically, the gray, blue, and black dots represent dry, wetted, and trapped surfaces, respectively. Since the algorithm is based on discretization of pore space (more specifically, discrete grain boundaries) into computational nodes with a certain mesh size, it is necessary to conduct a mesh refinement study. As the menisci advancement is fundamentally governed by the determined capillary pressures associated with different invasion modes, as a rule of thumb, the mesh size, i.e., the distance between adjacent nodes, should be much smaller than the grain size and local throat size. For simulations in this work, the mesh refinement test is conducted, and it is found that having, on average, more than 25 nodes per grain for the porous media in Fig. 4(B) and more than 50 nodes per grain for porous medium in Fig. 5(A) are sufficient, and a further increase in the number of nodes per grain does not lead to a significant change in the results.

There are two major differences between the current algorithm and the interface tracking algorithm (and its derivatives) originally proposed by Cieplak and Robbins [22, 33]. Firstly, the advancement of menisci in the current algorithm consists of a node-to-node movement, whereas the original one involves a pore-by-pore movement. The advantage of a node-to-node movement is its ability to capture the actual shape of a porous medium with complex geometrical features, including both convex and concave particles and irregular channels, whereas the original method is only applicable to porous media composed of perfect circular grains. Secondly, the sharp edge pinning effect is implemented such that pinning of the meniscus at corners can be naturally captured. This method provides a rigorous tool for studying two-phase flow in porous media, especially for investigating the combined impacts of topology, geometry, and wettability.

To validate the algorithm, we first conduct fluid displacement simulation in a single square junction under both drainage ($\theta_0 = 150^\circ$) and imbibition ($\theta_0 = 60^\circ$) conditions. This junction has a square pore body with edge length of unity, which is connected to four throats with widths of 0.1, 0.2, 0.3, and 0.4 on the top, right, bottom and left sides, respectively. For imbibition, to consider different pore filling mechanisms as in [39, 40], three types of boundary conditions are considered (Fig. 4(A)) with the number of inlets ranging from one to three. We also conduct simulations of multiphase flow in disordered porous media and compare the results with those from the literature [34]. The domain is filled with circular particles of different radii placed on a triangular lattice, and the fluid is injected from the
<figure><img src="image_5.png" /><figcaption>FIG. 3. **Details on the interface tracking algorithm.** (A) Flow chart of the algorithm. (B) A sub-region of grayscale image of a representation of a Berea sandstone [15]. Grains and pore space are in gray and black colors, respectively. (C) A snapshot of local menisci configuration during simulation. Dots denote the computational nodes. Specifically, gray, blue, and black colors represent dry, wetted, and trapped surfaces, respectively. Blue curves show the location of menisci. Blue arrows mark the menisci advancement directions.</figcaption></figure>

left boundary with top and bottom boundaries treated as impermeable walls (Fig. 4(B)). The simulations stop when the invading phase reaches the outlet (the right boundary).

To demonstrate the capability of the algorithm, we use the realistic porous medium described in [15], which is shown in Fig. 5(A). The network represents the distribution of pores and throats observed in a thin section of Berea sandstone, and the detailed procedure
is described in [41]. Here, we consider a geometry relevant to oil recovery processes: a point
inlet at the bottom-right corner as the injection well and a point outlet at the top-left corner
as the withdrawal well (similar to [42]). Again, the simulations stop when the invading phase
reaches the outlet. For multiphase flow with extreme wetting condition where the intrinsic
contact angle is smaller than 45°, the corner flow phenomenon may start to appear [14, 29],
allowing invading fluid to propagate without filling the pore bodies. This mechanism is
currently not captured in the algorithm and we thus limit our attention to displacement
processes under normal imbibition to strong drainage, i.e, intrinsic contact angles $\theta_0 \in$
[45°, 180°] at 5° increments. It is worth noting that recent experimental observations in
microfluidic systems show that the critical intrinsic contact angle that marks the emergence
of corner flow can be different from 45° [43].

### III. RESULTS AND DISCUSSION

#### A. Fluid displacement in a single junction

Micro-models with simple geometries have been used to understand fluid displacement
mechanisms at the pore scale [38, 39, 44–46]. Here, we adopt a pore body with square shape,
which is connected to four throats of different widths as shown in Fig. 4(A). The outlet is
indicated by red lines. During the drainage process (intrinsic contact angle $\theta_0 = 150°$), the
invading fluid fills the pore first before advancing to the widest throat. Trapping of the
defending phase at corners is also observed within the pore-body (see Movie S1). These
observations are consistent with past experiments [39, 44]. For imbibition, we consider
three scenarios with one, two, or three inlets, which will be referred to as $I_1$, $I_2$, and $I_3$
imbibition, respectively. During $I_1$ imbibition (Movie S2), the invading phase partially fills
the pore body, implying that it does not have access to all throats that are connected to
the pore. As shown in Fig. 4(A), the invading phase eventually advances towards the left
throat as the narrowest throat (top throat) is inaccessible to the invading fluid throughout
the displacement process. This limited accessibility to throats during imbibition process was
also identified previously [47]. These observations demonstrate the need to take into account
the pore shape for correctly predicting invasion paths. Only considering sizes of pores and
throats, as often done in conventional pore filling rules, is insufficient for prediction of the
<figure><img src="image_6.png" /><figcaption>FIG. 4. **Validation of the interface tracking algorithm.** (A) Two-phase flow in a single square junction for drainage case (top row, $\theta_0 = 150^\circ$) and imbibition (bottom three rows, $\theta_0 = 60^\circ$). Red lines represent outlet. The simulation stops when the invading fluid reaches an outlet. (B) Two-phase flow in disordered porous media with circular grains. The invading phase is injected from the left boundary, with top and bottom boundaries treated as impermeable walls. The left column is the result from Hu et al [34], and the right column is from the present algorithm.</figcaption></figure>

correct invasion paths [45]. During $I_2$ imbibition (Movie S3), the fluid in smaller throat moves first due to a smaller capillary pressure until it gets pinned at the sharp corner. Then the invading fluid from the other throat advances. However, the unpin event takes place first at the widest throat. This is due to the concave-convex transition (effective contact angle $\theta < 90^\circ$ within the tube, but $\theta > 90^\circ$ at the sharp corner due to the pinning effect) of the meniscus shape, i.e., the larger throat offers a smaller capillary pressure resistance for the unpin event. Then, the overlap, or cooperative pore filling event [27, 33], takes place and the invading phase moves upwards once it touches the top throat. For $I_3$ imbibition (Movie S4), an asymmetrical sequence of throat-filling and unpin events is observed again, i.e., throat-filling in smaller throat occurs first while the unpin event happens in the larger throat first. Moreover, it is observed that a snap-off trapping event takes place before the invading phase advances into the top throat. Similar entrapment phenomenon has also been observed in a recent study during the imbibition process [46]. In summary, our results of drainage
and different types of imbibition processes in a single junction demonstrate the capability
of the current algorithm to fully resolve the combined impact of wettability and geometry
at the pore scale. This is attributed to the implementation of different types of invasion
mechanisms especially the sharp edge pinning phenomenon associated with the unpin event.
Some recent studies have highlighted the importance of considering this mechanism using
different names, such as “capillary valve effect” in [38], and “capillary barrier” in [45].

## B. Fluid displacement in circle packing

Fig. 4(B) shows the invasion morphology at percolation in disordered porous media with $\theta_0 = \{55^\circ, 75^\circ, 105^\circ\}$. The left columns in Fig. 4(B) are from Hu et al [34] using the conventional interface tracking algorithm [33] assuming perfect spheres, and the right columns are results from the current algorithm. Despite potential uncertainties during image processing and discretization, the transition of invasion morphology in compact, crossover, and fingering regimes is captured faithfully, and the shape of the overall invasion morphology in the original work is reproduced.

## C. Fluid displacement in realistic porous media

Fig. 5(A) shows the pore structure generated based on a thin section of Berea sandstone [15, 41]. Clearly, irregular grains with both convex and concave surfaces are present. Fig. 5(B) shows the color map of the local corner angles ($\alpha$ in Fig. 1(A)) of a sub-region enclosed by the red-box in Fig. 5(A). The corner angle encompasses the geometrical information of the local surface feature, which directly impacts the local pinning strength (upper bound of effective contact angle) calculated by $\theta_0 + (180^\circ - \alpha)$. Fig. 5(C) shows the distribution of the local corner angles for the entire simulation domain, indicating a wide range of $\alpha$ values with a slight shift below 180°. This indicates the overall concavity of grain shapes in the porous medium. Note that in the case of perfect spheres, one would expect a single value of 180° for the corner angle distribution.

The invasion morphologies of fluid displacement at percolation for different contact angles are shown in Fig. 6(A), which are qualitative demonstrations of the impact of wettability on multiphase flow in porous media (also SI Movie 5-7). As the grains or matrix become less
<figure><img src="image_7.png" /><figcaption>FIG. 5. Porous medium constructed by 2D imaging of Berea sandstone as per Guo and Aryana [15] (A) Pore geometry with gray and white colors representing grain and pore space, respectively. The point inlet and outlet are located at bottom-right and top-left corners, respectively. The boundaries are solid (impermeable) walls. (B) Corner angle map (α as in Fig. 1(A)) for the sub-region marked by the red box in (A). (C) Corner angle distribution of the whole porous medium. A slight shift below 180° indicates the overall concavity of grain shapes.</figcaption></figure>

wetting to the invading fluid, the displacement pattern shifts from a relatively compact region and a smooth interface, i.e., stable displacement, to a dendritic pattern with many more trapped regions, i.e., capillary fingering. The proportion of pore-scale instabilities, including *burst*, *touch*, and *overlap* events, during the entire displacement process are reported in Fig. 6(B), with insets showing the schematic of each event. At $\theta_0 \approx 80^\circ$, it is observed that the frequency of *overlap* events starts to decrease with the emergence of *burst* events. This is consistent with the established understanding that the *overlap* events tend to smooth the invasion front, whereas *burst* events are associated with patterns described by the theory of invasion percolation [48, 49].

To provide quantitative insights into the drastic transitions in displacement patterns
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. Fluid-fluid displacement in porous media (A) Final invasion morphologies under different intrinsic contact angles $\theta_0$. (B) Occurrence proportion of pore-scale events as a function of intrinsic contact angle $\theta_0$.</figcaption></figure>

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

<figure><img src="image_11.png" /><figcaption>FIG. 7. Characterization of pore filling by the invading fluid. (A) Displacement efficiency, defined as the ratio of displaced defending fluid to the pore space initially occupied by the defending fluid. The insets show two final displacement patterns at $\theta_0 = \{100^\circ, 160^\circ\}$. Only the invading phase is shown. (B) Fractal dimension values for different $\theta_0$ values are calculated using the box-counting method. Inset shows an example of box-number vs. box-size plot for the case $\theta_0 = 45^\circ$.</figcaption></figure>
observed in Fig. 6(A), we first consider the statistics of pore space filling status. This is characterised by two metrics: displacement efficiency and fractal dimension. The displacement efficiency, also called the sweep efficiency, measures the proportion of the defending fluid that is displaced out of the porous medium and is of great importance in applications such as oil recovery and carbon sequestration. Fig. 7(A) shows the displacement efficiency across a wide range of contact angles with insets showing the invasion phase in blue color at percolation for $\theta_0 = \{100^\circ, 160^\circ\}$. Despite the general decreasing tendency in displacement efficiency as contact angle $\theta_0$ increases, relatively large fluctuations are observed, implying the significant inherent uncertainties due to heterogeneities in realistic porous media. The huge fluctuations in the relationship between wettability and displacement efficiency are also observed in a recent study using the lattice Boltzmann method [50]. The fractal dimension $D_f$, on the other hand, from the theoretical perspective provides a multi-scale characterization of space filling. Fractal dimension values are calculated using the box-counting method [51] and the results are shown in Fig. 7(B) with the inset showing the box-number vs. box-size in log-log plot for $\theta_0 = 45^\circ$, where the absolute value of the slope is the fractal dimension. It can be seen that $D_f$ decreases from a value of 1.97 to 1.84, corresponding to the regimes of compact displacement and capillary fingering, respectively. These values are consistent with documented values of 1.96 and 1.83 for compact growth and invasion percolation, respectively [13, 14, 29, 47, 48, 52]. Unlike the sweep efficiency, the fluctuations in $D_f$ become smaller as the contact angle increases, indicating that $D_f$ is less sensitive to the heterogeneity of the porous medium and has stronger a correlation with the wetting condition.

Another key parameter for quantifying immiscible fluid displacement processes is the fluid-fluid interfacial area (length in 2D), which can be of interest for consideration of subsequent chemical reactions between different phases or analysis of interfacial potential energy. Here, we report the specific interfacial length calculated as the ratio of total fluid-fluid interfacial length $L$ divided by the area occupied by the invading phase, $A$, at percolation. As can be seen in Fig. 8(A), as $\theta_0$ increases, $L/A$ is relatively constant for small values of $\theta_0$ and starts to monotonically increase at $\theta_0 \approx 80^\circ$. The intrinsic contact angle at which $L/A$ starts to increase coincides with the emergence of the *burst* events as in Fig. 6(B). The inset reveals the average (with error bars representing the standard deviation) of meniscus size shown with blue circles and total number of menisci shown with red squares. It can be seen
<figure><img src="image_12.png" /><figcaption>FIG. 8. **Fluid-fluid interfacial length.** (A) The specific interfacial length of the invading fluid, calculated as the total fluid-fluid interfacial length *L* divided by the total invaded area *A*. Inset: left y-axis: average meniscus size (blue circles). Error bars represent standard deviation; right y-axis: number of menisci (red squares) for different θ₀. (B) Percentage of menisci with both of the fluid-fluid-solid triple points belonging to the same grain. Insets show two sub-regions where trapping occurs at θ₀ = {45°, 180°}. The single-grain-trapped menisci are marked by red lines. Note that straight lines are drawn here for simplicity.</figcaption></figure>

that the average meniscus size decreases, implying the menisci tend to be stuck in narrower throats as θ₀ increases, and the main contributing factor that leads to an overall increase in the specific interfacial length is the increase in the number of menisci. In Fig. 8(B), we show the percentage of menisci that are associated with a single-grain trapping (SGT) where the trapped region can be regarded as being situated within a dead end. This type of trapping is different from the trapped defending phase in the pore body among several grains since it can be much harder to be further mobilized by, for example, momentum transfer when the injection velocity is increased. The insets in Fig. 8(B) show two zoomed-in snapshots of trapped defending phase at imbibition (θ₀ = 45°) and drainage (θ₀ = 180°) conditions with menisci that belong to SGT marked by red lines (the menisci shape is drawn with straight lines for simplicity). As expected, SGT tends to take place at locations where the grain is convex, and the likelihood of the occurrence of SGT in drainage is much greater than in imbibition. Particularly, SGT in drainage can be also thought as the Cassie-Baxter wetting state, where the trapped defending phase in roughness grooves are observed [25, 53, 54].
Fig. 8(B) shows a non-monotonic relationship between SGT and $\theta_0$. We note that the occurrence of SGT depends on the geometrical feature of grain shape and the value of contact angle. An investigation of obstacle shape effect on SGT can be found in [55]. The importance of incorporating SGT phenomenon under drainage condition (by modifying the criterion for overlap event) is examined in [29].

## D. Discussion

We have presented a pore-resolved interface tracking algorithm for immiscible fluid-fluid displacement simulations. Through discretization of pore space and incorporation of essential pore-scale invasion mechanisms, the method is able to accurately track the interface motion during multiphase flow and reproduce the macroscopic invasion patterns in arbitrarily structured porous media, including reconstructed pore geometries from 2D scans of rock samples. Potential future extensions include adding viscous effects by, for example, using the moving capacitor model [35]. Other modifications for simulating multiphase flow in porous media characterized by heterogeneous wettability, i.e., with spatial variations in intrinsic contact angles, or consideration of gravitational effects by adding height and density dependent values to the critical capillary pressure, can be easily implemented. The present algorithm is highly efficient compared with other conventional computational fluid dynamics methods, e.g., Volume of Fluid (VOF) method or lattice Boltzmann method (LBM) [24, 56]. Most simulations reported in this work take less than one minute to finish on a laptop.

## IV. CONCLUSIONS

In this work, we present a rigorous interface tracking algorithm for simulating fluid displacement processes under quasi-static conditions. Essential pore-scale meniscus advance-ment mechanisms, including *burst*, *touch*, *overlap*, *unpin*, and *trapping* have been incorpo-rated. The algorithm is applicable to porous media with arbitrary pore structure for both imbibition and drainage conditions. By conducting multiphase flow simulations on a single junction and disordered porous media, we demonstrate the capability and validity of this algorithm for accurately capturing invasion mechanisms at pore scale and for reproducing the macroscopic invasion morphology.
We then conduct systematic simulations of fluid displacement processes across a wide range of intrinsic contact angles in a realistic porous medium that is a representation of a Berea sandstone. The transition from compact displacement to capillary fingering is demonstrated qualitatively and quantitatively. Important metrics, including displacement efficiency, fractal dimension, and interfacial length, are analyzed. Compared with displacement efficiency, smaller fluctuations are observed in fractal dimension and specific interfacial length, indicating that the latter two metrics are less sensitive to the heterogeneity of the porous medium and have a stronger correlation with the wetting condition of the matrix. This work provides insights into characterization of displacement patterns in porous media and the presented interface tracking algorithm offers an efficient tool for further exploration of effects of wettability, geometry, and topology on multiphase flow in realistic porous materials, which may facilitate the prediction and control of fluid transport phenomena in various engineering applications.

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

## ACKNOWLEDGEMENT

This work was financially supported by Australian Research Council (Projects DP170102886) and The University of Sydney SOAR Fellowship. YG acknowledges the financial support of Labex MMCD(ANR-11-LABX-022-01) for his stay at Laboratoire Navier at ENPC.

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

[2] Yanyao Bao and Yixiang Gan. Roughness effects of gas diffusion layers on droplet dynamics in pemfc flow channels. *International Journal of Hydrogen Energy*, 45(35):17869-17881, 2020.

[3] Larry W. Lake, Russell Johns, Bill Rossen, and Gary Pope. *Fundamentals of Enhanced Oil Recovery*. Society of Petroleum Engineers, 2014.
[4] Martin Blunt, F.John Fayers, and Franklin M. Orr. Carbon dioxide in enhanced oil recovery. *Energy Conversion and Management*, 34(9):1197-1204, 1993. Proceedings of the International Energy Agency Carbon Dioxide Disposal Symposium.

[5] Michael L. Szulczewski, Christopher W. MacMinn, Howard J. Herzog, and Ruben Juanes. Lifetime of carbon capture and storage as a climate-change mitigation technology. *Proceedings of the National Academy of Sciences*, 109(14):5185-5189, 2012.

[6] Juerg M. Matter, Martin Stute, Sandra Ó. Snæbjörnsdottir, Eric H. Oelkers, Sigurdur R. Gislason, Edda S. Aradottir, Bergur Sigfusson, Ingvi Gunnarsson, Holmfridur Sigurdardottir, Einar Gunnlaugsson, Gudni Axelsson, Helgi A. Alfredsson, Domenik Wolff-Boenisch, Kiflom Mesfin, Diana Fernandez de la Reguera Taya, Jennifer Hall, Knud Dideriksen, and Wallace S. Broecker. Rapid carbon mineralization for permanent disposal of anthropogenic carbon dioxide emissions. *Science*, 352(6291):1312-1314, 2016.

[7] Roland Lenormand, Eric Touboul, and Cesar Zarcone. Numerical models and experiments on immiscible displacements in porous media. *Journal of Fluid Mechanics*, 189(165-187), 1988.

[8] R Lenormand. Liquids in porous media. *Journal of Physics: Condensed Matter*, 2(S):SA79-SA88, December 1990.

[9] Y. C. Yortsos, B. Xu, and D. Salin. Phase diagram of fully developed drainage in porous media. *Phys. Rev. Lett.*, 79:4581-4584, December 1997.

[10] Ryan T. Armstrong, Apostolos Georgiadis, Holger Ott, Denis Klemin, and Steffen Berg. Critical capillary number: Desaturation studied with fast x-ray computed microtomography. *Geo-physical Research Letters*, 41(1):55-60, 2014.

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

[12] Yang Ju, Wenbo Gong, Wei Chang, and Min Sun. Effects of pore characteristics on water-oil two-phase displacement in non-homogeneous pore structures: A pore-scale lattice Boltzmann model considering various fluid density ratios. *International Journal of Engineering Science*, 154:103343, 2020.

[13] Mathias Trojer, Michael L. Szulczewski, and Ruben Juanes. Stabilizing fluid-fluid displacements in porous media through wettability alteration. *Phys. Rev. Applied*, 3:054008, May 2015.
[14] Benzhong Zhao, Christopher W. MacMinn, and Ruben Juanes. Wettability control on multiphase flow in patterned microfluidics. *Proceedings of the National Academy of Sciences*, 113(37):10251-10256, 2016.

[15] Feng Guo and Saman A. Aryana. An experimental investigation of flow regimes in imbibition and drainage using a microfluidic platform. *Energies*, 12(7), 2019.

[16] Kamaljit Singh, Hagen Scholl, Martin Brinkmann, Marco Di Michiel, Mario Scheel, Stephan Herminghaus, and Ralf Seemann. The role of local instabilities in fluid invasion into permeable media. *Scientific Reports*, 7(1):444, Mar 2017.

[17] Senyou An, Hamidreza Erfani, Omar E. Godinez-Brizuela, and Vahid Niasar. Transition from viscous fingering to capillary fingering: Application of gpu-based fully implicit dynamic pore network modeling. *Water Resources Research*, 56(12):e2020WR028149, 2020.

[18] Anindityo Patmonoaji, Mushlih Muharrik, Yingxue Hu, Chunwei Zhang, and Tetsuya Suekane. Three-dimensional fingering structures in immiscible flow at the crossover from viscous to capillary fingering. *International Journal of Multiphase Flow*, 122:103147, 2020.

[19] D. J. Crisp and W. H. Thorpe. The water-protecting properties of insect hairs. *Discuss. Faraday Soc.*, 3:210-220, 1948.

[20] W.R. Purcell. Interpretation of Capillary Pressure Data. *Journal of Petroleum Technology*, 2(08):11-12, 08 1950.

[21] Geoffrey Mason and Norman R. Morrow. Effect of contact angle on capillary displacement curvatures in pore throats formed by spheres. *Journal of Colloid and Interface Science*, 168(1):130-141, 1994.

[22] Marek Cieplak and Mark O. Robbins. Influence of contact angle on quasistatic fluid invasion of porous media. *Phys. Rev. B*, 41:11508-11521, June 1990.

[23] Michael Jung, Martin Brinkmann, Ralf Seemann, Thomas Hiller, Marta Sanchez de La Lama, and Stephan Herminghaus. Wettability controls slow immiscible displacement through local interfacial instabilities. *Phys. Rev. Fluids*, 1:074202, November 2016.

[24] Zhongzheng Wang, Kapil Chauhan, Jean-Michel Pereira, and Yixiang Gan. Disorder characterization of porous media and its effect on fluid displacement. *Phys. Rev. Fluids*, 4:034305, March 2019.

[25] Zhongzheng Wang, Jean-Michel Pereira, and Yixiang Gan. Effect of wetting transition during multiphase displacement in porous media. *Langmuir*, 36(9):2449-2458, 2020. PMID: 32070092.
[26] Zhang Shi, Zhongzheng Wang, and Yixiang Gan. Effects of topological disorder in unsaturated granular media via a pore-scale lattice Boltzmann investigation. *Advances in Water Resources*, 149:103855, 2021.

[27] Ran Holtzman and Enrico Segre. Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore filling. *Phys. Rev. Lett.*, 115:164501, October 2015.

[28] Hu Ran, Wan Jiamin, Yang Zhibing, Chen Yi-Feng, and Tokunaga Tetsu. Wettability and flow rate impacts on immiscible displacement: A theoretical model. *Geophysical Research Letters*, 45(7):3077-3086, 2018.

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

[30] Zhongzheng Wang, Jean-Michel Pereira, and Yixiang Gan. Effect of grain shape on quasi-static fluid-fluid displacement in porous media. *Water Resources Research*, 57(4):e2020WR029415, 2021.

[31] Martin J. Blunt. Physically-based network modeling of multiphase flow in intermediate-wet porous media. *Journal of Petroleum Science and Engineering*, 20(3):117-125, 1998.

[32] Martin J. Blunt. Flow in porous media — pore-network models and multiphase flow. *Current Opinion in Colloid & Interface Science*, 6(3):197-207, 2001.

[33] Marek Cieplak and Mark O. Robbins. Dynamical transition in quasistatic fluid invasion in porous media. *Phys. Rev. Lett.*, 60:2042-2045, May 1988.

[34] Ran Hu, Tian Lan, Guan-Ju Wei, and Yi-Feng Chen. Phase diagram of quasi-static immiscible displacement in disordered porous media. *Journal of Fluid Mechanics*, 875:448-475, 2019.

[35] Bauyrzhan K. Primkulov, Amir A. Pahlavan, Xiaojing Fu, Benzhong Zhao, Christopher W. MacMinn, and Ruben Juanes. Signatures of fluid-fluid displacement in porous media: wettability, patterns and pressures. *Journal of Fluid Mechanics*, 875:R4, 2019.

[36] J. Willard Gibbs. *The scientific papers*, volume 1. Dover Publications, New York, 1961.

[37] J.F Oliver, C Huh, and S.G Mason. Resistance to spreading of liquids by sharp edges. *Journal of Colloid and Interface Science*, 59(3):568-581, 1977.

[38] Rui Wu, Abdolreza Kharaghani, and Evangelos Tsotsas. Two-phase flow with capillary valve effect in porous media. *Chemical Engineering Science*, 139:241-248, 2016.
[39] R. Lenormand, C. Zarcone, and A. Sarr. Mechanisms of the displacement of one fluid by another in a network of capillary ducts. *Journal of Fluid Mechanics*, 135:337-353, 1983.

[40] *Role Of Roughness And Edges During Imbibition In Square Capillaries*, volume All Days of SPE Annual Technical Conference and Exhibition, 09 1984.

[41] Mohammed Alaskar. *In-Situ Multifunctional Nanosensors for Fractured Reservoir Characterization*. Theses, Stanford University, 2013.

[42] Lincoln Paterson. Diffusion-limited aggregation and two-fluid displacements in porous media. *Phys. Rev. Lett.*, 52:1621-1624, Apr 1984.

[43] Saeed Golmohammadi, Yi Ding, Matthias Küchler, Danny Reuter, Steffen Schlüter, Mohd Amro, and Helmut Geistlinger. Impact of wettability and gravity on fluid displacement and trapping in representative 2d micromodels of porous media (2d sand analogs). *Water Resources Research*, 57(10):e2021WR029908, 2021.

[44] Ioannis Zacharoudiou, Emily M. Chapman, Edo S. Boek, and John P. Crawshaw. Pore-filling events in single junction micro-models with corresponding lattice Boltzmann simulations. *Journal of Fluid Mechanics*, 824:550-573, 2017.

[45] Saideep Pavuluri, Julien Maes, Jianhui Yang, Mohamed Regaieg, Arthur Moncorgé, and Florian Doster. Towards pore network modelling of spontaneous imbibition: contact angle dependent invasion patterns and the occurrence of dynamic capillary barriers. *Computational Geosciences*, 24(2):951-969, Apr 2020.

[46] Chia-Wen Tsao, Qun-Zhan Huang, Chang-Ye You, Markus Hilpert, Shao-Yiu Hsu, Krzysztof Lamorski, Liang-Cheng Chang, and Cezary Sławiński. The effect of channel aspect ratio on air entrapment during imbibition in soil-on-a-chip micromodels with 2d and 2.5d pore structures. *Lab Chip*, 21:385-396, 2021.

[47] Martin J. Blunt. *Multiphase Flow in Permeable Media: A Pore-Scale Perspective*. Cambridge University Press, 2017.

[48] D Wilkinson and J F Willemsen. Invasion percolation: a new form of percolation theory. *Journal of Physics A: Mathematical and General*, 16(14):3365-3376, October 1983.

[49] Roland Lenormand and Cesar Zarcone. Invasion percolation in an etched network: Measurement of a fractal dimension. *Phys. Rev. Lett.*, 54:2226-2229, May 1985.

[50] Peter Mora, Gabriele Morra, Dave A. Yuen, and Ruben Juanes. Influence of wetting on viscous fingering via 2d lattice Boltzmann simulations. *Transport in Porous Media*, 138(3):511-538, 2019.
Jul 2021.

[51] B. B. Mandelbrot. *The fractal geometry of nature*. W. H. Freeman and Comp., New York, 3 edition, 1983.

[52] Roland Lenormand and Cesar Zarcone. Capillary fingering: Percolation and fractal dimension. *Transport in Porous Media*, 4(6):599-612, December 1989.

[53] A. B. D. Cassie and S. Baxter. Wettability of porous surfaces. *Trans. Faraday Soc.*, 40:546-551, 1944.

[54] Bo He, Junghoon Lee, and Neelesh A. Patankar. Contact angle hysteresis on rough hydrophobic surfaces. *Colloids and Surfaces A: Physicochemical and Engineering Aspects*, 248(1):101-104, 2004.

[55] Hyundo Lee, Ankur Gupta, T. Alan Hatton, and Patrick S. Doyle. Creating isolated liquid compartments using photopatterned obstacles in microfluidics. *Phys. Rev. Applied*, 7:044013, Apr 2017.

[56] Aniket S. Ambekar, Sujoy Mondal, and Vivek V. Buwa. Pore-resolved volume-of-fluid simulations of two-phase flow in porous media: Pore-scale flow mechanisms and regime map. *Physics of Fluids*, 33(10):102119, 2021.