Stable and efficient time integration of a dynamic pore network model for two-phase flow in porous media

We study three different time integration methods for a dynamic pore network model for immiscible two-phase flow in porous media. Considered are two explicit methods, the forward Euler and midpoint methods, and a new semi-implicit method developed herein. The explicit methods are known to suffer from numerical instabilities at low capillary numbers. A new time-step criterion is suggested in order to stabilize them. Numerical experiments, including a Haines jump case, are performed and these demonstrate that stabilization is achieved. Further, the results from the Haines jump case are consistent with experimental observations. A performance analysis reveals that the semi-implicit method is able to perform stable simulations with much less computational effort than the explicit methods at low capillary numbers. The relative benefit of using the semi-implicit method increases with decreasing capillary number $\mathrm{Ca}$, and at $\mathrm{Ca} \sim 10^{-8}$ the computational time needed is reduced by three orders of magnitude. This increased efficiency enables simulations in the low-capillary number regime that are unfeasible with explicit methods and the range of capillary numbers for which the pore network model is a tractable modeling alternative is thus greatly extended by the semi-implicit method.


Introduction
Different modeling approaches have been applied in order to increase understanding of immicsible two-phase flow in porous media. On the pore scale, direct numerical simulation approaches using e.g. the volume of fluid method [1] or the level-set oscillations have been made by introduction of modifications to the model. Medici and Allen [18] used a scheme where water was allowed to flow in the forward direction only in order to study water invasion in fuel cell gas diffusion layers. While this approach led to interesting results, it has some downsides. First, the interface movement is artificially restricted, and certain dynamic effects can not be resolved. This includes e.g. invasion front retraction after a Haines jump. Second, the method can only be used in cases with transient invasion. Studies of steady-state flow, such as those performed by Knudsen et al. [21] and Savani et al. [19], are not possible.
Because the oscillations originate in the numerical methods, rigorous attempts to remove them should focus on these methods rather than the models themselves. Joekar-Niasar et al. [9] followed this avenue and achieved stabilization using a linearized semiimplicit method. Their work, however, concerned a different type of pore network model than that considered here.
In this work, we present three numerical methods that can be utilized to perform stable simulations of two-phase flow in porous media with pore network models of the Aker type. The stability problems previously observed are thus solved without the need to resort to model modifications that restrict interface movement or preclude steadystate flow simulations. Two explicit methods are discussed, the forward Euler method and the midpoint method. These are stabilized by a new time step criterion derived herein. The third method is a new semi-implicit method. Thorough verifications of all methods are performed, confirming correct convergence properties and stability. Finally, we compare the methods in terms of performance.
The rest of this paper is structured as follows. Section 2 contains background information on the pore network model. Section 3 presents briefly the nomenclature, used in subsequent sections to describe the time integration methods. In Section 4, we recapitulate how the forward Euler method is used to integrate the model and we present a new time step criterion that stabilizes both forward Euler and the midpoint method at low capillary numbers. We briefly review the midpoint method in Section 5. The new semi-implicit method is described in detail in Section 6. Some remarks about the numerical implementation are made in Section 7. Section 8 contains a description of the cases simulated. Numerical experiments, including a Haines jump case, that show convergence and stability are given in Section 9 and a comparison of the method performances are made in Section 10. Section 11 summarizes and concludes the paper.

Pore network model
We consider incompressible flow of two immiscible fluids in a porous medium, where one fluid is more wetting towards the pore walls than the other. We call the less wetting fluid non-wetting (n) and the more wetting fluid we call wetting (w). The porous medium is represented in the model by a network of N nodes connected by M links. Each node is given an index i ∈ [0, N − 1], and each link is identified by the indices of the two nodes it connects. An example pore network is shown in Figure 1. The nodes are points that have no volume and, consequently, all fluid is contained in the links. The links therefore represent both the pore and the throat volumes of the physical porous medium. In this respect, the pore network model studied here differ from most other pore network models [11]. Each fluid is assumed to fill the entire link cross section. The interface positions are therefore each represented in the model by a single number, giving its location along the link length. The flow in each link is treated in a one-dimensional manner, where the flow is averaged over the link cross section. As we consider flow in relatively small cross sections only, we neglect any inertial effects and the volume flow rate (m 3 s −1 ) from node j to node i through the link connecting then is given by [20] q i j = −g i j z i j p i − p j − c i j z i j . (1) Herein, p i (Pa) is the pressure in node i, g i j (m 3 s −1 Pa −1 ) is the link mobility, c i j (Pa) is the link capillary pressure and z i j (m) is a vector containing the positions of any fluid interfaces present in the link. Both the link mobility and the capillary pressure depend on the fluid interface positions in the link. If two nodes i and j are not connected by a link, then g i j = 0. Due to mass conservation, the net flow rate into every node i is zero While the mobilities are symmetric with respect to permutation of the indices, the capillary pressures are anti-symmetric, Introducing this into (1), we obtain the immediately intuitive result The cross-sectional area of link i j is denoted a i j (m 2 ). Interface positions are advected with the flow according to when they are sufficiently far away from the nodes. Near the nodes, however, the interfaces are subject to additional modeling to account for interface interactions in the pores. This is discussed further in Section 2.3. The form of the expressions for the mobilities and capillary pressures depends on the shape of the links, and many different choices and modeling approaches are possible. Here, we will use models similar to those previously presented and used by e.g. Knudsen et al. [21] and Aker et al. [10]. However, the treated time integration methods are more general and can be applied to other models as well.

Link mobility model
We apply a cylindrical link model when computing the mobilities, so that Here, r i j (m) is the link radius and L i j (m) is the link length. The viscosity µ i j (Pa s) is the volume-weighted average of the fluid viscosities and can be computed from the wetting and non-wetting fluid viscosities µ w and µ n and the wetting fluid saturation s i j ,

Capillary pressure model
In each link i j, there may be zero, one or more interfaces present. These are located at the positions specified in z i j . As the interfaces may be curved, there may be a discontinuity in pressure at these interface locations. The capillary pressure c i j is the sum of interfacial pressure discontinuities in the link i j. When computing the capillary pressures, we assume that the links are wide near each end, and therefore that interfaces located near a link end have negligible curvature and no pressure discontinuity, while the links have narrow throats in the middle. The link capillary pressures are thus modeled as The interfacial tension between the fluids is denoted σ wn (N m −1 ) and The χ i j -function ensures zones of length αr i j at both ends of each link with zero capillary pressure across any interface located there. Choosing α = 0 is equivalent to replacing χ i j with z/L i j in (9).

Fluid interface interaction models
The equations discussed so far in this section describe how the fluids and the fluid interfaces move through the links. In addition, we rely on models for how they behave close to the nodes. The purpose of these are to emulate interface interactions in the pore spaces.
The following is assumed about the fluid behavior near the nodes and is accounted for by the fluid interface interaction models.
• The mass of each fluid is conserved at every node. This means that at all times, all wetting and non-wetting fluid flowing into a node from one subset of its neighboring links must flow out into another disjoint subset of its neighboring links.
• The network nodes in the model have no volume. However, due to the finite size of the physical pore void spaces, wetting fluid flowing into a pore space must be able to flow freely past any non-wetting fluid occupying the node point if the non-wetting fluid does not extend far enough into the pore void space cut the wetting fluid off. An example is illustrated in Figure 2. We consider a link i j to be cut off from free outflow of wetting fluid if the non-wetting fluid continuously extends a length at least αr i j into the link. Non-wetting fluid may freely flow past wetting fluid, or not, the same manner.
• In each link i j, interfacial tension will prevent droplets with length smaller than αr i j from forming by separation from larger droplets. An example is illustrated in Figure 3. Figure 3: (a) Small non-wetting bubble (blue) whose volume is small compared to the link volumes and is prevented from splitting by interfacial tension. This limits the minimum size of non-wetting bubbles, which will either (b) be stuck or (c) move through one of the links without splitting.

Boundary conditions
We consider only networks where the nodes and links can be laid out in the twodimensional x-y plane. These networks will be periodic in both the x-and y-direction. However, the model is also applicable to networks that extend in three dimensions [22], and the presented numerical methods are also compatible both with networks in three dimensions and with other, non-periodic boundary conditions [23].
We will here apply two types of boundary conditions to the flow. With the first type, a specified pressure difference ∆P (Pa) will be applied across the network in the y-direction. This pressure difference will be equal to the sum of all link pressure differences in any path spanning the network once in the y-direction, ending up in the same node as it started. With the other type of boundary condition, we specify a total flow rate Q (m 3 s −1 ) across the network. This flow rate will be equal to the sum of link flow rates flowing through any plane drawn through the network normal to the y-axis.

Temporal discretization
In the following three sections, we describe the different time integration methods considered. These methods are applied to (6), where evaluation of the right hand side involves simultaneously solving the mass conservation equations (2) and the constitutive equations (1) to obtain all unknown link flow rates and node pressures.
The discretized times (s) are denoted with a superscript where n is the time step number, The time step ∆t (i) is the difference between t (i+1) and t (i) and the time t (0) is the initial time in a simulation. Similarly, quantities evaluated at the discrete times are denoted with time step superscripts, e.g.
Mobilities and capillary pressures with superscripts are evaluated using the interface positions at the indicated time step,

Forward Euler method
The forward Euler method is the simplest of the time integration methods considered here and is the one used most frequently in previous works, see e.g. [21] and [24]. We include its description here for completeness and to provide context for the proposed new capillary time step criterion that is introduced to stabilize the method at low capillary numbers.
The ordinary differential equation (ODE) (6) is discretized in a straightforward manner for each link i j using forward Euler, The flow rates are calculated by inserting (1), evaluated with the current known interface positions, into the mass conservation equations (2). This results in the a system of linear equations consisting of one equation, for each node i with unknown pressure. This linear system can be cast into matrix form, where the vector x contains the unknown node pressures, e.g.
The matrix elements are and the elements of the constant vector are The node pressures are obtained by solving this linear equation system. The flow rates are subsequently evaluated using (16) and the interface positions are then updated using (15) and the interface interaction models.

Time step restrictions
In previous works [10,21], the time step length was chosen from a purely advective criterion, The parameter C a corresponds to the maximum fraction of a link length any fluid interface is allowed to move in a single forward Euler time step. The value of C a must be chosen based on the level of accuracy desired from the simulation. However, selecting the time step based on the advective criterion only, often results in numerical instabilities at low capillary numbers, where viscous forces are small relative to the capillary forces. This is demonstrated in Section 9.2. The origins of the numerical instabilities can be identified by performing analysis on a linearized version of the governing equations. This is done in Appendix A. This analysis also leads to a new time step criterion, whereby the time step length is restricted by the sensitivity of the capillary forces to perturbations in the current interface positions, For the particular choice of capillary pressure model given by (9), we obtain According to the linear analysis, numerical instabilities are avoided if the parameter C c is chosen such that 0 < C c < 1. However, we must regard (23) as an approximation when we apply it to the full non-linear model in simulations and, consequently, we may have to chose C c conservatively to ensure stability for all cases. At each step in the simulation, the time step used is then taken as to comply with both the advective and the capillary time step criteria. The capillary time step restriction (23) is independent of flow rate. It therefore becomes quite severe, demanding relatively fine time steps, when flow rates are low.

Boundary conditions
The periodic boundary conditions, specifying a total pressure difference ∆P across the network, can be incorporated directly into the linear equation system (18). For each node i, a term g (n) i j ∆P is added to or subtracted from b i for any link i j that crosses the periodic boundary.
With the specified ∆P condition implemented, we can use it to obtain the node pressures and link flow rates corresponding to a specified total flow rate Q. Due to the linear nature of the model, the total flow rate is linear in ∆P [10], so that for some unknown coefficients C 1 and C 2 , that are particular to the current fluid configuration. We choose two different, but otherwise arbitrary, pressure drop values ∆P 1 and ∆P 2 and, using the above procedure, we solve the network model once for each pressure difference and calculate the corresponding total flow rates Q 1 and Q 2 . The coefficients C 1 and C 2 are then determined by, The pressure difference ∆P required to obtain the specified flow rate Q is determined by solving (26) for ∆P. Subsequently, the network model is solved a third time with pressure drop ∆P to obtain the desired node pressures and link flow rates.

Midpoint method
The forward Euler method is first-order accurate in time. To obtain smaller numerical errors, methods of higher order are desirable. We therefore include in our discussion the second-order midpoint method. This method is identical to that used by Aker et al. [10], except with respect to choice of time step length.
The ODE (6) is discretized as where q (n+1/2) i j is the flow rate at the midpoint in time between point n and n + 1. This flow rate is calculated in the same manner as described in Section 4. The interface positions at n + 1/2 are obtained by taking a forward Euler step with half the length of the whole time step,

Time step restrictions
Since the forward Euler stability region is contained within the stability region for the midpoint method, we use the same time step restrictions for the midpoint method as for forward Euler, see Section 4.1.

Boundary conditions
Both the specified ∆P and the specified Q boundary conditions are incorporated into the midpoint method by applying the procedures described in Section 4.2 for each evaluation of the right hand side of (6).

Semi-implicit method
To avoid both the numerical instabilities and the time step restriction (23), which becomes quite severe at low flow rates, we here develop a new semi-implicit time stepping method. Simulation results indicate that this method is stable with time steps determined by the advective criterion (22) only, and much longer time steps are therefore possible than with the forward Euler and midpoint methods at low capillary numbers.
The ODE (6) is now discretized according to The semi-implicit nature of this discretization comes from the flow rate used, Herein, the link mobility is evaluated at time step n, while the node pressures and the capillary pressure are evaluated time step n + 1.
The link mobilities could of course also have been evaluated at time step n + 1, thus creating a fully implicit backward Euler scheme. As is shown in Appendix A, we may expect backward Euler to be stable with any positive ∆t (n) . The backward Euler scheme may therefore seem like a natural choice for performing stable simulations with long time steps. However, to evaluate the mobilities at time step n + 1 complicates the integration procedure and was found to be unnecessary in practice. A semi-implicit alternative is therefore preferred.
To obtain the node pressures, we solve the mass conservation equations, Again, we have one equation for each node i with unknown pressure. However, because the capillary pressures now depend on the flow rates, the mass conservation equations are now a system of non-linear equations, rather than a system of linear equations. This system can be cast in the form where x contains the unknown pressures, e.g.
In order to solve (35) using the numerical method described in Section 7, it is necessary to have the Jacobian matrix of F. Details on how the Jacobian matrix is calculated are given in Appendix B.
The calculation of link flow rates from node pressures, and thus every evaluation of F and its Jacobian, involves solving one non-linear equation for each link flow rate, The procedure for updating the interface positions with the semi-implicit method may be summarized as follows. The non-linear equation system (35) is solved to obtain the unknown node pressures. In every iteration of the solution procedure, the flow rates are evaluated by solving (37) for each link. When a solution to (35) is obtained, the interface positions are updated using (31) and the interface interaction models.

Time step restrictions
We aim to select the time steps such that However, to solve the non-linear system (35) is challenging in practice and requires initial guess values for the link flow rates and node pressures that lie sufficiently close to the solution. For this purpose, we here use values from the previous time step. This turns out to be a sufficiently good choice for most time steps, but our numerical solution procedure does not always succeed. As the link flow rates and node pressures at two consecutive points in time become increasingly similar as the time interval between them is reduced, we may expect the guess values to lie closer to the solution if we reduce the time step. Thus, if our solution procedure is unable to succeed, our remedy is to shorten ∆t (n) . This will sometimes lead to time steps shorter than ∆t (n+1) a . If, for a given time step, ∆t (n) must be reduced to less than twice the time step length allowed by the explicit methods, we revert to forward Euler for that particular step. As we demonstrate in Section 10, however, this does not prevent the semi-implicit method from being much more efficient than the explicit methods at low capillary numbers.

Boundary conditions
As with the explicit methods, the specified ∆P boundary condition can be incorporated directly into the mass balance equation system, in this case (35). This is done by adding to or subtracting from the right hand sides of (32) and (37) a term g (n) i j ∆P for each link i j crossing the periodic boundary.
The specified flow rate boundary condition is incorporated by including ∆P as an additional unknown and adding an additional equation to the non-linear equation system (35). Herein, Ω is the set of links crossing the periodic boundary, with i being the node on the downstream side and j being the node on the upstream side. Thus, (40) is satisfied when the total flow rate through the network is equal to Q.

Implementation
The non-linear equation system (35) is solved using a Newton-type solution method that guarantees convergence to a local minimum of F · F, see Press et al. [25] pp. 477. However, a local minimum of F · F is not necessarily a solution to (35), and good initial guess values for the node pressures and link flow rates are therefore crucial. For this purpose, we use the values from the previous time step and reduce the length of the current time step if the solution method fails, as discussed in Section 6.1.
Solving (37) is done using a standard Newton solver [26]. For robustness, a bisection solver [26] is used if the Newton solver fails.
The Newton-type solver for non-linear systems and the explicit time integration methods require methods for solving linear systems of equations. We use the conjugate gradient method in combination with the LU preconditioner implemented in the PETSc library, see Balay et al. [27]. An introduction to solving systems of Kirchhoff-type equations numerically can be found in [28].

Case descriptions
In this section, we describe the two simulated cases. One is a test case where a single bubble is contained in a network consisting of links connected in series, while the other is designed to capture a single Haines jump in a small network where fluids flow at a specified rate.

Links-in-series test case
The verification will include comparison of results from the various numerical methods applied to a test case. The test case is chosen such that it can be set up as a single ODE with a closed expression for the right-hand side. An accurate reference solution can thus be easily obtained using a high-order Runge-Kutta method. As our  [30] test case, we consider a network consisting of M = 3 identical links connected in series. The network contains a single bubble of length (m) with center position z (m).
In the capillary pressure model, we choose α = 0. The ODE (6) can then be restated as an equivalent equation for the bubble position, where Q is the flow through the network and a is the link cross-sectional area. The model equations can be reduced to the following expression for flow rate.
Here, g is the mobility of a single link, L = 1.0 · 10 −3 m is the length of a single link and r = 1.0 · 10 −4 m is the link radius. The bubble has length = 4.8 · 10 −4 m and is initially located at z = 2.4 · 10 −4 m. The fluid parameters used in all simulations are given in Table 1. The pressure difference ∆P will be stated for each simulation.

Haines jump case
The Haines jump was first reported almost 90 years ago [13]. It refers to the sudden drops in driving pressure observed in drainage experiments when non-wetting fluid breaks through a throat and invades new pores. This process was studied experimentally and numerically by Måløy et al. [16] and, more recently, it was imaged directly and analyzed in detail by Armstrong and Berg [15] for flow in a micromodel and by Berg et al. [14] for flow in a sample of Berea sandstone. The Haines jump case simulated here captures one such break-through and subsequent pressure drop.
Among the findings in the study by Måløy et al. [16] was that pore drainage is a non-local event, meaning that as one pore is drained, imbibition occurs in nearby neck regions. This was also observed by Armstrong and Berg [15], and was explained as follows. When the imposed flow rates are low, the non-wetting fluid that fills the newly invaded pores needs to be supplied from nearby locations rather than the external feed. Armstrong and Berg [15] also found, for their range of investigated parameters, that pore drainage occurred on the same time-scale, regardless of the externally imposed flow rate.
We consider a hexagonal network consisting N = 24 nodes and M = 36 links. All links have length 1.0 · 10 −3 m, while the link radii are drawn randomly from a uniform distribution between 0.1 and 0.4 link lengths. In the capillary pressure model, we choose α = 1. The fluid parameters µ w , µ n and σ wn are the same as in the links-inseries test case, see Table 1. With these fluid parameters and network length scales, the case mimics the flow of water (w) and decane (n) in a Hele-Shaw cell filled with glass beads similar to those used in e.g. [16,31,32]. The linear dimensions are ∼ 10 times bigger in this network compared to the micromodel of Armstrong and Berg [15]. Initially, the fluids are distributed in the network as shown in Figure 4, with the nonwetting fluid in a single connected ganglion.
Simulations are run at different specified flow rates Q until a net fluid volume equivalent to 5 % of the total pore volume has flowed through the network. The flow dynamics will, of course, depend upon the specified flow rate. At low flow rates, however, the flow will exhibit some relatively fast fluid redistribution events and one relatively slow pressure build-up and subsequent Haines jump event. The Haines jump will occur as the non-wetting fluid breaks through the link connecting nodes 9 and 16 and invades node 16, see Figure 4.
It was mentioned by Armstrong and Berg [15] that the large local flow velocities that they observed as a pore was filled with non-wetting fluid during a Haines jump has implications for how such processes must be numerically simulated. Specifically, the time resolution of the simulation needs to be fine enough during these events to capture them. This poses a challenge when externally applied flow rates are low and there is thus a large difference in the large time scale that governs the overall flow of the system and the small time scale than governs the local flow during Haines jumps.

Verification
In this section, we verify that the time integration methods presented correctly solve the pore network model equations and that the time step criteria presented give stable solutions.

Convergence tests
All time integration methods presented should, of course, give the same solution for vanishingly small time steps. Furthermore, the difference between the solution obtained with a given finite time step and the fully converged solution should decrease as the time steps are refined, and should do so at a rate that is consistent with the order of the method. In this section, we verify that all three time integration methods give solutions that converge to the reference solution for the links-in-series test case and thus that the methods correctly solve the model equations for this case.
We choose the pressure difference to be ∆P = −3200 Pa. This value is large enough to overcome the capillary forces and push the non-wetting bubble through the links. We therefore expect a flow rate Q that varies in time, but is always positive.
As measures of the numerical error, we consider both the relative error in the flow rate Q and the relative error in bubble position z between the numerical solutions and reference solutions at the end of the simulation. Time integration is performed from t = 0 s to t = 0.001 44 s. To have control over the time step lengths, we ignore all time step criteria for now and instead set a constant ∆t for each simulation.
In Figure 5, flow rates are plotted for each of the time integration methods. Results using a coarse time step, ∆t = 4 · 10 −5 s, and a fine time step, ∆t = 1 · 10 −5 s, are shown along with the reference solution.
For the forward Euler and the semi-implicit method, there is considerable discrepancy between the numerical and the reference solution with the coarse time step. The flow rate obtained from forward Euler lags behind the reference solution, while that from the semi-implicit method lies ahead of it. This may be expected, however, since forward Euler at each time step uses current information in the right hand side evaluation, whereas the semi-implicit method uses a combination of current and future information. With the fine time step, there is less difference between the reference and the numerical solutions. With the more accurate midpoint method, the coarse-stepped numerical solution lies only marginally ahead of the reference solution while there is no difference between the fine-stepped numerical solution and the reference solution at the scale of representation.
The convergence of the numerical solutions to the reference solution upon time step refinement is quantified in Table 2, Table 3 and Table 4. Herein, the numerical errors and estimated convergence orders are given for the forward Euler, midpoint and semi-implicit method, respectively. For all methods considered, the numerical errors decrease when the time step is refined and do so at the rate that is expected. The forward Euler and the semi-implicit method exhibit first-order convergence, while the midpoint method shows second-order convergence. We note that the errors in both z and Q are similar in magnitude for the forward Euler and the semi-implicit method. The errors obtained with the midpoint method are smaller. The difference is one order of magnitude for ∆t = 1 · 10 −5 s.

Stability tests
In this section, we demonstrate that the proposed capillary time step criterion (23) stabilizes the forward Euler method and the midpoint method at low flow rates. We simulated two different cases and varied C c . Simulations run with low C c turned out to be free of spurious oscillations, indicating that the proposed criterion stabilizes the methods, while simulations run with C c significantly larger than unity produced oscillations, indicating that the proposed criterion is not unnecessarily strict.
First, consider the links-in-series test case with ∆P = 0 Pa. With no applied pressure difference, the flow is driven purely by the imbalance of capillary forces on the non-wetting bubble. Therefore, there should only be flow initially and the bubble should eventually reach an equilibrium position where both interfaces experience the same capillary force and the flow rate is zero. Simulations were run with C a = 0.1 and C c equal to 2.0, 1.0 and 0.5. Results from forward Euler are shown in Figure 6a and results from the midpoint method are shown in Figure 6b. In both figures, the reference solution is also shown.
The forward Euler results are stable and qualitatively similar to the reference solution with C c = 0.5. With C c = 1.0, there are some oscillations initially that are dampened and eventually vanish. From comparison with the reference solution, it is clear that such oscillations have no origin in the model equations and are artifacts of the numerical method. With C c = 2.0, the oscillations are severe and do not appear to be dampened by the method. Instead the non-wetting bubble keeps oscillating around its equilibrium position in a manner that is clearly unphysical.
The results from the midpoint method in Figure 6b follow a qualitatively similar trend as those from forward Euler with regard to stability. Results computed with C c = 0.5 are stable and results with C c = 2.0 exhibit severe oscillations. Still, the results from the midpoint method lie much closer to the reference solution than the results from the forward Euler method, as we would expect since the midpoint method is second-order. Both methods are, however, unstable with C c = 2.0, indicating that the while the midpoint method has improved accuracy over forward Euler, it is unable to take significantly longer time steps without introducing oscillations. This is consistent with the analysis in Appendix A, since the two methods have identical stability regions in real space.
Next, consider the Haines jump case with Q = 10 −9 m 3 s −1 , corresponding to Ca = 1.2 · 10 −5 . This case was run using the forward Euler method, C a = 0.1 and three different values of C c , equal to 4.0, 2.0 and 1.0. The required pressure difference to drive the flow at the specified rate is shown in Figure 7a. Figure 7b shows the pressure from Figure 7a in greater detail.
For all three values of C c , the main qualitative features of the flow are captured. We observe short transient pressure drops at t ≈ 0.08 s and t ≈ 0.20 s. These correspond to fluid redistribution events on the upstream side of the non-wetting ganglion, where the fluid rearranges itself to a more stable configuration with little change to the interface positions on the downstream side. The event at t ≈ 0.20 s is illustrated in Figure 8. The fluid redistribution is driven by capillary forces and less external pressure is therefore required to drive the flow during these events.
We also observe the slow pressure build-up from t ≈ 0.10 s to t ≈ 0.23 s, when the driving pressure becomes large enough to overcome the capillary forces and cause   break-through of non-wetting fluid in the link connecting nodes 9 and 16, and we observe the subsequent Haines jump. The fluid configurations before and after the Haines jump are shown in Figure 9. Notice also that non-wetting fluid at the downstream end of the moving ganglion retracts during the Haines jump in links near to where the break-through occurs. This is seen e.g. in the links downstream of nodes 10 and 14.
That such local imbibition occurs near the drained pore is in agreement with the observations of Armstrong and Berg [15], and shows that the model is able to capture the non-local nature of pore drainage events in a numerically stable manner when the new numerical methods are used. As in the links-in-series case, the solution exhibits oscillations for the values of C c that are larger than unity. With C c = 1.0, the results are free from oscillations and appear stable. This indicates that the stability criterion (23) is valid and not unnecessarily strict also for a network configuration that is much more complex than links in series.
Both the links-in-series case and the Haines jump case were simulated with the semi-implicit method and produced stable results with the advective time step criterion (22) only. The results from the links-in-series test case are shown in Figure 6c. For brevity, the results from the Haines jump case are omitted here. The reader is referred to Figure 10a in Section 10, where stable results are shown for a lower flow rate.
To summarize, both the forward Euler and midpoint methods produce stable results for the cases considered when the capillary time step criterion (23) is used in addition to (22) to select the time step lengths. By running simulations with different C c , we have observed a transition from stable to unstable results for values of C c near 1, in order of magnitude. In the Haines jump case, all methods presented are able to capture both the fast capillary-driven fluid redistribution events, and the slow pressure build-up before a Haines jump.

Performance analysis
In this section, we analyze and compare the performance of the time integration methods. In doing so, we consider the number of time steps and the wall clock time required to perform stable simulations of the Haines jump case with each of the methods at different specified flow rates Q. The flow rates simulated were 10 −7 m 3 s −1 , 10 −8 m 3 s −1 , 10 −9 m 3 s −1 , 10 −10 m 3 s −1 , 10 −11 m 3 s −1 and 10 −12 m 3 s −1 . The accuracy of the methods was studied Section 9.1, and will not be part of the performance analysis. Instead, stable simulations are considered sufficiently accurate.
First, we look more closely at the results for Q = 10 −11 m 3 s −1 , corresponding to Ca = 1.2 · 10 −7 . The pressure difference required to drive the flow is shown in Figure  10a, and the time step lengths used are shown in Figure 10b. From the latter Figure, we see that the semi-implicit method is able to take longer time steps than forward Euler for most of the simulation. During the pressure build-up phase, the difference is four orders of magnitude. During the fast capillary-driven fluid redistribution events, however, the length of the semi-implicit time steps drop to the level of those used by forward Euler. This is because we here have relatively large flow rates in some links, even though Q is low, and the advective time step criterion (22) becomes limiting for both the semi-implicit method and forward Euler.
It was mentioned by Armstrong and Berg [15] that any accurate numerical simulation on the pore scale must have a time resolution fine enough to capture the fast events. The semi-implicit method accomplishes this by providing a highly dynamic time resolution, which is refined during the fast events. The method is therefore able to resolve these events, while time resolution can be coarsened when flow is governed by the slow externally applied flow rate, saving computational effort.
The time duration of the Haines jump pressure drops for all except the two largest externally applied flow rates were around 10 ms. This is in qualitative agreement with the results presented by Armstrong and Berg [15]. They found that, for their investigated range of parameters, pores were drained on the millisecond time scale regardless of externally applied flow rate. However, we stress that although we consider the same fluids, the pore network used here was approximately one order of magnitude larger in the linear dimensions than that of [15].
The number of time steps and wall clock time required to simulate the Haines jump case at different specified flow rates Q are shown in Figure 11a and Figure 11b, respectively.
For the explicit methods, both the number of time steps and the wall time are proportional to Ca −1 at low capillary numbers. This is because the capillary time step criterion (23) dictates the time step at low capillary numbers (except during fast fluid redistribution events). The criterion depends on the fluid configuration, while it is independent of the flow rate. At low enough flow rates, the system will pass through roughly the same fluid configurations during the simulation, regardless of the applied Q. The speed at which the system passes through these configurations, however, will be inversely proportional to Q and therefore, so will the required wall time and number of time steps. As the forward Euler and the midpoint method are subject to the same time step criteria, these require roughly the same number of time steps at all considered flow rates. However, since the midpoint method is a two-step method, the wall time it requires is longer and approaches twice that required by the forward Euler for long wall times. For the semi-implicit method, on the other hand, the number of time steps required to do the simulation becomes effectively independent of the specified flow rate at capillary numbers smaller than approximately 10 −4 . The result is that low-capillary number simulations can be done much more efficiently than with the explicit methods, in terms of wall time required to perform stable simulations. This is seen in Figure 11b. At Ca ∼ 10 −5 , the computational time needed by all three methods are similar in magnitude. The relative benefit of using the semi-implicit method increases at lower capillary numbers. For the lowest capillary number considered, the difference in wall time between the explicit methods and the semi-implicit is three orders of magnitude.
The increased efficiency of the semi-implicit method over explicit methods at low capillary numbers means that one can use the semi-implicit method to perform simulations in the low capillary number regime that are unfeasible with explicit methods. Thus, the range of capillary numbers for which the pore network model is a tractable modeling alternative is extended to much lower capillary numbers. This includes e.g. simulations of water flow in fuel cell gas diffusion layers, where capillary numbers are can be 10 −8 [33].
Finally, to study the effect of an increase in network size on the wall time required by the semi-implicit method, the Haines jump case was run on three scaled-up versions of the network with N = 24 nodes considered so far, illustrated in Figure 4. All simulations were run at Ca ∼ 10 −7 . In Figure 12 the wall clock time time required is plotted against the number of nodes N for the different networks. The wall time is seen to increase proportionally to N 2 .

Conclusion
We have studied three different time integration methods for a pore network model for immiscible two-phase flow in porous media. Two explicit methods, the forward Euler and midpoint methods, and a new semi-implicit method were considered. The explicit methods have been presented and used in other works [10,21,24], and were reviewed here for completeness. The semi-implicit method was presented here for the first time, and therefore in detail.
The explicit methods have previously suffered from numerical instabilities at low capillary numbers. Here, a new time-step criterion was suggested in order to stabilize them and numerical experiments were performed demonstrating that stabilization was achieved.
It was verified that all three methods converged to a reference solution to a selected test case upon time step refinement. The forward Euler and semi-implicit methods exhibited first-order convergence and the midpoint method showed second-order convergence.
Simulations of a single Haines jump were performed. These showed that the all three methods were able to resolve both pressure build-up events and fluid redistribution events, including interfacial retraction after a Haines jump, which may occur at vastly different time scales when capillary numbers are low. The results from the Haines jump case were consistent with experimental observations made by Armstrong and Berg [15]. Fluid redistribution events cannot be properly captured when using solution methods that have previously been used at low capillary numbers that e.g. do not allow backflow [18].
A performance analysis revealed that the semi-implicit method was able to perform stable simulations with much less computational effort than the explicit methods at low capillary numbers. For the case considered, the computational time needed was approximately the same for all three methods at Ca ∼ 10 −5 . At lower capillary numbers, the computational time needed by the explicit methods increased inversely proportional to the capillary number, while the time needed by the semi-implicit method was effectively constant. At Ca ∼ 10 −8 , the computational time needed by the semi-implicit methods was therefore three orders of magnitude smaller than those needed by the explicit methods.
The superior efficiency of the new semi-implicit method over the explicit methods at low capillary numbers enables simulations in this regime that are unfeasible with explicit methods. Thus, the range of capillary numbers for which the pore network model is a tractable modeling alternative is extended to much lower capillary numbers. This includes e.g. simulations of water flow in fuel cell gas diffusion layers, where capillary numbers are can be 10 −8 [33].
In summary, use of Aker-type pore network models were previously restricted to relatively high capillary numbers due to numerical instabilities in the explicit methods used to solve them. With the new time step criterion presented here, these stability problems are removed. However, simulations at low capillary numbers still take a long time and the computational time needed increases inversely proportional to the capillary number. This problem is solved by the new semi-implicit method. With this method, the computational time needed becomes effectively independent of the capillary number, when capillary numbers are low.

Author contributions
To pursue low capillary number simulations with Aker-type pore network models was proposed by SK, MG and AH in collaboration. MV developed the particular variation of the pore network model used. MG developed the new numerical methods and performed the simulations. MG wrote the manuscript, aided by comments and suggestions from MV, SK and AH.
Consider a single link i j in a network and assume that p i and p j are given. Then the ODE (6) for the interface positions in the link is We further assume that the flow rate in this link is low. This means that the node and capillary pressures almost balance at the current interface positions z * i j , and thus q i j z * i j ≈ 0. Also, we neglect the dependence of g i j on the interface positions. Now rewrite (A.1) in terms ∆z i j = z i j − z * i j and linearize the right hand side around z * i j to get We can now read off the approximate ODE eigenvalue as Without loss of generality, we may assume that λ < 0. If this is not the case, we interchange the indices i and j and redefine our spatial coordinate so that z → −z to get an ODE with negative λ. We therefore write the eigenvalue as If the forward Euler method is to be stable on the linearized ODE, λ∆t must lie in the stability region of the forward Euler method [26], The criterion (23) is obtained by demanding that (A.8) be satisfied for all links in the network. If the advective criterion (22) is used by itself and the link flow rates are low, then (A.8) is not necessarily satisfied for all links and we must expect numerical instabilities from the forward Euler method.
As the midpoint method has the same real-space stability region as the forward Euler method (A.7), the above reasoning and the criterion (23) can be applied for the midpoint method also.
The backward Euler method, on the other hand, is stable if [26] λ∆t < 0, (A. 9) and, because λ is negative, it is stable with any positive ∆t for this linearized problem.
Appendix B. Jacobian matrix for the semi-implicit method In order to solve (35) using the numerical method described in Section 7, it is necessary to have the Jacobian matrix of F. This matrix may be written as for the specific choice of capillary pressure model given by (9).
As the pore network model is linear in the node pressures, it is intuitive that the effect on the link flow rate of increasing the pressure in the node at one end of a link is the same as decreasing it, by the same amount, in the node at the other end. Thus we may write ∂q (n+1) The addition of F m (40) to the non-linear system for the specified flow rate boundary condition, introduces some new terms in the Jacobian matrix of F. The derivatives of F m with respect to the node pressures are The additional terms corresponding to the derivatives with respect to ∆P of the mass balance equations for each node k with unknown pressures are (B.12)