# Fluid Meniscus Algorithms for Dynamic Pore-Network Modeling of Immiscible Two-Phase Flow in Porous Media

^{1}Beijing Computational Science Research Center, Beijing, China^{2}PoreLab, Department of Physics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

We present in detail a set of algorithms for a dynamic pore-network model of immiscible two-phase flow in porous media to carry out fluid displacements in pores. The algorithms are universal for regular and irregular pore networks in two or three dimensions and can be applied to simulate both drainage displacements and steady-state flow. They execute the mixing of incoming fluids at the network nodes, then distribute them to the outgoing links and perform the coalescence of bubbles. Implementing these algorithms in a dynamic pore-network model, we reproduce some of the fundamental results of transient and steady-state two-phase flow in porous media. For drainage displacements, we show that the model can reproduce the flow patterns corresponding to viscous fingering, capillary fingering and stable displacement by varying the capillary number and viscosity ratio. For steady-state flow, we verify non-linear rheological properties and transition to linear Darcy behavior while increasing the flow rate. Finally we verify the relations between seepage velocities of two-phase flow in porous media considering both disordered regular networks and irregular networks reconstructed from real samples.

## 1 Introduction

Flow of multiple immiscible fluids inside a porous medium shows a range of complex characteristics during transient as well as in steady state [1, 2]. A number of factors, such as the capillary forces at the fluid menisci, viscosity contrast, wettability and geometry of the system, make the properties of multiphase flow very different compared to single phase flow. When a non-wetting fluid displaces a wetting fluid the flow is called drainage whereas the opposite is called imbibition. During drainage, a less-viscous fluid displacing a more-viscous fluid in a porous medium creates a variety of fingering patterns, whereas a more viscous fluid displacing a less-viscous one shows a stable displacement front [3, 4]. The fingering patterns show different properties depending on whether the flow is dominated by capillary or viscous forces and correspondingly they are named as the capillary and viscous fingerings respectively [5, 6]. The capillary fingering patterns appear during slow displacement process and are well described by invasion percolation [7, 8], whereas the viscous fingering patterns appear during fast displacement and can be modeled by diffusion limited aggregation (DLA) model [9, 10]. If both the fluids are continuously fed into the porous medium, the initial transients will eventually die out and the system enters to a steady state when the macroscopic properties fluctuate around a steady average. It has been discovered that under steady-state conditions the two-phase flow rate does not obey the Darcy type linear relation between the total flow rate and pressure drop in the capillary dominated regime. Rather, it was found to have a power law dependence on the total pressure drop [11–16].

Two-phase flow in porous media has extensively been studied using laboratory experiments [3–6, 17, 18], statistical models [7, 9, 19] and numerical simulations [4, 20, 21]. There is also significant theoretical developments to describe the steady-state properties [22–25]. In recent years, due to the vast advancements in computer power and in high resolution scanning techniques for pore-space reconstruction, the computational methods became the leading tools for studying the two-phase flow. There are different modeling approaches. The direct numerical simulations, such as the volume-of-fluid method [26] and the level-set method [27, 28], perform discretization of the pore space into smaller cells and solve the Navier-Stokes equation. A popular voxel-based simulation technique is the lattice Boltzmann method (LBM) which solves the Boltzmann transport equations for different species of lattice gases at the discretized pore space [29–31]. These techniques all provide detailed information on the flow propagation at pore scale and are useful where the actual shape of the pores matter. However, due to the discretization of the pores, these models become computationally expensive with increasing system size when studying up-scaling problems.

In this regard, a computationally efficient method that can treat much larger systems is the pore-network modeling [32, 33]. In this modeling technique, a porous matrix is modeled by a network of links (pore throats) that are connected at nodes (pore bodies). The actual shapes of the pore space are replaced with simplified geometries and the average flow properties for each pore is considered to model the flow of the system. Unlike the voxel-based methods that discretize the pore space in smaller grids, in pore-network simulations one pore is the smallest computational unit. This simplification looses detailed description of fluid arrangements inside a pore, however it allows the pore-network models to apply for large systems and thereby studying their statistical properties. There are two main ingredients in a pore-network model: 1) solving the local pressure drops at pores and 2) displacing fluids accordingly. Based on these two ingredients, there are two major groups of pore-network models, quasi-static models and dynamic models. The quasi-static models are intended for the flow that is dominated by capillary forces, so that the viscous forces can be neglected. The displacement of fluids in the quasi-static models are performed by filling a whole pore at a time with invasion percolation-type algorithms [5, 7, 34, 35] where the filling of pores are decided by capillary entry pressures or by determining the stability of a meniscus for a given contact angle [36, 37]. Quasi-static models can successfully describe the equilibrium properties of two-phase flow at capillary dominated regime [38–40], however they are unable to capture the dynamic effects from the interaction between viscous and capillary forces at higher flow rates. This interaction between the viscous and capillary forces at the pore scale are taken into account in the dynamic pore-network models where the fluids inside pores are displaced under both the viscous and capillary pressure drops [41–43]. The viscous pressure drops are calculated by solving flow equations for fully developed viscous flow inside the pores and the capillary pressure drops are obtained from local fluid configurations inside a pore. There are many factors that make the dynamic models computationally more complex to implement compared to the quasi-static models and efficient algorithms are therefore necessary. One such factor is the mixing of fluids at the nodes and distributing them to the neighboring links while conserving the volumes of each fluid.

In this article, we present in detail a set of algorithms to model the displacements of two fluids in the links of a pore network and to distribute them to the neighboring links after they pass through nodes. The configuration of two fluids inside the pores in these algorithms at any instant are marked explicitly by the positions of all the menisci between the fluids. These menisci are displaced in small steps under the instantaneous viscous and capillary pressure drops at the pores. We call the algorithms as meniscus-dynamics algorithms which execute the transport of the fluids in the network. This is different from other pore-network models, such as in [41], where the meniscus positions are only available implicitly from the volumetric saturation of the fluid elements inside a pore, or from the model in [42] where a meniscus is moved through a whole volume element at each time step. Explicit positions of all the menisci in this meniscus-tracking model provide more detailed fluid configuration at any time. This not only provides straightforward calculation of capillary pressures from the meniscus positions but also resolves different dynamical events such as the retraction of invasion fronts after a Haines jump [44–47].

This *meniscus-tracking* approach to model the fluid transport in pore-network model was first introduced by Aker et. al. in the late nineties [43] to study transient two-phase flow – i.e. drainage – in a pore network. Those algorithms were restricted to regular networks with open boundary conditions and could only simulate drainage displacements, that is, when a non-wetting fluid invades a porous medium filled with a wetting fluid. Volumes of these two fluids as they pass through the nodes were not conserved in that model. Since then, the model and the algorithms were in continuous development to apply for different flow types and for different pore-networks. By implementing bi-periodic boundary conditions, steady-state flow was simulated by Knudsen et al. [48], however, the rules to transport the menisci were based on certain events at the nodes and can only be applied to a regular square network in two dimensions. Different driving conditions used in experiments, such as the injection of two fluids through alternating inlets at constant fractional flow [11] were also not implemented there. Alternate injection of fluids for steady-state flow was implemented in [50] and recently steady-state flow in irregular pore networks was simulated in [15, 49]. However, all the previous studies were performed for specific flow types or boundary conditions, and description of any universal set of meniscus algorithms were lacking in those works. This is the first time we present a complete model with a general set of meniscus algorithms which can be applied to simulate both the drainage displacements and the steady-state flow in different network geometries and driving conditions. Moreover, we present them here with all technical details so that the reader may reproduce the model.

The meniscus algorithms we present here do not need any further modifications when altering the network connectivity or topology and can be applied to regular or irregular networks. Both the transient and steady-state two-phase flow under different driving conditions, such as constant external pressure drop or constant flow rate, can be simulated. Furthermore, the algorithms also do not depend on different boundary conditions, such the injection of fluids in a network with open boundary or the flow in a network with periodic boundary. As we will see, these algorithms are simple and straight-forward to implement, however they are capable to capture the essential statistical properties of immiscible two-phase flow. This we will show by implementing these algorithms in a dynamic pore network model and reproducing a few fundamental results of transient and steady-state flow. When the statistical properties are concerned, fine details of dynamics are generally dropped out. This makes it possible to model the fluid displacements with simplified meniscus dynamics rules, while still preserving the fundamental statistical properties. Moreover, a set of universal algorithms that do not depend on the network geometry or the flow conditions are highly significant for the development of any commercial tool for multiphase flow simulations.

The article is structured as follows. In *Section 2*, we will present the governing equations that describe the two-phase flow at the pore level. We will then describe how we solve the equations to find the local pressures at the nodes at any time step. In *Section 3*, we present in detail the meniscus-dynamics algorithms which transport the fluids inside the links and distribute them to the next connected links after they pass a node. The algorithms are divided into three functions, which will be presented in the subsections. In *Section 4*, we will describe how different boundary conditions can be implemented in this model. As the exact computational code for a complete simulation will vary according to the system setup and boundary conditions, we will present the model in an algorithmic way. If necessary, the reader may contact the corresponding author to obtain a simplified code for a specific simulation. In order to provide the validation of the model, we will present a few examples in *Section 5* where we will reproduce a number of fundamental results of two-phase flow using the model. Both the transient and steady-state flow will be simulated and the corresponding results will be presented in *Subsections 5.1 and5.2* respectively. Finally, we will draw our conclusions in *Section 6*.

## 2 Flow Equations

We consider immiscible flow of two incompressible Newtonian fluids through a network of pores, where one of the fluids is more wetting than the other with respect to the pore walls. We denote the less and more wetting fluids as the non-wetting (n) and wetting (w), respectively. For a given network, two dimensionless macroscopic parameters that characterize the dynamics of two-phase flow are the capillary number (Ca) and the viscosity ratio (*M*). The capillary number is a measure of the ratio of viscous to capillary forces in the system. These parameters are defined as,

where *Q* is the total flow rate, γ is the interfacial tension between two fluids and *A* is the cross-sectional pore area of the network. **A**, we show a crop of a Hele-Shaw cell filled with a monolayer of glass beads [50] that is widely used as a two-dimensional porous medium in laboratory experiments. Such a porous medium may be modeled by a two-dimensional network of pores as shown in **B** with disorder in the link radii. A porous medium in 3D, a sample of sand-pack [49, 51], is shown in **C** and the reconstructed network from the sample is shown in **D**. Different techniques for reconstructing the pore-network from the scanned images of a sample can be found in a wide range of literature [53–57].

**FIGURE 1**. An overview of the network representation of two and three dimensional porous media. Figure **(A)** shows a **(B)**, where the distribution in link radii and node positions can be tuned according to the system properties. In **(B)**, the dark gray circles represent the beads and white represents the pores. The links are separated by dotted lines and the intersection of two such lines defines the position of a node. One link is colored by light gray. We like to point out that, **(B)** is not an exact reconstruction of the image shown in **(A)**, rather it is a simplified illustration. In **(C)**, we show a three dimensional Micro CT image of sand-pack (sample F42A in [51]) and the corresponding reconstructed pore network is shown in **(D)**.

In network representation, a pore is typically consists of two wider pore bodies that are connected by a narrower pore throat as shown in Figure 2A. In our modeling approach, we assign the total pore space to a *composite* link and the nodes do not contain any volume. The centers of the pore bodies correspond to the positions of nodes. This introduces a variation in the cross-sectional area along the length of the link which is modeled by a simplified hourglass shape as shown in Figure 2B. The interfacial pressure (

where *j* and *k*th meniscus inside the link. Here γ is the surface tension and *θ* is the contact angle between the meniscus and the pore wall. If *a* and *b* across the link, the instantaneous local flow rate *a* to *b* is proportional to the difference between the viscous pressure drop

where the sum is over all the menisci inside the link *j*, taking into account the direction of the capillary forces. *G*) defined as the ratio between the effective cross-sectional area of the pore and the square of its circumference. The effective cross-sectional areas of each of the three pore parts, the two wider pore bodies at the ends and the narrow throat in the middle, are calculated from the relation *ρ* is the radius of the inscribed circle in that pore part [60]. The mobility contribution from each part of the pore is then obtained from the relation

where

**FIGURE 2**. **(A)** Schematics of one pore in a reconstructed network which consists of three pore parts, two wider pore bodies at the ends and a narrow pore throat in the middle. **(B)** where the wetting and non-wetting fluids are colored by white and gray respectively. The variation in the interfacial pressure (

In order to move the fluids through the pores, the local pressures at nodes are needed to be solved at each time step. With the Kirchhoff equations, the net volume flux (*i*,

where the subscript *h* runs over the links connected to the node *i*. This sum, along with Eqs. 2, 3, constructs a set of linear equations. We solve these equations with conjugate gradient solver [63] or the Cholesky factorization method [64] and the solutions of which at any time step provides the local node pressures *Q* with open or periodic boundary conditions. More details about implementing different boundary conditions will be presented in *Section 4*.

## 3 Meniscus-Dynamics Algorithms

After the fluid velocities inside all the links are known from the solution of Eq. 3, the next step is to displace the two fluids inside the links and to distribute them to the next links after they mix at the nodes. These are done by updating the meniscus positions between the two fluids at every time step. In our implementation, the full algorithm for the meniscus displacements is sub-divided into three intermediate functions, namely the (a) meniscus_move, (b) meniscus_create and (c) meniscus_merge, which are applied to each link and/or node at every time step. The meniscus_move function makes lateral displacements of the menisci in proportional to the fluid velocities in the respective links and calculates the amount of fluids that exit from each link to the connected node. The meniscus_create function then distribute these fluids from each node to the connected links with outward flow. The total pore space in our model is assigned to the composite links of the network where the competition between the viscous and capillary forces at any time step is taken care by Eq. 3. The nodes in our system do not contain any physical volume and this distribution of fluids from node to links is only an intermediate step. We therefore considered this distribution rule to be democratic, that is, it does not treat the two fluids differently. The volumes of the two fluids distributed to the links are therefore proportional to the ratio of the velocities in the outward links and to the ratio of the incoming volumes. These democratic rules are also possible to apply without any further change for a mixed wet system where wettability varies from pore to pore [65, 66], as all the changes in the contact angles can be accounted by Eq. 3 and no change in the meniscus algorithms are necessary. These two functions interpret the volume conservations of the two fluids together with the Kirchhoff laws. Finally, the third function meniscus_merge is used to merge bubbles inside any link and it limits the maximum number of menisci (*Section 3.3*.

In the following we present with all the technical details how these three functions can be implemented in any pore network. They are illustrated in Figure 3. We will verify a few identities in order to check whether they introduce any unphysical effect on fundamental flow properties.

**FIGURE 3**. The three intermediate steps for the meniscus dynamics, the meniscus_move, meniscus_create and meniscus_merge, at every time step are illustrated here. Here five links, numbered as *i*. However the above illustrations are also valid for any number of links. Wetting and non-wetting fluids inside the links are colored by white and dark gray respectively. All of the pore space are distributed to the links in this model and the area around the node *i*, colored by light gray, does not contain any real volume. **(A)** shows an example configuration of two fluids at a time step. Solution of the Kirchhoff equations as described in *Section 2* provides the local flow rates *j* are moved toward (in the incoming links) or away from (in the outgoing links) the node by a distance **(B)**. This is performed by the function meniscus_move. The total volume of the fluids (*i* is now to be placed at the beginning of the links 4 and 5, which will fill the volume **(C)**, and the other is to place the non-wetting fluid first and the wetting next as shown in **(D)**. We adopt **(C)** and **(D)** alternatively at each consecutive time steps. **(E)** shows the meniscus_merge function for the link 5 when there are four allowed menisci. There link 5 has total six menisci before merging and the two bubbles around the two nearest menisci marked by the arrows are merged. The merging process are illustrated in more detail in Figure 4.

### 3.1 meniscus_move

The pore network consists of *i* for a node (*j* for a link (*k*. Inside the code, one also needs to mark the type of the menisci, for example, whether it is a meniscus at the beginning of a non-wetting bubble or at the end. This is necessary in order to calculate the direction of the capillary forces and also to measure the amount of two individual fluids that enter from a link to a node. Whether it is a regular network with constant coordination number for each node or an irregular network where the nodes have different coordination numbers, such as a reconstructed network from a real sample, it does not change any of the rules described in the following. In Figure 3, we illustrate them for an example of five links connected to a node *i*. At a time step, when all the pressures (*j*. The time step is then calculated as

where

After deciding the time step, menisci inside each link *j* are to be moved in the direction of

as shown in Figures 3A,B. Here, the links have uniform cross-sectional area in terms of the mobility, so all the menisci inside the same link move by the same distance. The positions of all the menisci inside any link (*incoming links* and the *outgoing links* for that node. In the example shown in Figure 3 the links 1, 2 and 3 are the incoming links and the links 4 and 5 are the outgoing links for the node *i* at that time step. Due to the displacement of the menisci by an amount *i* receives a certain volume (*j*, given by *j*. The total volume of fluids (*i* from all of its incoming links is therefore,

where *i*. Here *i*. We point out that, in order to apply this meniscus_move function for the whole network, one can loop through all links *j*, move the menisci by the distance

### 3.2 meniscus_create

After the function meniscus_move is performed over all the links, the list of the volumes *i* to an outgoing link *j* is given by, *i*. However, we still need to find out the volumes of each individual fluid for any outgoing link, that is, how much of the wetting (*j*. As illustrated in *Section 3*, we adopt a “democratic” rule to calculate this, that means both the wetting and non-wetting fluids get the same preference and the volumes depend on the flow rates *j*,

which imply that the ratio of the volumes of a fluid among all the outgoing links are the same as the ratio between flow rates among those links. Distributing the fluids in this way also preserves the volume conservation of each individual fluid, that is,

In each link, the wetting and non-wetting bubbles can be placed in two different ways, the non-wetting fluid at the beginning and the wetting fluid at the next, as shown in Figure 3C or in the other way as shown in Figure 3D. Here we adopt **C** or **D** alternatively at every consecutive time steps. One can also chose **C** or **D** randomly at every time step. This is equivalent of assuming that at some time step the wetting fluid coming from the incoming nodes pass the node before the non-wetting fluid enters the node, and the situation is the opposite in the other case.

As we mentioned before, the competition between the viscous and capillary forces only appear in Eq. 3 in our model. The meniscus displacement algorithms are therefore democratic and do not treat the wetting and the non-wetting fluids differently. This means, when the surface tension is set to zero, the capillary forces disappear and one will obtain [69],

where

where, *Section 5*.

### 3.3 meniscus_merge

After the two functions meniscus_move and meniscus_create are executed for each link, we look for any link in which the number of menisci exceeds the maximum number **A** (merge_back), where we identified two nearest non-wetting bubbles inside a link. Among these two non-wetting bubbles, only one, say the one in the front toward the flow direction, is moved back toward the other bubble and then merged. The length of this merged non-wetting bubble is the sum of the two bubbles before merging which maintains the volume conservation. This reduces the meniscus count by two inside the link and the meniscus positions are updated accordingly. This is illustrated in Figure 4A. However, when we measure the non-wetting fractional flow for zero surface tension with equal viscosities of the fluids, that is for two identical fluids, we find a deviation from Eq. 12. This is shown in Figure 5A where a small but systematic deviation from the diagonal **B** (merge_cm), where instead of merging two nearest non-wetting bubbles, we identified any two nearest menisci and merge the two bubbles across them. With this process, both the non-wetting or wetting bubbles can be merged whichever are the nearest. Moreover, instead of moving only one bubble toward the other, here we moved both the bubbles toward each other by such a distance so that the center of mass of these two bubbles does not change after merging. We illustrate this in Figure 4B. This rule shows satisfactory **B**), we see that the curve deviates from a typical S-shape and approaches to the diagonal line at higher saturation values (**B**. It seems that, as we moved both the bubbles toward each other while merging, small bubbles connected to different links at the nodes got disconnected and started flowing. Moreover, the results also seem to be sensitive on the meniscus count. We therefore updated the merging scheme further on. In the third and final rule **C** (merge_cmnn), we added one additional criteria compared to **B**. There we made sure that any fluid bubble that is in contact to a node, does not get disconnected from the node during the merging process. In order to do so, if one of the two nearest bubbles are connected to a node, we did not move that bubble during the merging process and only moved the other one. Everything else in this rule **C** are the same as rule **B**. With this merging scheme, we found exact match of fractional flow with Eq. 12 at zero surface tension and also obtained the expected qualitative behavior for non-zero capillary pressure. These are shown in Figure 5C and in the bottom row of Figure 6 respectively. Moreover, while changing the maximum bubble count with this merging scheme, we find no noticeable difference in the qualitative behavior of

**FIGURE 4**. Illustration of the meniscus_merge function for maximum meniscus count **(A)** we show the simplest possible try, namely the merge_back scheme. The link had six menisci before merging among which the menisci 2 and 3 are the nearest. We merge them by moving only one bubble toward the other, so the meniscus 1 does not change its position, only the meniscus 4 is moved toward 1 by a distance that is equal to the distance between 2 and 3. In **(B)** we show the next possible try, namely the merge_cm. Here instead of moving the meniscus 4 toward 1, we move both 1 and 4 toward each other by such distances so that the center of mass of the two bubbles do not change after merging. The third and final meniscus merging scheme is merge_cmnn, which is same as **(B)** with an additional criteria as illustrated in **(C)**. In the situation when one of the two nearest bubbles is connected to a node, such as the non-wetting bubble connected to the meniscus 5 in **(C)**, it will get disconnected from the node according to the merging scheme **(B)**. We avoid this disconnection of any bubble from a node by moving only the other bubble in such a situation as shown in **(C)**. So to merge the nearest menisci 4 and 5, we only move the meniscus 3 toward the node. The merge_cmnn scheme is therefore a combination of **(B)** and **(C)** which we adopt here.

**FIGURE 5**. Plot of steady-state non-wetting fractional flow *F*_{n} as a function of non-wetting saturation *S*_{n} for zero capillary pressure at the menisci. The three plots correspond to the three merging schemes **(A)** merge_back, **(B)** merge_cm and **(C)** merge_cmnn as described in Section 3.3. A small but systematic deviation from the diagonal *F*_{n} = *S*_{n} line can be observed for **A**.

**FIGURE 6**. Steady-state fractional flow at finite capillary numbers **B**: merge_cm (top row) and **C**: merge_cmnn (bottom row). The left and right plots for each capillary number correspond to **B**,

**FIGURE 7**. We show the effect of changing the maximum number of menisci (*N*_{max}) with the merging scheme **C**, the merge_cmnn. The steady‐state non-wetting fractional flow is plotted as a function of non-wetting saturation for *N*_{max}= 2, 3, 4 and 5. With such a change in the meniscus count, no noticeable difference is observed in the qualitative behavior in *F*_{n}.

## 4 Boundary Conditions

Simulations of different types of flow need proper boundary conditions to be implemented. The drainage displacement simulations can be performed with open boundary conditions (OB) where two opposite edges can be used as the *inlets* and *outlets* and the other edges are kept closed. Depending on the setup, all or some of the nodes and links at the inlet edge can be considered to inject a fluid. Depending on whether the system is driven under constant pressure drop or constant flow rate, either the node pressures (

The open boundary conditions can also be used for steady-state flow for the setup that is generally used in laboratory experiments [12]. There, instead of injecting one fluid, two fluids are injected simultaneously through alternate inlets. For this setup, all we need to do is to set the inlet node fluxes accordingly, all the *Q*) and the fractional flow

Another way to simulate the steady-state flow is by implementing periodic boundary (PB) conditions in the direction of overall flow, which is generally not possible in the experiments. In the open boundary setup, the injection of two fluids via alternate inlets in the experiments creates boundary effects and makes spatially homogeneous steady-state regime smaller. With the periodic boundaries, such boundary effects can be avoided and relatively smaller systems may be used for similar statistics. For this, the inlet and outlet edges are connected so that the fluids leaving at the outlets, enter the system again through the inlets. The links that connect the inlets and outlets are called as the *ghost links*. The network therefore becomes a closed system and the two fluids keep flowing through the system. The saturation **B**. For the reconstructed network in 3D, the opposite inlet and outlet edges are not identical. Therefore, to apply to periodic boundary conditions to this network we double the system by a mirror copy of the network in the direction of the flow, as shown by the two cuboids in **D**. This makes the inlets and outlets at the opposite edges identical, which we then connect with the ghost links to make the system periodic. With the periodic boundary, the inlets are seen as the neighbors of the outlets for the meniscus-dynamics algorithms. Here, the global pressure drop

**FIGURE 8**. Implementation of periodic boundary conditions for 2D and 3D networks. A 2D regular network is shown in **(A)** where periodic boundary conditions are applied in both the directions. Here the overall flow is in the upward direction and the gray links at the top row represent the ghost links. The links are drawn as cylindrical tubes for the simplicity of drawing, but they are hourglass shaped in terms of the capillary pressure. The nodes are marked by small dots where the red dots at the left and right edges represent identical nodes. This makes the flow essentially on the surface of a torus as shown in **(B)** where the arrow represents the effective flow direction. For 3D, a network reconstructed from Beria sandstone is shown in **(C)** where the overall flow is from left to right. Periodic boundary conditions are applied in the direction of flow by adding a mirror copy of the the network as shown in **(D)**. The ghost links at the right are colored by gray.

The structure of the complete simulation with all the meniscus dynamics functions is the following:

(1) Network: construct or read

(2) Define: boundary conditions

(3) Initialize: random or sequential fluid distribution

(4) **for****do**

(5) Solve the pressure field

(6) **for****do**

(7) meniscus_move(*j*)

(8) **end for**

(9) **for****do**

(10) meniscus_create(*j*)

(11) meniscus_merge(*j*)

(12) **end for**

(13) Calculate measurable quantities at *t*

(14) **end for**

During initialization, the initial positions of all the menisci in each link need to be defined depending on the saturation and how we want to start the simulation. Then after solving the pressure field, the meniscus_move function is to be performed on all the links which will generate the array of the node fluxes *i*. These arrays will then be used as the input to the meniscus_create function which will create the new menisci in the links. The code for the first three steps in the above method will vary according to different requirements of the simulation.

## 5 Applications and Validation

As stated earlier, the meniscus algorithms in this pore-network model has the flexibility to be applied for different network topologies and boundary conditions. Moreover, we can study both the transient and the steady-state properties. In our simulations, we consider a network of links forming a tilted square lattice at an angle

### 5.1 Drainage Displacements

As described in the introduction, when a non-wetting fluid is invaded into a porous medium filled with a wetting fluid, it generates different types of invading flow patterns depending on the capillary number and viscosity ratio. A less viscous fluid displacing a high viscous fluid creates viscous fingering patterns at a high Ca [6, 17] and capillary fingering patterns at low Ca [8]. The capillary and viscous fingering patterns resemble with the invasion percolation [7] and diffusion limited aggregation (DLA) models [9] respectively. Alternatively, when a high viscous fluid displaces a low viscous fluid, a stable displacement front is observed. In order to verify whether our pore-network model can reproduce these different flow patterns during the drainage, we run different simulations for different values of the capillary number *M*. In our simulations, we set *Q* and the viscosities of the two fluids to set the values of Ca and *M*, and no modification in the meniscus dynamics algorithms were made between different simulations.

**FIGURE 9**. Development of transient flow patterns during the drainage simulations in a disordered square network of 64 × 64 links. The network is initially filled with wetting fluid (gray) and the non-wetting fluid (blue) is injected through four inlet nodes at the bottom edge of the network, shown by the arrow, with a constant flow rate. The top edge of the network is kept open which works as the outlet. Periodic boundary conditions are applied in the horizontal direction, therefore the left and right edges are connected together. Only the capillary number Ca and the viscosity ratio M are altered during these simulations and all other flow parameters and the algorithms for meniscus dynamics are exactly the same. The flow patterns show the different regimes of transient two-phase flow, namely **(A)** the viscous fingering **(B)** the stable displacement and **(C)** the capillary fingering. We like to point out that, though we adopted “democratic” rules for meniscus algorithms, it is the solution of the flow equations that lead the fluids to generate these different flow patterns.

### 5.2 Steady State

Steady-state flow can be simulated by implementing either open or periodic boundary conditions as shown in Figure 10. In the open boundary simulation shown in the first row, the wetting and non-wetting fluids are injected through alternate injection points at the bottom edge of the network. The top edge is kept open through which fluids leave the system. The two edges on the sides are connected through periodic boundary which can also be kept closed if necessary. The simulations are performed at constant fractional flow by setting the flow rates at the inlet links. Due to the traces of the injections near the inlet edge, a long enough system is necessary for the open boundary simulations to obtain a region of spatially homogeneous steady-state flow away from the inlets. The second row of Figure 10 shows simulation with periodic boundary conditions in both directions. Here the system is closed and the control parameter is the fluid saturation. In the third row, we show the simulation with a three dimensional reconstructed Berea network where the periodic boundary condition is applied in the direction of the overall flow. The total flow rate is controlled in these simulations and we measure the global pressure drop (

**FIGURE 10**. Steady-state simulations at

In the following, we present two sets of simulation results with this dynamic pore network model to show the agreement with the steady-state properties of two-phase flow. First, we will present the effective rheological properties where the total flow rate shows non-linear dependence on the pressure drop in the capillary dominated regime [12–15]. Next, we will measure seepage velocities and will verify the relations between them [25]. We will also describe in detail how to measure the different quantities from the simulation, such as the flow rates and the seepage velocities of different fluid components.

##### Effective Rheology and Crossover from Linear to Non‐linear Flow Regime

When two Newtonian fluids flow through a porous media, it was found experimentally [11–13, 15], numerically [14, 15] and with mean field calculations [14] that, in the regime when capillary forces compete with the viscous forces, that is, in the low capillary number regime, the total flow rate *Q* in the steady state does not follow the linear Darcy relation and varies quadratically with the excess pressure drop. This can be written as,

where (*Q* to increase faster than the increase in

In order to verify this non-linear rheological behavior with our model, we performed steady-state simulations with the 2D square network and the 3D Berea network. The details of the networks were described in *Section 5*. Simulations are performed at constant total flow rate *Q* and the average global pressure drop *Q* as a function of (*y*-axis correspond to the values of

**FIGURE 11**. Variation of the total volumetric flow rate *Q* (*y*-axis correspond to the values of threshold pressures (

##### Relation between Seepage Velocities

We will now measure the seepage velocities of the fluids in the steady state with this model and will verify the relations between them. When the system is driven under a constant pressure drop *v* and the fluid saturations can be derived using the Euler homogeneity property of the total flow rate *Q* in the steady state [25]. These relations necessitate a new velocity function, namely the co-moving velocity (

respectively, where *x* direction. The length of the systems in this direction is *L*. A cross section normal to the *x* direction is shown by an orange straight line for the 2D network and by an orange grid plane for the 3D network. Orthogonal views of these cross sections are shown underneath where the gray and blue patches show the pore areas occupied by the wetting and non-wetting fluids respectively. The sum of the areas of individual colors correspond to the wetting and non-wetting pore areas *x*. The total pore area *A* related to all the fluids is therefore given by *Q* of the two fluids is the sum of the wetting and non-wetting flow rates given by *v* associated with the total flow rate *Q* is defined by,

and we can find,

by using the relations mentioned above.

**FIGURE 12**. Description of the system to measure the flow rates (*Q*, *A*, *v*, *x* direction which is the direction of overall flow. A random cross section of the system normal to the direction of overall flow is shown by the orange line for 2D and by the orange grid plane for 3D. The normal view of the cross sections are shown underneath where the gray and blue patches show the occupation by wetting and non-wetting fluids in the cross section. The total gray and blue areas correspond to the wetting and non-wetting pore areas respectively and the sum of them correspond to the total pore area along this cross section. The averages of these areas over all the possible cross sections lead to the measurement of *A*,

In our network simulations, we have the information about the local flow rates *i* at any time step. When the system is driven under constant pressure drop (

where *L*, i. e. 64 here. Similarly, we can calculate the cross-sectional areas as,

where *j*, projected into the plane normal to the flow direction. The wetting saturations of the links

For an irregular network that is considered here in 3D, the links are of different lengths and oriented in different directions. In that case, measurement of the flow rates and the areas by summing over all the links and dividing by the number of rows using Eqs. 17, 18 will lead to wrong results. In this case, we measure these quantities in the following way. Let us consider an orthogonal cross-sectional plane at a random position through which we like to measure the flow rates (see Figure 12). The probability that any link *j* will pass through this plane will be proportional to *j* in the *x* direction, the direction of the overall flow. After the link is selected, the probability that the plane will pass through the wetting fluid inside the link will be proportional to the local wetting saturation *Q*,

where *L* is the length of the network in the *x* direction. The areas can be measured in the similar way given by,

For the regular 2D network, *v*,

The calculation of *v*, *v* where we used Eqs. 17, 18 for 2D and Eqs. 19, 20 for 3D. The plots show exact match of Eq. 16 for the whole range of parameters.

**FIGURE 13**. Plots of *v* for 2D and 3D are shown in the top. The flow rates and pore areas to calculate *v* are measured using Eqs. 17, 18 for 2D and using **19** and **20** for 3D. The velocities are in the unit of

The total flow rate *Q* in the steady state is a homogeneous function of order one of the pore areas *Q* scales as, *Q* leads to a new set of equations between the seepage velocities. Complete derivations of the equations can be found in Ref. [25] and here we will present them in brief and will use them to validate our model. Taking the derivative of the homogeneity equation of *Q* with respect to λ and then setting

These two partial derivatives in the above equation have the units of velocity and correspondingly they define two thermodynamic velocities

With these definitions Eq. 21 becomes,

which has the similar form of Eq. 16. However, this does *not* imply that the thermodynamic velocities

which fulfill both Eqs. 16, 23. This velocity function

and

In order to verify whether our pore-network model with the set of meniscus algorithms described here do satisfy these equations, we perform a large number of simulations with a wide range of parameters for the 2D square network and the 3D Berea network. Five viscosity ratios, *M*. For each value of *M*, three different values of the surface tension, *M* and γ, we considered a set of pressure drops

The seepage velocities are measured for all the simulations and the derivatives with respect to the saturations are measured with central difference technique. We calculate the values of

Eqs. 16, 26 transform the wetting and non-wetting velocities (*v*, *v*,

With these equations we can verify the measured values of wetting and non-wetting seepage velocities

**FIGURE 14**. Plot of the wetting and non-wetting seepage velocities (mm/sec) calculated using Eq. 27 against the measured values of

Finally, we plot the co-moving velocity

By fitting all the data points for the whole set of simulations, we find *v* and the variation of *v* while changing the saturation. Therefore it is not enough to specify only the *v* to determine *a*, *b* and *c* characterize this functional dependency. However, we have considered only a few samples in our simulations to explore the whole parameter space which lead to higher statistical errors and therefore more sample averages may be necessary to improve the statistics.

**FIGURE 15**. Plot of the co-moving velocity

## 6 Summary

We presented a detailed description of a set of algorithms for transporting fluids in a dynamic pore-network model of two-phase flow in porous media. The displacements of the fluids in this model are monitored by updating the positions of all the menisci with time. The fundamental concept of the algorithms is simple, at every time step all the menisci are displaced according to the velocity in the corresponding link and then all the fluids arriving at a node from the incoming links are distributed to the outgoing links. There, the distributed volumes are proportional to the ratio between the velocities in the outgoing links and the ratio between the incoming fluid volumes. We presented these algorithms with all the technical details so that it is possible for the reader to reproduce this model. We have illustrated that this pore-network model and the meniscus algorithms are applicable to both the regular and irregular network topologies in two and three dimensions. By reproducing some of the fundamental results of two-phase flow, we have shown that the model can be used to simulate both the transient and steady-state flow. The model is able to generate different drainage displacement patterns by altering the capillary number and viscosity ratios. In steady state, the model successfully reproduces the linear to non-linear transition in the effective rheological properties as well as the relations between the seepage velocities. For all the results presented here, we used the same meniscus dynamics algorithms without any modification.

Apart from the standard applications of two-phase flow, this model can also be very useful to study the effect of the pore-scale dynamics at low capillary numbers, such as the retractions of menisci after a Haines jump and the long-range capillary effects, which are not possible to capture with the quasi-static or percolation-based models as observed experimentally in [46]. There is also a number of possibilities for further developments. In a recent paper by Zhao et al. [4] ten different groups with different approaches to model two-phase flow in porous media were invited to reproduce fluid injection in a circular Hele-Shaw cell at different capillary numbers and wetting properties, ranging from drainage to strong imbibition, i.e., imbibition where film flow dominates the process. The conclusion of that work was, whereas all the different approaches were able to reproduce the drainage processes well, none succeeded in reproducing strong imbibition. An important mechanism during imbibition is wetting film flow along pore walls. Our presented work accounts the Darcy-like creep flow and modeling the film flow mechanisms within the same framework can also be done as demonstrated in the first attempt in [73]. Another possible development can be replacing the time integration with a Monte Carlo algorithm as proposed by Savani et al. [74] that can further improve the computational efficiency of this model. However, that method needs to be tested at low capillary numbers and was only been implemented for regular lattices. Other important phenomena such as wetting angle hysteresis is straight-forward to implement. Further applications such as the wettability alteration due to changes in the composition of the immiscible fluids can also be investigated with this model as demonstrated in [65, 66].

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

The meniscus algorithms in their latest form was developed by SS in collaboration with MG, MV, and AH. SS performed all the numerical simulations and the results were analyzed by all the authors. The analytical part related to the relation between the seepage velocities was developed by AH. SS wrote the first draft of the manuscript which then updated to its final form by all the authors.

## Funding

SS was supported by the National Natural Science Foundation of China under grant number 11750110430. This work was partly supported by the Research Council of Norway through its Center of Excellence funding scheme, project number 262644.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We thank our co-workers, past and present, who have been involved in developing this model through its various stages, Eyvind Aker, G. George Batrouni, Morten Grøva, Henning A. Knudsen, Knut Jørgen Måløy, Pål Eric Øren, Thomas Ramstad, Subhadeep Roy, Isha Savani and Glen Tørå.

## References

1. Dullien FAL. *Porous media: fluid transport and pore structure*. San Diego, CA: Academic Press (1992)

3. Løvoll G, Méheust Y, Toussaint R, Schmittbuhl J, Måløy KJ. Growth activity during fingering in a porous Hele-Shaw cell. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2004) 70:026301. doi:10.1103/PhysRevE.70.026301

4. Zhao B, MacMinn CW, Primkulov BK, Y Chen, Valocchi AJ, Zhao J, et al. Comprehensive comparison of pore-scale models for multiphase flow in porous media. *Proc Natl Acad Sci USA* (2019) 116:13799. doi:10.1073/pnas.1901619116

5. Lenormand R, Touboul E, Zarcone C. Numerical models and experiments on immiscible displacements in porous media. *J Fluid Mech* (1988) 189:165. doi:10.1017/S0022112088000953

6. Måløy KJ, Feder J, Jøssang T. Viscous fingering fractals in porous media. *Phys Rev Lett* (1985) 55:2688. doi:10.1103/PhysRevLett.55.2688

7. Wilkinson D, Willemsen JF. Invasion percolation: a new form of percolation theory. *J Phys Math Gen* (1983) 16:3365. doi:10.1088/0305-4470/16/14/028

8. Lenormand R, Zarcone C. Invasion percolation in an etched network: measurement of a fractal dimension, *Phys Rev Lett* (1985) 54:2226. doi:10.1103/PhysRevLett.54.2226

9. Witten TA, Sander LM. Diffusion-limited aggregation, a kinetic critical phenomenon. *Phys Rev E* (1981) 47:1400. doi:10.1103/PhysRevLett.47.1400

10. Paterson L. Diffusion-limited aggregation and two-fluid displacements in porous media. *Phys Rev Lett* (1984) 52:1621. doi:10.1103/PhysRevLett.52.1621

11. Tallakstad KT, Knudsen HA, Ramstad T, Løvoll G, Måløy KJ, Toussaint R, Flekkøy EG, et al. Steady-state two-phase flow in porous media: statistics and transport properties. *Phys Rev Lett* (2009) 102:074502. doi:10.1103/PhysRevLett.102.074502

12. Tallakstad KT, Løvoll G, Knudsen HA, Ramstad T, Flekkøy EG, Måløy KJ. Steady-state, simultaneous two-phase flow in porous media: an experimental study. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2009) 80:036308. doi:10.1103/PhysRevE.80.036308

13. Rassi EM, Codd SL, Seymour JD. Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow. *New J Phys* (2011) 13:015007. doi:10.1088/1367-2630/13/1/015007

14. Sinha S, Hansen A. Effective rheology of immiscible two-phase flow in porous media. *Europhys Lett* (2012) 99:44004. doi:10.1209/0295-5075/99/44004

15. Sinha S, Bender AT, Danczyk M, Keepseagle K, Prather CA, Bray JM, et al. Effective rheology of two-phase flow in three-dimensional porous media: experiment and simulation. *Transport Porous Media* (2017) 119:77. doi:10.1007/s11242-017-0874-4

16. Gao Y, Lin Q, Bijeljic B, Blunt MJ. Pore-scale dynamics and the multiphase Darcy law. *Phys Rev Fluids* (2020) 5:013801. doi:10.1103/physrevfluids.5.013801

17. Chen JD, Wilkinson D. Pore-scale viscous fingering in porous media, *Phys Rev Lett* (1985) 55:1892. doi:10.1103/PhysRevLett.55.1892

18. Saffman PG, Taylor GI. The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. *Proc. Royal Soc. A: Math. Phys. Eng. Sc.* (1958) 245:312. doi:10.1098/rspa.1958.0085

19. Bensimon D, Kadanoff LP, Liang S, Shraiman BI, Tang C. Viscous flows in two dimensions. *Rev Mod Phys* (1986) 58:977. doi:10.1103/RevModPhys.58.977

20. Koplik J, Lasseter TJ. Two-phase flow in random network models of porous media. *Soc Petrol Eng J* (1985) 25:89. doi:10.2118/11014-PA

21. Gjennestad MA, Winkler M, Hansen A. Pore network modeling of the effects of viscosity ratio and pressure gradient on steady-state incompressible two-phase flow in porous media. *Transp. Porous Med* (2020) 132 (2):355–379. doi:10.1007/s11242-020-01395-z

22. Hassanizadeh SM, Gray WG. Toward an improved description of the physics of two-phase flow. *Adv Water Resour* (1993) 16:53. doi:10.1016/0309-1708(93)90029-F

23. Hilfer R. Macroscopic equations of motion for two-phase flow in porous media. *Phys Rev E* (1998) 58:2090. doi:10.1103/PhysRevE.58.2090

24. Gray WG, Miller CT. *Introduction to the thermodynamically constrained averaging theory for porous medium systems*. Berlin, Germany: Springer (2014)

25. Hansen A, Sinha S, Bedeaux D, Kjelstrup S, Gjennestad MA, Vassvik M. Relations between seepage velocities in immiscible, incompressible two-phase flow in porous media. *Transport Porous Media* (2018) 125:565. doi:10.1007/s11242-018-1139-6

26. Raeini AQ, Blunt MJ, Bijeljic B. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. *J Comput Phys* (2012) 231:5653. doi:10.1016/j.jcp.2012.04.011

27. Jettestuen E, Helland JO, Prodanovic M. A level set method for simulating capillary-controlled displacements at the pore scale with nonzero contact angles. *Water Resour Res* (2013) 49:4645. doi:10.1002/wrcr.20334

28. Gjennestad MA, Munkejord ST. Modelling of heat transport in two-phase flow and of mass transfer between phases using the level-set method. *Energy Proc* (2015) 64:53. doi:10.1016/j.egypro.2015.01.008

29. Gunstensen AK, Rothman DH, Zaleski S, Zanetti G. Lattice Boltzmann model of immiscible fluids. *Phys Rev A* (1991) 43:4320. doi:10.1103/PhysRevA.43.4320

30. Ramstad T, Øren PE, Bakke S. Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method. *SPE J* (2010) 15:917. doi:10.2118/124617-PA

31. Aursjø O, Løvoll G, Knudsen HA, Flekkøy EG, Måløy KJ. A direct comparison between a slow pore scale drainage experiment and a 2D lattice Boltzmann simulation. *Transport Porous Media* (2011) 86:125. doi:10.1007/s11242-010-9611-y

32. Blunt MJ. Flow in porous media-pore-network models and multiphase flow. *Curr Opin Colloid Interface Sci* (2001) 6:197. doi:10.1016/S1359-0294(01)00084-X

33. Joekar-Niasar V, Hassanizadeh SM. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: a review. *Crit Rev Environ Sci Technol* (2012) 42:1895. doi:10.1080/10643389.2011.574101

34. Chandler R, Koplik J, Lerman K, Willemsen JF. Capillary displacement and percolation in porous media. *J Fluid Mech* (1982) 119:249. doi:10.1017/S0022112082001335

35. Blunt MJ. Physically-based network modeling of multiphase flow in intermediate-wet porous media. *J Petrol Sci Eng* (1999) 20:117. doi:10.1016/S0920-4105(98)00010-2

36. Cieplak M, Robbins MO. Dynamical transition in quasistatic fluid invasion in porous media. *Phys Rev Lett* (1988) 60:2042. doi:10.1103/PhysRevLett.60.2042

37. Cieplak M, Robbins MO. Influence of contact angle on quasistatic fluid invasion of porous media. *Phys Rev B Condens Matter* (1990) 41:11508. doi:10.1103/PhysRevB.41.11508

38. Øren P-E, Bakke S, Arntzen OJ. Extending predictive capabilities to network models. *SPE J* (1998) 3:324. doi:10.2118/52052-PA

39. Blunt MJ, Jackson MD, Piri M, Valvatne PH. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. *Adv Water Resour* (2002) 25:1069. doi:10.1016/S0309-1708(02)00049-0

40. Primkulov BK, Talman S, Khaleghi K, Shokri AR, Chalaturnyk R, Zhao B, et al. Quasistatic fluid-fluid displacement in porous media: invasion-percolation through a wetting transition. *Phys. Rev. Fluids* (2018) 3:104001. doi:10.1103/PhysRevFluids.3.104001

41. Joekar-Niasar V, Hassanizadeh SM, Dahle H. Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling. *J Fluid Mech* (2010) 655:38. doi:10.1017/S0022112010000704

42. Hammond PS, Unsal E. A dynamic pore network model for oil displacement by wettability-altering surfactant solution. *Transport Porous Media* (2012) 92:789. doi:10.1007/s11242-011-9933-4

43. Aker E, Måløy KJ, Hansen A, Batrouni GG. A two-dimensional network simulator for two-phase flow in porous media. *Transport Porous Media* (1998) 32:163. doi:10.1023/A:1006510106194

44. Haines WB. Studies in the physical properties of soil. v. the hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith. *J Agric Sci* (1930) 20:97. doi:10.1017/S002185960008864X

45. Berg S, Ott H, Klapp SA, Schwing A, Neiteler R, Brussee N Real-time 3D imaging of Haines jumps in porous media flow. *Proc Natl Acad Sci USA* (2013) 110:3755. doi:10.1073/pnas.1221373110

46. Armstrong RT, Berg S. Interfacial velocities and capillary pressure gradients during haines jumps. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2013) 88:043010. doi:10.1103/PhysRevE.88.043010

47. Måløy KJ, Furuberg L, Feder J, Jøssang T. Dynamics of slow drainage in porous media. *Phys Rev Lett* (1992) 68:2161. doi:10.1103/PhysRevLett.68.2161

48. Knudsen HA, Aker E, Hansen A. Bulk flow regimes and fractional flow in 2D porous media by numerical simulations. *Trannsp. Porous Med* (2002) 47:99. doi:10.1023/A:1015039503551

49. Ramstad T, Hansen A, Øren PE. Flux-dependent percolation transition in immiscible two-phase flows in porous media. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2009) 79:036310. doi:10.1103/PhysRevE.79.036310

50. Erpelding M, Sinha S, Tallakstad KT, Hansen A, Flekkøy EG, Måløy KJ. History independence of steady state in simultaneous two-phase flow through two-dimensional porous media. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2013) 88:053004. doi:10.1103/PhysRevE.88.053004

51.https://www.imperial.ac.uk/earth-science/research/research-groups/perm/research/pore-scale-modelling/micro-ct-images-and-networks/

52. Vogel HJ, Roth K. Quantitative morphology and network representation of soil pore structure. *Adv Water Resour* (2001) 24:233. doi:10.1016/S0309-1708(00)00055-5

53. Okabe H, Blunt MJ. Pore space reconstruction using multiple-point statistics. *J Petrol Sci Eng* (2005) 46:121. doi:10.1016/j.petrol.2004.08.002

54. Øren P-E, Bakke S. Process based reconstruction of sandstones and prediction of transport properties. *Transport Porous Media* (2002) 46:311. doi:10.1023/A:1015031122338

55. Øren P-E, Bakke S. Reconstruction of berea sandstone and pore-scale modeling of wettability effects. *J Petrol Sci Eng* (2003) 39:177. doi:10.1016/S0920-4105(03)00062-7

56. Dong H. *Micro-CT imaging and pore network extraction*. [PhD thesis]. London (UK): Imperial College (2007).

57. Dong H, Blunt MJ. Pore-network extraction from micro-computerized-tomography images. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2009) 80:036307. doi:10.1103/PhysRevE.80.036307

58. Sinha S, Hansen A, Bedeaux D, Kjelstrup S. Effective rheology of bubbles moving in a capillary tube. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2013) 87:025001. doi:10.1103/PhysRevE.87.025001

60. Mason G, Morrow NR. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. *J Colloid Interface Sci* (1991) 141:262. doi:10.1016/0021-9797(91)90321-X

62. Jia P, Dong M, Dai L, Yao J. Slow viscous flow through arbitrary triangular tubes and its application in modelling porous media flows. *Transport Porous Media* (2008) 74:153. doi:10.1007/s11242-007-9187-3

63. Batrouni GG, Hansen A. Fourier acceleration of iterative processes in disordered systems. *J Stat Phys* (1988) 52:747. doi:10.1007/BF01019728

64. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. *Numerical recipes in C: The art of scientific computing*. 2nd ed. New York, NY: Cambridge University Press (1992)

65. Sinha S, Grøva M, Ødegården TB, Skjetne E, Hansen A. Local wettability reversal during steady-state two-phase flow in porous media. *Phys Rev E - Stat Nonlinear Soft Matter Phys* (2011) 84:037303. doi:10.1103/PhysRevE.84.037303

66. Flovik V, Sinha S, Hansen A. Dynamic wettability alteration in immiscible two-phase flow in porous media: effect on transport properties and critical slowing down. *Front Phys* (2015) 3:86. doi:10.3389/fphy.2015.00086

67. Garstecki P, Fuerstman MJ, Stone HA, Whitesides GM. Formation of droplets and bubbles in a microfluidic T-junction-scaling and mechanism of break-up. *Lab Chip* 6, 437 (2006) doi:10.1039/b510841a

68. Gjennestad MA, Vassvik M, Kjelstrup S, Hansen A. Stable and efficient time integration of a dynamic pore network model for two-phase flow in porous media. *Front Phys* (2018) 6:56. doi:10.3389/fphy.2018.00056

69. Sinha S, Gjennestad MA, Vassvik M, Winkler M, Hansen A, Flekkøy EG. Rheology of high-capillary number two-phase flow in porous media. *Front Phys* (2019) 7:65. doi:10.3389/fphy.2019.00065

70. Craig FF *The reservoir engineering aspects of waterflooding*. Richardson, TX: Monograph Series, Society of Petroleum Engineers (1971)

72. Roy S, Hansen A, Sinha S. Effective rheology of two-phase flow in a capillary fiber bundle model. *Front Phys* (2019) 7:92. doi:10.3389/fphy.2019.00092

73. Tørå G, Øren PE, Hansen A. A dynamic network model for two-phase flow in porous media. *Transport Porous Media* (2012) 92:145. doi:10.1007/s11242-011-9895-6

Keywords: pore-network modeling, two-phase flow, porous media, interface dynamics, numerical simualtion

Citation: Sinha S, Gjennestad MA, Vassvik M and Hansen A (2021) Fluid Meniscus Algorithms for Dynamic Pore-Network Modeling of Immiscible Two-Phase Flow in Porous Media. *Front. Phys.* **8**:548497. doi: 10.3389/fphy.2020.548497

Received: 02 April 2020; Accepted: 09 November 2020;

Published: 11 March 2021.

Edited by:

José S. Andrade Jr, Federal University of Ceara, BrazilReviewed by:

Tom Bultreys, Ghent University, BelgiumAdam Jan Gadomski, University of Science and Technology (UTP), Poland

Copyright © 2021 Sinha, Gjennestad, Vassvik and Hansen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Santanu Sinha, santanu@csrc.ac.cn