A Dynamic Network Simulator for Immiscible Two-Phase Flow in Porous Media

TL;DR

This paper introduces a comprehensive dynamic pore-network model with algorithms capable of simulating various two-phase flow patterns in porous media, capturing both transient and steady-state behaviors with different boundary conditions.

Contribution

The paper develops a general set of algorithms for dynamic pore-network modeling of two-phase flow, capable of reproducing complex flow patterns and rheological transitions in porous media.

Findings

Reproduces viscous, capillary fingering, and stable displacement patterns.

Verifies the transition from linear to non-linear rheological behavior.

Satisfies relations between seepage velocities in steady-state flow.

Abstract

We present in detail a set of algorithms to carry out fluid displacements in a dynamic pore-network model of immiscible two-phase flow in porous media. The algorithms are general and applicable to regular and irregular pore networks in two and three dimensions with different boundary conditions. Implementing these sets of algorithms, we describe a dynamic pore-network model and reproduce some of the fundamental properties of both the transient and steady-state two-phase flow. During drainage displacements, we show that the model can reproduce the flow patterns corresponding to viscous fingering, capillary fingering and stable displacement by altering the capillary number and the viscosity ratio. In steady-state flow, the model verifies the linear to non-linear transition of the effective rheological properties and satisfy the relations between the seepage velocities of two-phase flow in porous media.

Equations

$$\text{Ca}=\frac{{\mu }_{\mathrm{e}}Q}{\gamma A}\text{and}M=\frac{{\mu }_{\mathrm{n}}}{{\mu }_{\mathrm{w}}}$$

$$|p_c\left(x_k\right)|=\frac{2\gamma \cos \theta }{r_j}\left[1-\cos \left(\frac{2\pi x_k}{l_j}\right)\right]$$

$$q_j=-\frac{g_j}{l_j{\mu }_j}\left[\Delta p_j-\sum _k{p}_c(x_k)\right]$$

$$\frac{l_j}{g_j}=\frac{{\lambda }_1}{g_1}+\frac{{\lambda }_2}{g_2}+\frac{{\lambda }_3}{g_3}$$

$$f_i=\sum _h{q}_h=0,$$

$$\Delta t=\frac{0.1l_j}{c_{j,\max }},$$

$$\Delta x_j=\frac{q_j}{a_j}\Delta t$$

$$V_i=\sum _{j\in \mathrm{I}}{\varphi }_{i,j}=V_i^{\mathrm{w}}+V_i^{\mathrm{n}},$$

$${\varphi }_{i,j}^{\mathrm{w}}=q_j{V}_i^{\mathrm{w}}\Delta t/V_i\text{and}{\varphi }_{i,j}^{\mathrm{n}}=q_j{V}_i^{\mathrm{n}}\Delta t/V_i$$

$$\sum _{j\in \mathrm{O}}{\varphi }_{i,j}^{\mathrm{w}}=V_i^{\mathrm{w}}\text{and}\sum _{j\in \mathrm{O}}{\varphi }_{i,j}^{\mathrm{n}}=V_i^{\mathrm{n}}.$$

$$Q_{\mathrm{w}}(S_{\mathrm{w}},M)=Q_{\mathrm{n}}(1-S_{\mathrm{w}},1/M)$$

$$F_{\mathrm{w}}=S_{\mathrm{w}}$$

$$Q$$

$$\sim (\Delta P-P_t{)}^2$$

$$v_{\mathrm{w}}=\frac{Q_{\mathrm{w}}}{A_{\mathrm{w}}}$$

$$v=\frac{Q}{A}$$

$$v=S_{\mathrm{w}}{v}_{\mathrm{w}}+S_{\mathrm{n}}{v}_{\mathrm{n}}$$

$$Q=\frac{1}{N_L}\sum _j{q}_j,Q_{\mathrm{w}}=\frac{1}{N_L}\sum _j{s}_{\mathrm{w},j}{q}_j\mathrm{and}{Q}_{\mathrm{w}}=\frac{1}{N_L}\sum _j(1-s_{\mathrm{w},j})q_j$$

$$A=\frac{1}{N_L}\sum _j{a}_j,A_{\mathrm{w}}=\frac{1}{N_L}\sum _j{s}_{\mathrm{w},j}{a}_j\mathrm{and}{A}_{\mathrm{n}}=\frac{1}{N_L}\sum _j(1-s_{\mathrm{w},j})a_j$$

$$Q=\frac{1}{L}\sum _j{l}_{x,j}{q}_j,Q_{\mathrm{w}}=\frac{1}{L}\sum _j{l}_{x,j}{s}_{\mathrm{w},j}{q}_j\mathrm{and}{Q}_{\mathrm{n}}=\frac{1}{L}\sum _j{l}_{x,j}(1-s_{\mathrm{w},j})q_j$$

$$A=\frac{1}{L}\sum _j{l}_{x,j}{a}_j,A_{\mathrm{w}}=\frac{1}{L}\sum _j{l}_{x,j}{s}_{\mathrm{w},j}{a}_j\mathrm{and}{A}_{\mathrm{n}}=\frac{1}{L}\sum _j{l}_{x,j}(1-s_{\mathrm{w},j})a_j.$$

$$v=S_{\mathrm{w}}{\left(\frac{\partial Q}{\partial A_{\mathrm{w}}}\right)}_{A_{\mathrm{n}}}+S_{\mathrm{n}}{\left(\frac{\partial Q}{\partial A_{\mathrm{n}}}\right)}_{A_{\mathrm{w}}}.$$

$${\hat{v}}_{\mathrm{w}}={\left(\frac{\partial Q}{\partial A_{\mathrm{w}}}\right)}_{A_{\mathrm{n}}}$$

$$v=S_{\mathrm{w}}{\hat{v}}_{\mathrm{w}}+S_{\mathrm{n}}{\hat{v}}_{\mathrm{n}},$$

$${\hat{v}}_{\mathrm{w}}=v_{\mathrm{w}}+S_{\mathrm{n}}{v}_m$$

$$\frac{\mathrm{d}v}{\mathrm{d}{S}_{\mathrm{w}}}=v_{\mathrm{w}}-v_{\mathrm{n}}+v_m$$

$$S_{\mathrm{w}}\frac{\mathrm{d}{v}_{\mathrm{w}}}{\mathrm{d}{S}_{\mathrm{w}}}+(1-S_{\mathrm{w}})\frac{\mathrm{d}{v}_{\mathrm{n}}}{\mathrm{d}{S}_{\mathrm{n}}}=v_m.$$

$$v_{\mathrm{w}}=v+S_{\mathrm{n}}\left(\frac{\mathrm{d}v}{\mathrm{d}{S}_{\mathrm{w}}}-v_m\right)$$

$$v_m=a{S}_{\mathrm{w}}+b\frac{\mathrm{d}v}{\mathrm{d}{S}_{\mathrm{w}}}+c.$$

Full text

††journal: Transport in Porous Media

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 d92 ; b72 . A number of factors, such as the capillary forces at the interfaces, viscosity contrast between the fluids, wettability and geometry of the system, make the multiphase flow properties very different compared to single phase flow. When a non-wetting fluid displaces a wetting fluid the flow is called as 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 lmtsm04 ; vhj19 . 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 ltz88 ; mfj85 . The capillary fingering patterns appear during slow displacement process and are well described by invasion percolation ww83 ; lz85 , whereas the viscous fingering patterns appear during fast displacement and can be modeled by diffusion limited aggregation (DLA) model ws81 ; p84 . If both the fluids are continuously fed into the porous medium, the initial transients will eventually die out and the system enters a steady state. It has been discovered that under steady-state conditions, the two-phase flow rate do not obey the linear Darcy law in the capillary dominated regime. Rather, it was found to have a power law dependence on the total pressure drop tkrlmtf09 ; tlkrfm09 ; rcs11 ; sh12 ; sbdk17 .

The problem of two-phase flow in porous media has extensively been studied using laboratory experiments lmtsm04 ; vhj19 ; ltz88 ; mfj85 ; cw85 ; st58 , statistical models ws81 ; ww83 ; bkl86 and numerical simulations vhj19 ; kl85 . There have been significant theoretical developments in describing the steady-state properties hg93 ; h98 ; gm14 ; hsbkgv18 . Due to the vast advancement in the computer power and development in scanning instruments for pore space reconstruction, the computational modeling techniques play a leading role in the study of two-phase flow problems. There are different modeling approaches. The direct numerical simulations, such as the volume-of-fluid method rb12 and the level-set method jh13 ; gm15 , 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 gr91 ; rob10 ; alk11 . 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 when studying systems with large number of pores, such as for the scale-up problems.

Pore-network modeling b01 ; jh12 , in this regard, proves a computationally efficient method that can treat much larger systems. In pore-network modeling, a porous matrix is modeled as 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, here one pore is the smallest computational unit in the pore-network simulations. This simplification looses the detailed evolution of the fluid arrangements inside a pore, however it allows the pore-network models to apply for large systems and thereby study their statistical properties. There are two main ingredients for transporting fluids in a pore-network model: (1) the local pressure drops in the pores and (2) the resulting displacements of the fluids. 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 ww83 ; ltz88 ; ck82 ; b98 where the filling of pores are decided by capillary entry pressures or by determining the stability of an interface for a given contact angle cr88 ; cr90 . Quasi-static models can successfully describe the equilibrium properties of two-phase flow at capillary dominated regime ob98 ; bj02 ; pt18 , 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 jh10 ; hu12 ; ak98 . 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 the local fluid configuration 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 implement the displacements of two fluids in the links of a pore network and to distribute them to the neighboring links after they pass the network nodes. Together with these algorithms, we describe an interface-tracking dynamic pore network model in which the location of two fluids inside the pores are marked explicitly by the positions of all the interfaces at any time step. These interfaces are displaced in small steps under the instantaneous viscous and capillary pressure drops at the pores. We therefore call the algorithms as interface-dynamics algorithms. The displacements of all the interfaces with time describe the transport of the fluids through the network. This is different from other pore-network models, such as in jh10 , where the interface positions are only available implicitly from the volumetric saturation of the fluid elements inside a pore, or from the model in hu12 where an interface is moved through a whole volume element at each time step. Explicit positions of all the interfaces in this interface-tracking model provides a detailed picture of the fluid configuration at any time. Calculation of the capillary forces thus become straightforward from the interface positions and different dynamical events, such as the retraction of invasion fronts after a Haines jump h90 ; bo13 ; ab13 ; mf92 , can easily be resolved.

The interface-tracking model we describe here was first introduced by Aker et al. in late nineties ak98 to study transient two-phase flow – i.e. drainage – in a pore network. The interface algorithms in that model were developed for the drainage in a regular pore network with open boundaries, where a non-wetting fluid invades the network filled with wetting fluid. There, the algorithms were made in such a way that the fluid-volumes were not conserved as they pass across the nodes. Later, the model was extended for steady-state flow by implementing periodic boundary conditions in the flow direction ka02 . There the interface-dynamics algorithms of the previous model were updated to maintain the volume conservation. However, the rules to transport the interfaces through nodes described in the model were based on different events at the nodes of a regular square network in two dimension and were not straightforward to implement in an irregular network topology, such as the networks reconstructed from real samples. Since then, the model has been updated over the years and efforts have been made to generalize the rules for interface dynamics that can be applied for any type of network topology and spatial dimension, as well as for both the transient and steady-state flow sbdk17 . However any detailed description of the new algorithms were absent. Therefore, one of the main goals of this article is to describe these new interface algorithms with sufficient details so that the reader may reproduce the model. The interface algorithms we present here are universal in terms of the network connectivity and topology and can be used for different networks without any further modification. Though these algorithms are simple and straight-forward to implement, they are capable to capture the essential physical properties of both the transient as well as of steady-state immiscible two-phase flow. We will show this by presenting a few fundamental results of transient and steady-state two-phase flow in porous media using this model. 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 interface dynamics rules, while still preserving the fundamental statistical properties of the flow.

We like to point out that, complex flow mechanisms such as the wetting films along the pore walls can also be included in the model as presented in toh12 . However, we do not consider them in this article and restrict ourselves only to the simplest case of Darcy-like creep flow. Moreover, here we will use an explicit Euler integration method to do the time steps, whereas semi-implicit methods can be used for computational efficiency at very low flow rates as illustrated in gvkh18 . A Monte-Carlo algorithm was also developed recently for this model to achieve the steady state in less computational time sshbkv16 . We also note that important phenomena such as wetting angle hysteresis is straight-forward to implement but it has not been included so far. However, wetting angle alteration due to changes in the composition of the immiscible fluids has been considered in the past sgosh11 ; fsh15 .

The article is structured as follows. In section 2, we will present the governing equations describing the two-phase flow at the pore level and will describe how we solve the equations to find the local pressures at the nodes at any time step. Then we present in detail in section 3 the interface-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 described how different boundary conditions can be implemented in this model. In order to validate the model, we will present a few examples in section 5 where the model can successfully reproduce some fundamental results of two-phase flow. Both the transient and steady-state flow will be simulated and the corresponding results will be presented in the subsections 5.1 and 5.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,

[TABLE]

where $Q$ is the total flow rate, $\gamma$ is the interfacial tension between two fluids and $A$ is the cross-sectional pore area of the network. $\mu_{\rm n}$ and $\mu_{\rm w}$ are the viscosities of the two phases. In case of transient studies $\mu_{\rm e}$ is the viscosity of the invading phase, whereas for steady-state flow $\mu_{\rm e}$ is considered as the saturation-weighted effective viscosity of the two phases. Hydraulic properties of a pore-network depend on the geometrical shape of the individual pores, as well as on the network topology, that is, the connectivity and spatial organization of the nodes and links of the network vr01 . With our model, we can consider pore networks in two (2D) or three (3D) dimensions with different topologies as illustrated in figure 1. In (a), we show a crop of a Hele-Shaw cell filled with monolayer of glass beads esth13 that is widely used as a two-dimensional model porous medium in laboratory experiments. Such a porous medium may modeled as a two-dimensional network of disordered pores as shown in (b). A porous medium in 3D, a sample of sand-pack rh09 ; mbweb , is shown in (c) and a network that is reconstructed from the sample is shown in (d). There are different techniques for reconstructing the pore-network from scanned images of the sample, which can be found in a wide range of literature ob05 ; ob02 ; ob03 ; d07 ; db09 .

A pore typically consists two wider pore bodies that are connected by a narrower pore throat. In our network representation, the centers of the pore bodies provide the positions of nodes which are connected by links. In this model, we assign all the volume of the pore space to the links of the network and the nodes do not contain any pore volume. This introduces a variation in the cross-sectional area along the length of the link due to the wider pore bodies at the ends and the narrow pore throat in between. We model the cross section of the link by a simplified hourglass shape as shown in figure 2. The interfacial pressure ($p_{c}$) at an interface therefore depends on the position of the interface as it moves along the link. The functional dependence of $p_{c}$ on the position, obtained from Young-Laplace equation, takes the form shbk13 ,

[TABLE]

where $r_{j}$ is the average radius of a link $j$ and $x_{k}\in[0,l_{j}]$ is the position of the $k$th interface inside the link. Here $\gamma$ is the surface tension and $\theta$ is the contact angle between the interface and the pore wall. If $p_{a}$ and $p_{b}$ are the local pressures at the two nodes $a$ and $b$ across the link, the instantaneous local flow rate $q_{j}$ inside the link from node $a$ to $b$ is proportional to the difference between the viscous pressure drop $\Delta p$ ($=p_{b}-p_{a}$) across the link and the total capillary pressure drop due to all the interfaces inside the link. This can be calculated by w21 ,

[TABLE]

where the sum is over all the interfaces inside the link $j$, taking into account the direction of the capillary forces. $\mu_{j}$ is the saturation-weighted viscosity of the fluids present inside the link at that instant, given by $\mu_{j}=s_{j,{\rm n}}\mu_{\rm n}+s_{j,{\rm w}}\mu_{\rm w}$. Here $s_{j,{\rm n}}=l_{j,{\rm n}}/l_{j}$ and $s_{j,{\rm w}}=l_{j,{\rm w}}/l_{j}$ are the fractions of the link length occupied by the non-wetting and wetting fluids respectively, so that $s_{j,{\rm n}}+s_{j,{\rm w}}=1$. The term $g_{j}$ is the mobility of the link, which depends on the cross section of the link. This model deals with piston-like creep flow at low Reynolds number without any corner flow or film flow. In that case, the effect on $q_{j}$ due to a different cross-sectional shape can be taken into account by the link mobility term $g_{j}$. For the regular network in 2D we chose the links to be cylindrical with circular cross section for which $g_{j}=a_{j}r_{j}^{2}/8$ for Hagen-Poiseuille flow d92 , where $a_{j}=\pi r_{j}^{2}$ is the cross sectional area. In the reconstructed 3D network, the pores are triangular in shape and characterized by a shape factor ($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 $\alpha=\rho^{2}/(4G)$, where $\rho$ is the radius of the inscribed circle in that pore part mm91 . The mobility contribution from each part of the pore is then obtained from the relation $g=3\rho^{2}\alpha/20$ l64 ; jddy08 and the mobility term $g_{j}$ for the total link is then calculated from the harmonic average of the contribution from the three individual terms given by,

[TABLE]

where $\lambda_{1,2,3}$ and $g_{1,2,3}$ are the lengths and mobility of the each pore parts respectively as shown in figure 3.

In order to move the fluids through the pores, the local pressures at nodes are needed to be solved at each time step. The net volume flux ($f_{i}$) through every node at each time step will be zero with the Kirchhoff equations, therefore for any node $i$,

[TABLE]

where the subscript $h$ runs for the links connected to the node $i$. This sum, along with the equations 2 and 3, constructs a set of linear equations. We solve these equations with conjugate gradient solver bh88 or the Cholesky factorization method numrec and the solutions of which at any time step provides the local node pressures $p_{i}$ at that step. The system can be driven by a constant global pressure drop $\Delta P$ or a total flow rate $Q$ with open or periodic boundary conditions. More details about implementing different boundary conditions will be presented in section 4.

3 Interface-dynamics algorithms

We specify the locations and displacements of the two fluids inside the pore space by the positions of all the interfaces. Transport of the fluid bubbles through the links, the coalescence and the snap off of the bubbles – all are treated by the interface-dynamics algorithms by the displacement, creation and deletion of the interfaces. The underlying idea of these algorithms are very simple, at every time step we calculate the volumes of fluids arriving at each node from incoming links and then distribute them to the outgoing links with the volumes proportional to the flow rates in them. In our implementation, the full algorithm for the interface dynamics is sub-divided into three intermediate functions, we name them (a) interface_move, (b) interface_create and (c) interface_merge. Each of these functions is applied at every time step. The purpose of interface_move is to move every interface in each link according to the local link flow rate and to measure the amount of fluids that enter to a node from all the connected links with incoming flow. Next, the interface_create function displaces all the accumulated fluids from each node to the connected links that have the outgoing flow and creates new interfaces at the entrance of those links. Finally the interface_merge function controls the maximum number of interfaces in any link, which can be viewed as an equivalent of the merging and mixing of fluids at the nodes. These three intermediate steps are applied at every links and/or nodes at each time step which together manage the whole process of fluid displacements. We describe them in the following and illustrate them in Fig 4. When describing the algorithms, we will show simulation results at limiting values of parameters to make sure that the rules do not produce any unphysical results. This also helps to select the best possible rule for an algorithm when there can be a few different choices.

3.1 interface_move:

The pore network consists of $N_{\rm L}$ number of links that are connected to each other by $N_{\rm N}$ number of nodes. We will use the subscript $i$ for a node ($i=1,...,N_{\rm N}$) and the subscript $j$ for a link ($j=1,...,N_{\rm L}$). The interface number inside a link is denoted by $k$. Inside the code, one also needs to mark the type of the interfaces, for example, whether it is an interface 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 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 4, we illustrate them for an example of five links connected to a node $i$. At a time step, when all the pressures ($p_{i}$) at the nodes for an existing fluid configuration are known from the solutions of Kirchhoff equations, the flow rates $q_{j}$ inside the links are obtained from equation 3. A time interval ($\Delta t$) is then decided in such a way that an interface inside any link does not move more than $10\%$ of corresponding link length at that time step. In order to do so, the link-velocities $c_{j}$ are calculated from local flow rates $q_{j}$ by $c_{j}=q_{j}/a_{j}$, where $a_{j}$ is the cross sectional area of the link $j$. The time step is then calculated as

[TABLE]

where $c_{j,\rm max}$ is the largest fluid velocity among all the links and $l_{j}$ is the length of that link. For the time evolution, we use an explicit Euler-type procedure. This works well for a large range of capillary numbers with only the time step criterion mentioned in equation 6. In order to ensure numerical stability at low capillary numbers, however, the sensitivity of interfacial pressure jumps to perturbations in interface positions must also be taken into account when choosing the time step, see gvkh18 . This puts a limit on the maximum time step size that is independent of flow rate and thus becomes a severe criterion when flow rates are low. In gvkh18 , computational efficiency was improved at low capillary numbers by a semi-implicit method. Here, however, we consider only capillary numbers that are large enough for equation 6 to be sufficient to ensure numerical stability.

After deciding the time step, interfaces inside each link $j$ are to be moved in the direction of $q_{j}$ by a distance,

[TABLE]

as shown in figure 4(b). Here, the links have uniform cross-sectional area in terms of the mobility, so all the interfaces inside the same link move by the same distance. The positions of all the interfaces inside any link ($x_{j,k}$) are then updated, $x_{j,k}=x_{j,k}+\Delta x_{j}$. By doing this, if any of the interfaces moves outside any link ($x_{j,k}>l_{j}$), it is deleted from the list of interfaces. The set of $p_{i}$s at the nodes decide which links, among all the links connected to a node, have fluids that flow towards the node and have fluids that flow away from the node at that time step. We name these links respectively as the incoming links and the outgoing links for that node. A typical example is shown in figure 4(a) for five links connected to a node $i$. Due to the displacement of the interfaces by an amount $\Delta x_{j}$, a node $i$ receives a certain volume ($\phi_{i,j}$) of wetting and/or non-wetting fluids from an incoming link $j$, given by $\phi_{i,j}=a_{j}\Delta x_{j}=q_{j}\Delta t$ that leaves from the end of the incoming link $j$. The total volume of fluids ($V_{i}$) received by a node $i$ from all of its incoming links is therefore,

[TABLE]

where $V_{i}^{\rm w}$ and $V_{i}^{\rm n}$ respectively are the total volume of wetting and non-wetting fluids arrived from all the incoming links to the node $i$. Here $j\in\rm I$ denotes the set of all the incoming links connected to $i$. We point out that, in order to apply this interface_move function for the whole network, one can loop through all links $j$, move the interfaces by the distance $\Delta x_{j}$ and add up the amount of fluids that exit from the link to the connected node in the direction of the link-flow. This will produce the arrays of $V_{i}^{\rm w}$ and $V_{i}^{\rm n}$ the end of the loop which contain the incoming volume flux at each the nodes at that time step. The marking of the interface types, as mentioned before, will allow here to calculate the individual terms $V_{i}^{\rm w}$ and $V_{i}^{\rm n}$. This is illustrated in figure 4.

3.2 interface_create:

After the function interface_move is performed over all the links, the list of volumes $V_{i}^{w}$ and $V_{i}^{n}$ of incoming fluids for all the nodes are generated for that time step. As all of the pore space in this model is assigned to the links of the network and the nodes do not contain any physical volume, the total volume of incoming fluids ($V_{i}$) to a node are to be injected at the beginning of the outgoing links at the same time step. This will create new bubbles and interfaces at the beginning of the outgoing links. The total volume of fluids that will enter from a node $i$ to an outgoing link $j$ is given by, $\phi_{i,j}=q_{j}\Delta t$, and due to the balancing of the Kirchhoff equations it ensures that $\sum_{j\in\rm I}\phi_{i,j}=\sum_{j\in\rm O}\phi_{i,j}$ where $j\in\rm I$ and $j\in\rm O$ respectively denote the set of all incoming and outgoing links for the node $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 ($\phi_{i,j}^{\rm w}$) and the non-wetting ($\phi_{i,j}^{\rm n}$) volumes will enter into $j$. 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 only on the flow rate $q_{j}$ of the respective outgoing link. For any outgoing link $j$, $\phi_{i,j}^{\rm w}$ and $\phi_{i,j}^{\rm n}$ are therefore calculated as,

[TABLE]

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,

[TABLE]

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 4(c) or in the other way as shown in figure 4(d). 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.

These “democratic” rules are symmetric in terms of the wetting and the non-wetting fluids. Therefore when the surface tension is set to zero, the capillary forces will disappear and one should obtain sg19 ,

[TABLE]

where $Q_{\rm n}$ and $S_{\rm n}$ are the total non-wetting flow rate and the non-wetting saturation. Moreover, if we further set $\mu_{\rm w}=\mu_{\rm n}$ when the surface tension is zero, we should obtain,

[TABLE]

Where, $F_{\rm n}=Q_{\rm n}/Q$, the non-wetting fractional flow in the steady state. These conditions can be used as preliminary tests for the interface functions and to verify their implementations in the code. The precise way to measure the fractional flow and other measurable quantities will be presented in section 5.

3.3 interface_merge:

Creating new interfaces in the outgoing links using the interface_create function can increase the number of bubbles inside a link indefinitely with time. In a real system, coalescence of bubbles occur inside the pores which prevents the number of interfaces from increasing indefinitely. The limit depends on the pore geometry such as the aspect ratio as well as on the flow parameters such as the capillary number and surface tension gf06 . To limit the number of interfaces inside a link in our model, we use an interface_merge function which merges two bubbles when the number of interfaces inside a link exceeds a maximum number ($N_{\rm max}$). After the two functions interface_move and interface_create are executed for each link, we look for any link in which the number of interfaces exceeds the maximum number $N_{\rm max}$. Notice that, limiting the number of interfaces inside a link does not necessarily impose any restriction on the minimum or maximum size of a bubble. The merging of interfaces can be viewed as the coalescence of fluids during they pass the nodes. However, it is a crucial step to formulate how exactly we merge bubbles inside a link so that the volume of the two fluids remain conserved and it does not introduce any artificial effect on the flow properties.

We start with a simple rule 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 towards the flow direction, is moved back towards the other bubble and then merged. The length of this merged non-wetting bubble is then sum of the two non-wetting bubbles before merging. This reduces the interface count by two inside the link and the other interface positions are updated accordingly. This is illustrated in figure 5 (a). However, when we measure the non-wetting fractional flow for zero surface tension with equal viscosities of the fluids, we find a deviation from equation 12. This is shown in figure 6 (a) where a small but systematic deviation from the diagonal $F_{\rm n}=S_{\rm n}$ line can be observed. This effect appears due to displacing the non-wetting bubble opposite to the flow which introduces a decrease in the non-wetting fractional flow. We then tried the next rule B (merge_cm), where instead of merging two nearest non-wetting bubbles, we identified any two nearest interfaces 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 towards the other, here we moved both the bubbles towards 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 5 (b). This rule shows satisfactory $F_{\rm n}=S_{\rm n}$ behavior at the zero capillary pressure as shown in figure 6 (b). However, when we measured $F_{\rm n}$ at finite capillary numbers ${\rm Ca}=0.1$ and $0.01$, we see some discrepancy in the qualitative behavior. The results are shown in the top row of figure 7 with different maximum bubble counts $N_{\rm max}=2$ and $3$. The $F_{\rm n}$ versus $S_{\rm n}$ plot generally shows an S-shaped curve for finite Ca. There the fluids need to overcome the capillary barriers at the narrow pore throats and they flow with different velocities. The phase with larger volume fraction wins and as a result the non-wetting fractional flow curve stays under the diagonal for lower non-wetting saturation and stays above the diagonal for higher non-wetting saturation. Due to the asymmetry between the wetting and non-wetting phases, the curve does not cross the diagonal at the middle. This makes fractional flow curve S-shaped when plotted against the saturation. When the saturation of a fluid is very high, it percolates through the system by trapping the other fluid in small clusters and the fractional flow of the percolating fluid become close to $1$. However, for the merging scheme merge_cm (B), we see that at higher saturation values ($S_{\rm n}>0.8$), the curve again approaches to the diagonal line. We do not see any physical reason for this behavior other than an artificial effect introduced by the merging scheme B. It seems that, as we moved both the bubbles towards each other while merging, small bubbles connected to different links at the nodes got disconnected and started flowing. 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 equation 12 at zero surface tension and also obtained the expected qualitative behavior for non-zero capillary pressure. These are shown in figure 6 (c) and in the bottom row of figure 7 respectively. Moreover, while changing the maximum bubble count with this merging scheme, we find no noticeable difference in the qualitative behavior of $F_{\rm n}$ with the change in $N_{\rm max}$. This is shown in figure 8 for $N_{\rm max}=2$, $3$, $4$ and $5$. We therefore finally adopt the merge_cmnn as the merging scheme.

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 ($p_{i}$) or the link flow rates ($q_{j}$) at the inlets are to be set externally. We also need to set the node fluxes ($V^{\rm w}_{i}$ and $V^{\rm n}_{i}$) at the inlets depending upon the type of fluid injected. For the injection of non-wetting fluid, all the $V^{\rm n}_{i}$s are set to be one and the $V^{\rm w}_{i}$s are set to be zero for all the inlet nodes. With these, we solve the equations for pressures and flow rates for all the other nodes and links inside the network.

The open boundary conditions can also be used for steady-state flow for the setup that is generally used in laboratory experiments tlkrfm09 . 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 $V^{\rm n}_{i}$s are set to one for the non-wetting inlets and all the $V^{\rm w}_{i}$s are set to one for the wetting inlets. The inlet flow rates define the total flow rate ($Q$) and the fractional flow $F_{\rm n}=Q_{\rm n}/(Q_{\rm n}+Q_{\rm w})$. Here, $Q_{\rm n}=\sum^{\prime}_{\rm n}q_{j}$ and $Q_{\rm w}=\sum^{\prime}_{\rm w}q_{j}$ where $\sum^{\prime}_{\rm n}$ and $\sum^{\prime}_{\rm w}$ indicate the sum over all non-wetting and wetting inlets respectively. The fractional flow $F_{\rm n}$ is an external parameter in this setup and the saturation $S_{\rm n}$ is decided by the system.

There is a better way to run steady-state flow in the simulations which is not possible in case of the experiments. This is by implementing the periodic boundary (PB) conditions in the direction of overall flow. For this, we connect the inlet and outlet edges so that the fluids leaving at the outlets, enter the system again through the inlets. The links that connect the inlets and outlets, we call ghost links. With this, the network becomes a closed system and the two fluids keep flowing through the system. The saturation $S_{\rm n}$ therefore is an external parameter here and the fractional flow is decided by the system. This is illustrated in figure 9 for 2D and 3D networks. In case of 2D, periodic boundary conditions are applied in both directions which makes the flow equivalent to the flow on the surface of a torus as shown in (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 interface-dynamics algorithms. Here, the global pressure drop $\Delta P$ to drive the flow needs to be added with the node pressure drops across the ghost links while solving the equations 3 and 5.

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

1:  Network: construct or read

2:  Define: boundary conditions

3:  Initialize: random or sequential fluid distribution

4:  for $t=1$ to $timesteps$ do

5:     Solve the pressure field

6:     for $j=1$ to $totallinks$ do

7:        interface_move($j$)

8:     end for

9:     for $j=1$ to $totallinks$ do

10:        interface_create($j$)

11:        interface_merge($j$)

12:     end for

13:     Calculate measurable quantities at $t$

14:  end for

During initialization, the initial positions of all the interfaces in each link will be defined depending on the saturation and how we want to start the simulation. Then after solving the pressure field, the interface_move function is performed on all the links and that will generate the array of the node fluxes $V^{w}_{i}$ and $V^{n}_{i}$ for each node $i$. These arrays will then be used as input to the interface_create function that will create the new interfaces in the links.

5 Applications and validation

As stated earlier, the interface 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 $45^{\circ}$ with the direction of applied pressure drop in 2D. The length of each link ($l_{j}$) is $1\,{\rm mm}$ and their radii ($r_{j}$) are taken from a uniform distribution of random numbers in the range from $0.1\,{\rm mm}$ to $0.4\,{\rm mm}$. In 3D, we consider a network that is reconstructed from a real sample of Berea sandstone by using micro-CT scanning rh09 ; mbweb . The reconstructed network contains $2274$ links that are connected via $1163$ nodes and has a physical dimension of $1.8^{3}\,{\rm mm}^{3}$. We doubled this network by adding a mirror image of itself in the direction of the applied pressure drop in order to apply periodic boundary conditions. In the following, we present some of the fundamental results of transient and steady-state of two-phase flow in porous media by using this pore-network model. All the following results are obtained only by changing the parameters values and the boundary conditions as necessary, while using the exactly same code for the interface dynamics algorithms.

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 cw85 ; mfj85 and capillary fingering patterns at low Ca lz85 . The capillary and viscous fingering patterns resembles with the invasion percolation ww83 and diffusion limited aggregation (DLA) models ws81 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 ${\rm Ca}$ and viscosity ratio $M$. In our simulations, we set $\gamma\cos\theta=0.03$ N/m. The results are highlighted in figure 10. Here $M=\mu_{\rm n}/\mu_{\rm w}$ where $\mu_{\rm n}$ and $\mu_{\rm w}$ are the viscosities of the invading blue fluid and the wetting gray fluid. In the top row the low viscous non-wetting fluid displaces the more viscous fluid where one can observe the development of the viscous fingering pattern. In the second row $M=10^{2}$, and a more viscous blue fluid is displacing the less viscous gray fluid and a compact and stable displacement front is observed ltz88 . Flow patterns in the third row correspond to capillary fingering which are generated with the simulations at a very low capillary number ${\rm Ca}=10^{-5}$. These fingering patterns are more fractal than the viscous fingering and depend strongly on the system disorder. We like to point out that, these three different transient regimes are generated only by altering the values of the flow rate $Q$ and the viscosities of the two fluids to set the values of Ca and M, and no modification in the interface dynamics algorithms were made between different simulations.

5.2 Steady state

Steady-state flow can be simulated by implementing either open or periodic boundary conditions. In figure 11, we show the evolution of the system to reach the steady state at a constant flow rate with different boundary conditions. The top row shows simulations with open boundary condition, where wetting (gray) and non-wetting (blue) fluids are injected through alternate injection points at the bottom edge of a two-dimensional pore network consisting $64\times 100$ links. The radii ($r_{j}$) of the links are drawn from a uniform random number distribution in the range $0.1l$ to $0.4l$, where $l_{j}$ is the length of the links. The top edge of the network is kept open, through which fluids leave the system. The two side edges are connected with periodic boundary here but can also be kept closed if necessary. We can control the flow rates of the individual fluids which sets the fractional flow. Here the simulations are shown for $F_{\rm w}=0.5$. As there will be traces of injection of the alternate fluids near the inlet edge, one consider a system that is long enough in the direction of the overall flow so that a region of spatially homogeneous steady-state flow can be achieved away from the inlets. At the right, we plot the average global pressure drops as a function of the injected pore-volumes which show the evolution to the steady state when $\Delta P$ fluctuate around an average value. The second row of figure 11 shows simulations with periodic boundary conditions in a disordered square lattice of $64\times 64$ links. Here the system is closed and the control parameter is the fluid saturation. The simulations are shown for $S_{\rm w}=0.5$. In the bottom row, we show the simulation with a three dimensional pore network reconstructed from a real sample of Berea sandstone for $S_{\rm w}=0.5$.

In the following, we present two examples of simulation results with this dynamic pore network model to show excellent 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 tlkrfm09 ; rcs11 ; sh12 ; sbdk17 . Next, we measure seepage velocities and will verify the relations between them hsbkgv18 . We will also describe in detail how to measure different quantities of flow, such as the flow rates and the seepage velocities of the different fluid components.

5.2.1 Effective rheology and crossover from linear to

non-linear flow regime

Experiments have shown that the two-phase flow in the steady state do not follow the linear Darcy relation between the total flow rate and the pressure drop in the steady state in capillary dominated regime tkrlmtf09 ; tlkrfm09 ; rcs11 . The capillary pressures at the interfaces in two phase flow create threshold barriers. The distribution of these thresholds over the system depends on the pore sizes as well and on the interface configurations in the steady state, which in turn depend on the flow dynamics. These threshold pressures create a global yield threshold pressure ($P_{t}$) below which there is no flow through the system. Above $P_{t}$, the pores start opening while increasing the pressure drop, due to which the total flow rate no longer vary linearly with the pressure drop $\Delta P$ and shows non-linear behavior. This non-linear behavior depends on the distribution of the threshold pressure as well as on the system geometry rhs19 . For two Newtonian fluids flowing through a porous media, it was found experimentally tkrlmtf09 ; tlkrfm09 ; rcs11 ; sbdk17 , numerically sh12 ; sbdk17 and with mean field calculations sh12 that, in the regime when capillary pressures compete with the viscous forces, that is, in the low capillary number regime, the total flow rate $Q$ in the steady state vary quadratically with the excess pressure drop. We can write as the following,

[TABLE]

The quadratic regime corresponds to when the individual pores keep opening with the increasing pressure drop. Due to this, more and more flow paths appear in the system which makes the flow rate to increase faster than $\Delta P$. When the capillary number is high enough, there is a transition from the quadratic regime and the flow rate becomes linear with the pressure drop. To verify whether our model can produce this cross-over, we performed steady-state simulations at constant flow rate $Q$. The results are presented in figure 12. The threshold pressures $P_{t}$ are calculated by plotting $\Delta P$ versus $\sqrt{Q}$ as shown in the insets. When the results follow equation 13, it will produce straight lines for the lower values of $\Delta P$, where the intercepts of the straight lines at $y$-axis correspond to the values of $P_{t}$. Using these values of $P_{t}$, we plot $Q$ versus $(\Delta P-P_{t})$ in the log-log scale the disordered 2D network and the reconstructed 3D network. Results are shown in figure 12 which show two distinct regimes, a quadratic regime with slope $\approx 2$ at the low pressure drops and a linear regime with slope $\approx 1$ at higher pressure drops. A detailed numerical study on this non-linear effective rheological properties using the interface algorithms presented here can be found in sbdk17 , where the results are also compared with experiments.

5.2.2 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 $\Delta P$, a set of equations relating the wetting and non-wetting seepage velocities ($v_{\rm w}$, $v_{\rm n}$) to the total seepage velocity $v$ and the fluid saturations can be derived using the Euler homogeneity property of the total flow rate $Q$ in the steady state hsbkgv18 . These relations necessitate a new velocity function, namely the co-moving velocity ($v_{m}$), which is a characteristic of the porous medium. The seepage velocities for the wetting and non-wetting fluids are defined as

[TABLE]

respectively, where $Q_{\rm w}$ and $Q_{\rm n}$ are the volumetric flow rates of the two fluids in the direction of the applied pressure drop. Quantitatively, $Q_{\rm w}$ and $Q_{\rm n}$ are defined as the volume of the wetting and non-wetting fluids that pass through any cross section of the system, perpendicular to the overall flow direction, per unit time. $A_{\rm w}$ and $A_{\rm n}$ are the wetting and non-wetting pore areas defined as the areas occupied by the wetting and non-wetting fluids along any orthogonal cross section through the system. This is illustrated in figure 13 where $\Delta P$ is applied in the positive $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 $A_{\rm w}$ and $A_{\rm n}$ along this cross section. For a homogeneous porous medium the average values of $A_{\rm w}$ and $A_{\rm n}$ remain same for any orthogonal cross section of the system in the steady state. Here we measure $A_{\rm w}$ and $A_{\rm n}$ by averaging the pore areas over all possible cross sections along $x$. The total pore area $A$ related to all the fluids is therefore given by $A=A_{\rm w}+A_{\rm n}=\phi A_{s}$ where $\phi$ is the porosity and $A_{s}$ is the average cross-sectional area of the total system including the pore space and the solid. With this, we can express the fluid saturations $S_{\rm w}$ and $S_{\rm n}$ in terms of the pore areas by $S_{\rm w,n}=V_{\rm w,n}/V_{p}=(A_{\rm w,n}L)/(AL)=A_{\rm w,n}/A$. The total flow rate $Q$ of the two fluids is the sum of the wetting and non-wetting flow rates given by $Q=Q_{\rm w}+Q_{\rm n}=A_{\rm w}v_{\rm w}+A_{\rm n}v_{\rm n}$. Correspondingly, the total seepage velocity $v$ associated with the total flow rate $Q$ is defined by,

[TABLE]

and we can find,

[TABLE]

by using the relations mentioned above.

In our network simulations we have the information about the local flow rates $q_{i}$ and the interface positions for each link $i$ at any time step. However, calculating the average flow rates and the pore areas along any orthogonal cross section of the network from those quantities may not be straight forward, specially in case of an irregular network. For a regular network considered in 2D, all the links are of the same length $l_{j}=1\,{\rm mm}$ and they are oriented along the same angle ($45^{\circ}$) with respect to the overall flow direction. The sum of the local flow rates through each row of links normal to the flow direction is therefore the same for any row and the flow rates can therefore be measured by summing over all the links of the network and then dividing by the number of rows. This given by,

[TABLE]

where $N_{L}$ is the number of rows along $L$, i. e. $64$ here. Similarly, we can calculate the cross-sectional areas as,

[TABLE]

where $a_{j}$ is the cross-sectional area of the link $j$, projected into the plane normal to the flow direction. The wetting saturations of the links $s_{{\rm w},j}$ are provided by the interface positions. Here the individual terms corresponding to the wetting and non-wetting phases are multiplied with the corresponding link saturations as the probability that a cross section through a link will pass through the wetting or non-wetting phase is proportional to the link saturation of that phase.

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 the equations 18 leads 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 (figure 13). The probability that any link $j$ will pass through this plane will be proportional to $l_{x,j}/L$ where $l_{x,j}$ is the length of the link $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 $s_{{\rm w},j}$ of the link. Considering these probabilistic terms, the total flow rates $Q$, $Q_{\rm w}$ and $Q_{\rm n}$ through a random orthogonal cross section can therefore be calculated from the sum of the local flow rates over the links which pass through this cross section,

[TABLE]

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

[TABLE]

For the regular 2D network, $l_{x,j}$ are the same for all the links ($=l$) and we recover equations 17 and 18 by using $N_{L}=L/l$. After measuring the flow rates and the pore areas, the seepage velocities $v$, $v_{\rm w}$ and $v_{\rm n}$ are calculated using the equations 15 and 27. Results are averaged over time in the steady state.

The calculation of $v$, $v_{\rm w}$ and $v_{\rm n}$ from the measurements of the flow rates and the pore areas should satisfy equation 16. In figure 14, we plot $S_{\rm w}v_{\rm w}+S_{\rm n}v_{\rm n}$ against the total seepage velocity $v$ where we used the equations 17 and 18 for 2D and the equations 19 and 20 for 3D. Figure 14 shows an exact match of equation 16 for the whole range of parameters.

The total flow rate $Q$ in the steady state is a homogeneous function of order one of the pore areas $A_{\rm w}$ an $A_{\rm n}$, that is, if we scale the three areas by $A\to\lambda A_{\rm w}$, $A_{\rm n}\to\lambda A_{\rm n}$ and $A_{s}\to\lambda A_{s}$ by keeping the porosity $\phi$ constant, the volumetric flow rate $Q$ scales as, $Q(\lambda A_{\rm w},\lambda A_{\rm n})=\lambda Q(A_{\rm w},A_{\rm n})$. This property of $Q$ leads to a new set of equations between the seepage velocities. Complete derivations of the equations can be found in reference hsbkgv18 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 $\lambda$ and then setting $\lambda=1$ we get,

[TABLE]

These two partial derivatives in the above equation have the units of velocity and correspondingly they define two thermodynamic velocities $\hat{v}_{\rm w}$ and $\hat{v}_{\rm n}$ given by,

[TABLE]

With these definitions equation 21 becomes,

[TABLE]

which has the similar form of equation 16. However, this does not imply that the thermodynamic velocities $\hat{v}_{\rm w}$ and $\hat{v}_{\rm n}$ are the same as the seepage velocities $v_{\rm w}$ and $v_{\rm n}$ that we measure. These two types of velocities can be related by a new velocity function is $v_{m}$ given by,

[TABLE]

which fulfill both the equations 16 and 23. This velocity function $v_{m}$ is a function of the saturation $S_{\rm w}$ and called as the co-moving velocity, which is a property of the pore-network. With this definition of $v_{m}$, we can derive two equations that are related to the variation of saturation,

[TABLE]

and

[TABLE]

In order to verify whether our pore-network model with the set of interface 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=0.5$, $1$, $2$, $5$ and $10$ are considered where the wetting viscosity is chosen as $\mu_{\rm w}=0.1\,{\rm Pa\,s}$. The non-wetting viscosity $\mu_{\rm n}$ is then chosen accordingly in order to set the value of $M$. For each value of $M$, three different values of the surface tension, $\gamma=0.02$, $0.03$ and $0.04\;{\rm N/m}$, are considered. For each set of $M$ and $\gamma$, we considered a set of pressure drops $|\Delta P/L|=0.16$, $0.20$, $0.40$, $0.80$, $1.0$ and $2.0\;{\rm MPa/m}$ for 2D and $10$, $20$, $40$, $80$ and $160\;{\rm MPa/m}$ for 3D. These values of pressure drops are chosen in order to get capillary numbers ${\rm Ca}$ in a range around $10^{-3}$ to $10^{-1}$. Specifically, we find ${\rm Ca}$ in the range of $0.0007-0.2985$ for 2D and $0.001-0.271$ for 3D. For any set of parameters, saturations are varied in the range of [math] to $1$ in the steps of $0.05$ which correspond to $21$ saturation values. This led to a total of $1890$ independent simulations for 2D and $1575$ simulations for 3D in steady state.

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 $v_{m}=(\mathrm{d}v/\mathrm{d}S_{\rm w})-v_{\rm w}+v_{\rm n}$ (equation 25) and $v_{m}=S_{\rm w}(\mathrm{d}v_{\rm w}/\mathrm{d}S_{\rm w})+S_{\rm n}(\mathrm{d}v_{\rm n}/\mathrm{d}S_{\rm n})$ (equation 26) and plot in figure 15. The results show good agreement with the equation 25 and 26 for both square and Berea networks for the whole range of the capillary numbers as indicated by the color scale. The few data points that are outside the straight line mostly correspond to the simulations near $S_{\rm w}=0$ or $1$ where the system undergoes from two-phase to single phase regime that creates a jump in the derivatives.

Equations 16 and 26 transforms the wetting and non-wetting velocities ($v_{\rm w}$, $v_{\rm n}$) to ($v$, $v_{m}$) while varying the saturation. By inverting the velocity transformation, it is possible to find equations to transform the total velocity and the co-moving velocity ($v$, $v_{m}$) to ($v_{\rm w}$, $v_{\rm n}$), which are given by,

[TABLE]

With these equations we can verify the measured values of wetting and non-wetting seepage velocities $v_{\rm w}$ and $v_{\rm n}$ against the ones calculated from these equation. In figure 16, we plot $v+S_{\rm n}(\mathrm{d}v/\mathrm{d}S_{\rm w}-v_{m})$ and $v-S_{\rm w}(\mathrm{d}v/\mathrm{d}S_{\rm w}-v_{m})$ against the measured value of $v_{\rm w}$ and $v_{\rm n}$ respectively. While calculating $v_{\rm w}$ and $v_{m}$ using above equations, we used the values of $v_{m}$ that are obtained from equation 26. For the whole range of the capillary numbers, good match with equations 27 can be observed for both the square and the Berea networks.

Finally, we plot the co-moving velocity $v_{m}$ as a function of both $S_{\rm w}$ and $\mathrm{d}v/\mathrm{d}S_{\rm w}$ in figure 17. The co-moving velocity is a property of the porous material and a function of the saturation $S_{\rm w}$, total seepage velocity $v$ and the variation of $v$ while changing the saturation. It is not enough to specify only the $S_{\rm w}$ and $v$ to determine $v_{m}$, as $\mathrm{d}v/\mathrm{d}S_{\rm w}$ depends on how the external parameters are controlled while varying the $S_{\rm w}$. Here $v_{m}$ was calculated using equation 25 in figure 17. The data for $v_{m}$ roughly shows a planer form given by,

[TABLE]

By fitting all the data points for the whole set of simulations, we find $a=-6.36\pm 0.25$, $b=0.94\pm 0.01$ and $c=5.00\pm 0.13$ for the square network and $a=-12.94\pm 0.62$, $b=0.88\pm 0.01$ and $c=10.10\pm 0.32$ for the Berea network. The planes using these parameters are shown in the respective figures with grid lines. Interestingly, the values of $b$ are close to $1$ which leads the data points to fall around the $y=x$ straight line while plotting $v_{m}$ against $\mathrm{d}v/\mathrm{d}S_{\rm w}$ as shown the figure 17.

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 interfaces with time. The basic concept of the algorithms are simple, at every time step all the fluids arriving at a node from the incoming links are displaced to all the outgoing links, and the volumes of the fluids are distributed according to the ratio of the fluxes of the outgoing links. Our aim in this article is to present 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 interface algorithms are applicable for both the regular and irregular network topologies as well as for both two and three dimensional pore networks. Moreover, by reproducing some of the fundamental results of two-phase flow, we have also shown that the model can be used to simulate both the transient and steady-state flow. We have shown different drainage flow patterns that can be generated with this model when a fluid displaces the other in a porous medium. 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.

In a recent paper by Zhao et al. vhj19 ten different groups with different approaches to modeling 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. Film flow is an important mechanism during imbibition and in a first attempt, we expanded our model in toh12 to include films. However, this work has so far not been followed up by us.

The time integration procedure is time consuming. In sshbkv16 , Savani et al. proposed a Monte Carlo algorithm to replace the time integration. This approach promises a large increase in efficiency of the model. However, the method needs to be tested at low capillary numbers. Moreover, it has so far only been implemented for regular lattices and needs to be generalized to realistic pore networks.

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, M. 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å. This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. SS was supported by the National Natural Science Foundation of China under grant number 11750110430.