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

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.


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][12][13][14][15][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 steadystate properties [22][23][24][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][30][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 voxelbased 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 porenetwork 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 twophase flow at capillary dominated regime [38][39][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][42][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 porenetwork 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][45][46][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 flowi.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 twophase 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 and 5.2 respectively. Finally, we will draw our conclusions in Section 6.

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 twophase 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, Ca μ e Q cA and M μ n μ w (1) 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. μ n and μ w are the viscosities of the two phases. In case of transient studies μ e is the viscosity of the more viscous phase, whereas for steady-state flow μ 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 [52]. 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 a monolayer of glass beads [50] that is widely used as a twodimensional 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][54][55][56][57].
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 crosssectional area along the length of the link which is modeled by a simplified hourglass shape as shown in Figure 2B. The interfacial pressure (p c ) therefore depends on the position of the meniscus as it moves along the link. The functional dependence of p c on the position, obtained from Young-Laplace equation, takes the form [58], where r j is the average radius of a link j and x k ∈ [0, l j ] is the position of the kth meniscus inside the link. Here γ is the surface tension and θ is the contact angle between the meniscus 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 Δp ( p b − p a ) across the link and the total capillary pressure drop due to all the menisci inside the link. This can be calculated by [59], where the sum is over all the menisci inside the link j, taking into account the direction of the capillary forces. μ j is the saturationweighted viscosity of the fluids present inside the link at that instant, given by μ j s j,n μ n + s j,w μ w . Here s j,n l j,n /l j and s j,w l j,w /l j are the fractions of the link length occupied by the non-wetting and wetting fluids respectively, so that s j,n + s j,w 1. The term g j is the mobility of the link that depends on the cross section. We model piston-like creep flow at low Reynolds number without any film flow and different cross-sectional shapes can therefore be taken into account by this 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 2 j /8 for Hagen-Poiseuille flow [1], where a j πr 2 j is the cross-sectional area. In the reconstructed 3D network, the pores are triangular in shape that is 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 α ρ 2 /(4G), where ρ 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 g 3ρ 2 α/20 [61,62] 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, where λ 1,2,3 and g 1,2,3 are the lengths and mobilities of the pore parts respectively as shown in Figure 2A.
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 (f i ) through every node at each time step will be zero, that is, any node 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 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. λ 1,2,3 and g 1,2,3 are the lengths and mobility contributions from these three parts. The total length of the link is given by l j λ 1 + λ 2 + λ 3 and the total mobility of the link is obtained from Eq. 4. This total pore space is modeled by an hour-glass shaped link as shown in (B) where the wetting and non-wetting fluids are colored by white and gray respectively. The variation in the interfacial pressure (p c ) due to this shape is modeled by Eq. 2 as shown above. Here p c is modeled by a cosine function, however one may consider some other functional form with a maximum in the middle that corresponds to the minimum of the pore radii at the narrow pore throat.
Frontiers in Physics | www.frontiersin.org March 2021 | Volume 8 | Article 548497 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 Δ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.

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 (N max ). This function is necessary as when fluids enter in the links from the nodes, new menisci are created at the entrance of the links which can make the meniscus number to increase indefinitely during the time stepping. In a real system, coalescence of bubbles occur inside the pores which limits the number of menisci. There are experimental studies to investigate the formation of bubbles and to find their characteristic size during flow through pore junctions [67], however, to our knowledge, there is no clear information available on the functional dependency of the maximum number of menisci or bubbles on the pore geometry and aspect ratio as well as on other flow parameters such as the fluid velocity, contact angle and surface tension. In order to get a detailed arrangement of fluids inside a pore, one needs to solve the flow equations at a much smaller scale, such as in Lattice Boltzmann simulations. In pore-network modeling, we are interested in the large-scale flow properties and not in the detailed configurations of fluids inside a pore. We therefore kept this number (N max ) as a tuning parameter in our algorithm which may be calibrated based on experimental observations. We constructed a merging algorithm that keeps the center of mass of the bubbles unchanged so that the overall fluid displacement does not get affected by the merging. We have performed a set of sensitivity test on the maximum meniscus number which we will present in 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.

meniscus_move
The pore network consists of N L number of links that are connected to each other by N N number of nodes. We will use the subscript i for a node (i 1, . . . , N N ) and the subscript j for a link (j 1, . . . , N L ). The meniscus number inside a link is denoted by 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 (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 Eq. 3. A time interval (Δt) is then decided in such a way that a meniscus 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 where c j,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 Eq. 6. In order to ensure numerical stability at low capillary numbers, however, the sensitivity of interfacial pressure jumps to perturbations in meniscus positions must also be taken into account when choosing the time step, see [68]. 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 [68] computational efficiency was improved at low capillary numbers by a semi-implicit method for a different type of meniscus algorithms. Here, however, we consider only capillary numbers that are large enough for Eq. 6 to be sufficient to ensure numerical stability. We will show some simulation results in Section 5.1 for Ca 10 − 5 and in our tests, we could run simulations at capillary numbers as low as Ca 10 − 6 , though the simulations were highly time consuming. Simulations with capillary numbers lower than that are not doable within reasonable computational time. . This is performed by the function meniscus_move. The total volume of the fluids (a 1 Δx 1 + a 2 Δx 2 + a 3 Δx 3 ) which left the links 1, 2 and 3 into the node i is now to be placed at the beginning of the links 4 and 5, which will fill the volume (a 4 Δx 4 + a 5 Δx 5 ) and will create new menisci in the outgoing links as described in meniscus_create. There are two different ways of creating the new menisci, one is to place the wetting fluid first and the non-wetting next as shown in (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. In (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. After deciding the time step, menisci inside each link j are to be moved in the direction of q j by a distance, as shown in Figures 3A,B. Here, the links have uniform crosssectional 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 (x j,k ) are then updated, x j,k x j,k + Δx j . By doing this, if any of the menisci moves outside any link (x j,k > l j ), it is deleted from the list of menisci. The set of p i s at the nodes decide which links, among all the links connected to a node, have fluids that flow toward 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. 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 Δx j , a node i receives a certain volume (ϕ i,j ) of wetting and/or non-wetting fluids from an incoming link j, given by ϕ i,j a j Δx j q j Δt that leaves from the end of an incoming link j. The total volume of fluids (V i ) received by a node i from all of its incoming links is therefore, where V w i and V n i respectively are the total volume of wetting and non-wetting fluids arrived from all the incoming links to the node i. Here j ∈ I denotes the set of all the incoming links connected to 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 Δ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 w i and V n i at the end of the loop which will contain the incoming volume flux at each of the nodes at that time step. The marking of the meniscus types, as mentioned before, will allow here to calculate the individual terms V w i and V n i . This is illustrated in Figure 3.

meniscus_create
After the function meniscus_move is performed over all the links, the list of the volumes V w i and V n i 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 menisci 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, ϕ i,j q j Δt, and due to the balancing of the Kirchhoff equations it ensures that j ∈ I ϕ i,j j ∈ O ϕ i,j where j ∈ I and j ∈ 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 (ϕ w i,j ) and the non-wetting (ϕ n i,j ) volumes will enter into 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 q j of the respective outgoing links. For any outgoing link j, ϕ w i,j and ϕ n i,j are therefore calculated as, 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 Q n and S n are the total non-wetting flow rate and the nonwetting saturation. If we further set μ w μ n when the surface tension is zero, the two fluids essentially become identical which will flow with equal velocity. This leads to, where, F n Q n /Q, the non-wetting fractional flow in the steady state. These conditions can be used as preliminary tests for the meniscus functions whether they are implemented properly in the code. The precise way to measure the fractional flow and other measurable quantities will be presented in Section 5.

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 N max . Notice that, restricting the number of menisci inside a link does not necessarily impose any restriction on the minimum or maximum size of a bubble. In order to find out the best possible algorithm and to test its sensitivity, we performed numerical simulations with a network of 64 × 64 links forming a diamond lattice in 2D with disorder in the link radii. We start with a simple merging 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 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 nonwetting 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 F n S n line can be observed which is unphysical. This effect is a consequence of 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 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 F n S n behavior at the zero capillary pressure as shown in Figure 5B. We then measured F n at finite capillary numbers Ca 0.1 and 0.01, and check the sensitivity and the qualitative behavior of the results. They are shown in the top row of Figure 6 for maximum bubble counts N max 2 and 3. The F n vs. S n plot generally shows an S-shaped behavior for finite Ca in experiments and simulations [70,71]. There are capillary barriers associated with the narrow pore throats which the fluids need to overcome. Consequently, the fluids flow with different velocities and the phase with larger volume fraction wins. The non-wetting fractional flow curve therefore stays under the diagonal for lower non-wetting saturation and above the diagonal for higher nonwetting saturation. Due to the direction of capillary forces between the wetting and non-wetting phases, the curves do 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 the curve deviates from a typical S-shape and approaches to the diagonal line at higher saturation values (S n > 0.8). This seems an artificial effect introduced by the merging scheme 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 F n with the change in N max . This is shown in Figure 7 for N max 2, 3, 4 and 5. We therefore finally adopt the merge_cmnn as the merging scheme.   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 w i and V n i ) at the inlets depending upon the type of fluid injected. For the injection of non-wetting fluid, all the V n i s are set to be one and the V 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.

BOUNDARY CONDITIONS
The open boundary conditions can also be used for steadystate 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 V n i s are set to one for the non-wetting inlets and all the V 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 n Q n /(Q n + Q w ). Here, Q n ′ n q j and Q w ′ w q j where Σ ′ n and Σ ′ w indicate the sum over all non-wetting and wetting inlets respectively. The fractional flow F n is an external parameter in this setup and the saturation S n is decided by the system.
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 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.
Frontiers in Physics | www.frontiersin.org March 2021 | Volume 8 | Article 548497 system and the two fluids keep flowing through the system. The saturation S n therefore is a control parameter here and the fractional flow is decided by the system. This is illustrated in Figure 8 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 meniscus-dynamics algorithms. Here, the global pressure drop ΔP to drive the flow needs to be added with the node pressure drops across the ghost links while solving Eqs. 3, 5.
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 t 1 to timesteps do (5) Solve the pressure field (6) for j 1 to totallinks do (7) meniscus_move(j) (8) end for (9) for j 1 to totallinks 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 V w i and V n i for each node 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.

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 45 + with the direction of applied pressure drop in 2D. The length of each link (l j ) is 1 mm and their radii (r j ) are taken from a uniform distribution of random numbers in the range from 0.1 mm to 0.4 mm. In 3D, we consider a network that is reconstructed from a real sample of Berea sandstone by using micro-CT scanning [49,51]. The reconstructed network contains 2,274 links that are connected via 1,163 nodes and has a physical dimension of 1.8 3 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 parameter values and the boundary conditions as necessary, while using the exactly same algorithms for meniscus displacements.

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 Ca and viscosity ratio M. In our simulations, we set c cos θ 0.03 N/m. The results are highlighted in Figure 9. Here M μ n /μ w where μ n and μ w are the viscosities of the invading blue non-wetting fluid and the defending gray wetting 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 [5]. Flow patterns in the third row correspond to capillary fingering which are generated with the simulations at a very low capillary number Ca 10 − 5 . These fingering patterns are more fractal than the viscous fingering and depend strongly on the system disorder. For a more quantitative measurement of the different regimes, one can measure the fractal dimensions of the fingering patterns. 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 meniscus dynamics algorithms were made between different simulations.

Steady State
Steady-state flow can be simulated by implementing either open or periodic boundary conditions as shown in Figure 10.  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 (ΔP). By plotting ΔP as a function of the injected pore-volumes as shown in the right of each row, we identify the steady states, that is when ΔP fluctuates around an average value.
In the following, we present two sets of simulation results with this dynamic pore network model to show the agreement with the steadystate 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][13][14][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 Nonlinear Flow Regime
When two Newtonian fluids flow through a porous media, it was found experimentally [11][12][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 (P t ) is a global yield threshold pressure (P t ) below which there is no flow through the system. Study also shows that the non-linear behavior depends on the distribution of the threshold pressures as well as on the system geometry [72]. This is related to the capillary pressures at the menisci that create threshold barriers at the pore throats. These pores start opening while increasing the pressure drop, which leads the total flow rate Q to increase faster than the increase in ΔP. At sufficiently high pressure drop, all the pores start flowing and there is a crossover to a linear regime. 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 ΔP is measured in the steady state. Here the steady-state was achieved by using periodic boundary conditions. For 2D, the data is averaged over 20 different samples. The results are presented in Figure 11 where we plot the overall flow rate Q as a function of (ΔP − P t ). The threshold pressures P t are calculated by plotting ΔP vs. Q √ as shown in the insets. When the results follow Eq. 13, it will produce straight lines for the lower values of ΔP, where the intercepts of the straight lines at y-axis correspond to the values of P t . The results show two distinct regimes, a quadratic regime with slope ≈ 2 satisfying Eq. 13 at the low pressure drops and a linear regime with slope ≈ 1 at high pressure drops. More simulation results on this non-linear rheological properties for different flow parameters and network geometries can be found in [15], where the meniscus algorithms were used for steady-state flow and the results were compared with experiments.

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 ΔP, a set of equations relating the wetting and non-wetting seepage velocities (v w , v 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 [25]. These relations necessitate a new velocity function, namely the comoving velocity (v m ), which is a characteristic of the porous medium. The seepage velocities for the wetting and non-wetting fluids are defined as v w Q w A w and v n Q n A n (14) respectively, where Q w and Q n are the volumetric flow rates of the two fluids in the direction of the applied pressure drop. Quantitatively, Q w and Q 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 w and A n are the wetting and nonwetting 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 12 where Δ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 nonwetting pore areas A w and A n along this cross section. For a homogeneous porous medium the average values of A w and A n remain same for any orthogonal cross section of the system in the steady state. Here we measure A w and A 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 FIGURE 11 | Variation of the total volumetric flow rate Q (mm 3 /s) with the overall pressure drop ΔP (kPa) in the steady state for square and Berea networks. In the inset, ΔP is plotted against Q √ for the low Ca regime where the intercepts at the y-axis correspond to the values of threshold pressures (P t ). For 2D, we find P t 3.65kPa and 4.09 kPa for S n 0.3 and 0.4 respectively. For the Berea network in 3D, we find P t 3.54kPa and 0.32kPa for S n 0.3 and 0.4 respectively. Using these values, we plot log 10 Q vs. log 10 (ΔP − P t ) which shows two distinct regimes at low and high pressures. For the linear regime, the slopes are obtained as 0.99 ± 0.01 and 1.00 ± 0.01 for 2D and 1.03 ± 0.01 and 1.04 ± 0.01 for 3D for the two saturations respectively. For the quadratic regime, the slopes are obtained as 1.96 ± 0.02 and 1.98 ± 0.03 for 2D and 1.98 ± 0.03 and 1.99 ± 0.04 for 3D.
A w + A n ϕA s where ϕ 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 w and S n in terms of the pore areas by S w,n V w,n /V p (A w,n L)/(AL) A w,n /A. The total flow rate Q of the two fluids is the sum of the wetting and nonwetting flow rates given by Q Q w + Q n A w v w + A n v n . Correspondingly, the total seepage velocity v associated with the total flow rate Q is defined by, and we can find, by using the relations mentioned above. In our network simulations, we have the information about the local flow rates q i and the meniscus positions for each link i at any time step. When the system is driven under constant pressure drop (ΔP), calculating the average flow rates and the pore areas along any orthogonal cross section of the network from the local flow rates are not straight forward, specially in case of an irregular network. For a regular network in 2D, all the links are of the same length l j 1 mm and they are oriented along the same angle (45 + ) 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,

Q
1 N L j q j , Q w 1 N L j s w,j q j and Q n 1 N L j 1 − s w,j q j (17) where N L is the number of rows along L, i. e. 64 here. Similarly, we can calculate the cross-sectional areas as, 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 w,j are provided by the meniscus positions. Here the individual terms corresponding to the wetting and nonwetting 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 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 FIGURE 12 | Description of the system to measure the flow rates (Q, Q w , Q n ), the pore areas (A, A w , A n ) and the seepage velocities (v, v w , v n ) for the 2D and 3D networks. The global pressure drop ΔP is applied in the 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, A w and A n . In this figure, the gray and blue patches shown in the normal view of the cross sections do not reflect the actual occupations of the fluids in the above networks and are given as illustration purpose only.
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 w,j of the link. Considering these probabilistic terms, the total flow rates Q, Q w and Q 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, 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, l x,j are the same for all the links ( l) and we recover Eqs 17, 18 by using N L L/l. After measuring the flow rates and the pore areas, the seepage velocities v, v w and v n are calculated using Eqs. 14, 15. Results are averaged over time in the steady state. The calculation of v, v w and v n from the measurements of the flow rates and the pore areas should satisfy Eq. 16. In Figure 13 (top row), we plot S w v w + S n v n against the total seepage velocity 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.
The total flow rate Q in the steady state is a homogeneous function of order one of the pore areas A w an A n , that is, if we scale the three areas by A → λA w , A n → λA n and A s → λA s by keeping the porosity ϕ constant, the volumetric flow rate Q scales as, Q(λA w , λA n ) λQ(A w , A 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 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 λ 1 we get, v S w zQ zA w An + S n zQ zA n Aw . These two partial derivatives in the above equation have the units of velocity and correspondingly they define two thermodynamic velocities v w and v n given by, v w zQ zA w An and v n zQ zA n Aw .
With these definitions Eq. 21 becomes, which has the similar form of Eq. 16. However, this does not imply that the thermodynamic velocities v w and v n are the same as the seepage velocities v w and v n that we measure. These two types of velocities can be related by a new velocity function is v m given by, which fulfill both Eqs. 16, 23. This velocity function v m is a function of the saturation S 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, dv 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.  25) and v m S w (dv w /dS w ) + S n (dv n /dS n ) (Eq. 26) and plot in Figure 13 (bottom row). The results show good agreement with Eqs. 25, 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 w 0 or 1 where the system undergoes from two-phase to single phase regime that creates a jump in the derivatives.
Eqs. 16, 26 transform the wetting and non-wetting velocities (v w , v 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 w , v n ), which are given by, With these equations we can verify the measured values of wetting and non-wetting seepage velocities v w and v n against the ones calculated from these equation. In Figure 14, we plot v + S n (dv/dS w − v m ) and v − S w (dv/dS w − v m ) against the measured value of v w and v n respectively. While calculating v w and v m using above equations, we used the values of v m that are obtained from Eq. 26. For the whole range of the capillary numbers, good match with Eq. 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 w and dv/dS w in Figure 15.
Here v m was calculated using Eq. 25 in Figure 15. The data for v m roughly shows a planer form given by, v m aS w + b dv dS w + c.
By fitting all the data points for the whole set of simulations, we find a −6.36 ± 0.25, b 0.94 ± 0.01 and c 5.00 ± 0.13 for the square network and a −12.94 ± 0.62, b 0.88 ± 0.01 and c 10.10 ± 0.32 for the Berea network. The planes using these parameters are shown in the respective figures with grid lines. The co-moving velocity is a property of the porous material and a function of the saturation S w , total seepage velocity v and the variation of v while changing the saturation. Therefore it is not enough to specify only the S w and v to determine v m , as dv/dS w depends on how the external parameters are controlled while varying S w . The parameters 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.

SUMMARY
We presented a detailed description of a set of algorithms for transporting fluids in a dynamic pore-network model of twophase 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.