Pore network modeling of the effects of viscosity ratio and pressure gradient on steady-state incompressible two-phase flow in porous media

We perform more than 6000 steady-state simulations with a dynamic pore network model, corresponding to a large span in viscosity ratios and capillary numbers. From these simulations, dimensionless quantities such as relative permeabilities, residual saturations, mobility ratios and fractional flows are computed. Relative permeabilities and residual saturations show many of the same qualitative features observed in other experimental and modeling studies. However, while other studies find that relative permeabilities converge to straight lines at high capillary numbers we find that this is not the case when viscosity ratios are different from 1. Our conclusion is that departure from straight lines occurs when fluids mix rather than form decoupled flow channels. Another consequence of the mixing is that computed fractional flow curves, plotted against saturation, lie closer to the diagonal than they would otherwise do. At lower capillary numbers, fractional flow curves have a classical S-shape. Ratios of average mobility to their high-capillary number limit values are also considered. These vary, roughly, between 0 and 1, although values larger than 1 are also observed. For a given saturation and viscosity ratio, the mobilities are not always monotonically increasing with the pressure gradient. While increasing the pressure gradient mobilizes more fluid and activates more flow paths, when the mobilized fluid is more viscous, a reduction in average mobility may occur.


Introduction
A number of different modeling approaches have been applied to study two-phase flow in porous media. These include direct numerical simulations (DNS), which employ, e.g., the volume-of-fluid method (Raeini et al. 2012) or the level-set method (Jettestuen et al. 2013;Gjennestad and Munkejord 2015) to keep track of the fluid interfaces, lattice-Boltzmann methods (Ramstad et al. 2012;Armstrong et al. 2016) and pore network models. Recently, several of these methods were compared in a benchmark study 1 3 by Zhao et al. (2019), where participants were asked to reproduce experimentally studied transient fluid displacement processes at different capillary numbers and wettability conditions. This benchmark study, and the bulk of works in the literature, focuses on transient processes. Less attention has been given to pore-scale modeling and experiments on steady-state flow. In steady state, fluids may flow and interfaces move on the pore scale. However, the total flow rates and average pressure drop for a sufficiently large representative elementary volume fluctuate around well-defined time-averaged values (Erpelding et al. 2013;Rücker et al. 2015;Hansen et al. 2018). Our focus here is on such states and their corresponding time-averaged steady-state properties.
On the modeling side, part of the explanation for the limited focus on steady state is probably that steady-state simulations require longer simulation times compared to transient processes. While breakthrough of the invading phase typically happens for simulation times corresponding to much less than one pore volume of flow in transient cases, several pore volumes may be required to obtain decent time averages of steady-state quantities.
In spite of this, some studies on steady-state two-phase flow have been carried out. Avraam and Payatakes (1995) performed quasi-2D micro-model experiments, varied the capillary number, the viscosity ratio and the flow rate ratio, and found four different flow regimes. They also studied relative permeabilities. Steady-state simulations with a pore network model of the Aker type (Aker et al. 1998b) have also been performed by , , Ramstad and Hansen (2006), Tørå et al. (2009), Sinha et al. (2017) and Sinha et al. (2019b). In particular,  performed simulations with equal viscosities and one value for the interfacial tension, and studied the effect of changing total flow rate on fractional flow and relative permeabilities. Results for equal viscosities are interesting and applicable in some cases, e.g., for oil and water (Oak et al. 1990). For other applications, e.g., sequestration of supercritical CO 2 (Bennion and Bachu 2005) and gas-liquid flows such as in fuel cells, viscosity contrast should be taken into account.
The aim of this work is to shed light on how different steady-state flow properties behave as pressure gradients are increased from values corresponding to moderate capillary numbers around 10 −3 -10 −4 to the high capillary number limit. Furthermore, we aim to assess the impact of viscosity ratio in this context. To this end, we perform steady-state simulations using a dynamic pore network model of the Aker type (Aker et al. 1998b;Sinha et al. 2019a) to represent a block of porous material. In each simulation, the time evolution of the fluid configurations in the network is resolved, yielding a time series of total flow rates and average pressure drops for the entire network. These time series are subsequently timeaveraged to obtain steady-state values that are used to calculate quantities such as relative permeabilities, fractional flow and capillary number. We give results from more than 6000 such steady-state simulations that cover a large range of viscosity ratios and capillary numbers. To aid further research, the simulation data are published along with this article.
In the simulations, we utilize a new time step criterion in the numerical solution method  to perform numerically stable simulations at low and moderate capillary numbers. The new methodology has an important effect on capillary numbers below 10 −3 . In addition, we make extensive use of results from a recent study of the high capillary number regime (Sinha et al. 2019b) in the analysis of the results.
The discussion is restricted to capillary numbers above 10 −4 , where history-dependence of the steady-state quantities is negligible Erpelding et al. 2013). At lower capillary numbers, steady-state quantities are harder to define and calculate. To allow for a discussion which is as general as possible, and which allows for comparison with other studies of slightly different systems, we focus on dimensionless steady-state quantities, such as relative permeabilities, mobility ratios and fractional flow. 1 The chosen pore network model is of the Aker type (Aker et al. 1998b). It has several properties that are advantageous when computing steady-state quantities. First, it is dynamic and thus captures the effects of both viscous and capillary forces. Second, it can be solved in a numerically stable manner at arbitrarily low capillary numbers ). Third, it is possible to apply periodic boundary conditions, keeping the saturation constant and eliminating effects of saturation gradients. Furthermore, it is computationally cheap, making the study of large enough systems over long enough times possible.
In spite of these advantages, however, the model also has some limitations. In particular, film flow is not accounted for. An extension of the model that includes film flow has been developed (Tørå et al. 2012), but it is computationally more demanding than the present model and is therefore not used here. While film flow effects could, in principle, also be captured, e.g., by a DNS or lattice-Boltzmann simulations, very high spatial resolution is required to resolve such films properly (Zhao et al. 2019). This makes such an approach prohibitively expensive for steady-state calculations, especially when a large number of them are desired.
The rest of the paper is structured as follows. In Sect. 2, we describe the system under consideration, define some important steady-state flow properties and discuss the high capillary number limit. In Sect. 3, we describe the pore network model used, the numerical methods used to solve it and the procedures used to obtain steady-state time averages in some detail. Results are presented and discussed in Sect. 4, and concluding remarks are given in Sect. 5.

Steady-State Flow
In this section, we define the system under study, some quantities that will be used to describe steady-state flow and discuss the high capillary number limit.
The system we consider is a block of porous material, as illustrated in Fig. 1. It has cross-sectional area A and thickness Δx in the direction of flow (the x-direction). The volume of the block is The pore space volume in the block is V p , so that the porosity is The pore space is occupied by two fluids, where one is more wetting toward the pore walls than the other. In the following, we will call the more wetting fluid wetting ( w ) and the less wetting fluid non-wetting ( n ). The fluids are assumed to be incompressible and S w is the wetting fluid saturation, i.e., the fraction of the pore space volume occupied by the wetting fluid.
1 These quantities are used to provide a familiar framework of dimensionless quantities in which results are presented and discussed. However, other quantities could also, in principle, be used to convey the same information. One example is the velocities presented by Hansen et al. (2018).
A pressure difference, either constant or fluctuating, exists across the porous block. This causes the wetting and non-wetting fluids to flow at time-dependent rates Q w (t) and Q n (t) , respectively. In a steady-state Q w and Q n may fluctuate, but do so around a well-defined time-averaged value (Erpelding et al. 2013;Rücker et al. 2015;Hansen et al. 2018). We shall be most concerned with such time averages, which we call the steady-state flow properties. The steady-state wetting flow rate is defined as where t 1 , t 2 is the time period of averaging. The steady-state non-wetting flow rate and the steady-state pressure difference Δp are defined analogously. 2 The steady-state total flow rate is, and the fractional flow of wetting fluid is The volume-averaged fluid velocity in the pore space, the seepage velocity, is and the average mobility is In the literature, description of two-phase flow in porous media is typically done in terms of the relative permeabilities r w and r n . These are related to the steady-state properties Q w , Q n and Δp by Herein, w and n are the viscosities of the wetting and non-wetting fluids, respectively, and is the absolute permeability of the porous medium.
When measuring relative permeabilities (Oak et al. 1990;Bennion and Bachu 2005) and when using relative permeability models to do continuum-scale calculations, it is often only their dependence on S w which is considered. It is, however, well established that variation with viscosity ratio M = n ∕ w and capillary number Ca cannot, in general, be neglected (Avraam and Payatakes 1995;Bardon and Longeron 1980;Datta et al. 2014;Armstrong et al. 2016;Guo et al. 2015).

The High Capillary Number Limit
Even though the relative permeabilities are dependent on capillary number, there seems to be general agreement in the literature that this dependence disappears at high capillary numbers (Bardon and Longeron 1980;Avraam and Payatakes 1995;Whitson et al. 2003;Ramstad et al. 2012;Schechter and Haynes 1992;Sinha et al. 2019b). We call this the high-Ca limit.
In particular, Sinha et al. (2019b) studied the high-Ca limit obtained when the two fluids are driven through the porous medium by a large pressure gradient while retaining their immiscibility. For a pore network model and 2D lattice-Boltzmann simulations, they showed that the fluid velocity in this limit could be described by with an effective viscosity The exponent depended on the degree of intermixing of the fluids, induced by the flow through the porous medium. Intermixing here means that the fluids flow as many small disconnected structures when well intermixed, while they flow as larger structures when less well intermixed. For highly intermixed flows containing small droplets, it was found that = 1 , while for less well intermixed flows was shown to decrease. For a flow regime with lubrication layers, where the two fluids flowed in parallel through the pores, = −0.5 was obtained.
These results for may be illustrated by considering the following two examples. First, an exponent = −1 would be obtained if the porous medium were modeled by a bundle of identical capillary tubes, where each tube contains only one or the other of the two fluids. In this case, the fluids are not intermixed at all and their flows are completely decoupled. Second, consider the case where the tubes in the model each contain a train of bubbles, with the same saturation in each tube. In this case, the fluids must be considered well intermixed and there is strong coupling between their flows. They are in fact forced to have same velocity. The corresponding exponent would be = 1 . The degree of intermixing between the flow of the two fluids, resulting in the coupling between their flows being strong or weak, thus determines . In a real porous medium and in a network model, one might expect partial coupling of the flows and an exponent somewhere in between −1 and 1. This was indeed what was found by Sinha et al. (2019b). For the pore network model and porous medium studied here, it was shown that = 0.6.
Many other works in the literature refer to the more general idea of a high-Ca limit (Avraam and Payatakes 1995;Whitson et al. 2003;Ramstad et al. 2012;Armstrong et al. 2016). However, the prevailing view seems to be that the relative permeabilities approach straight lines at high capillary numbers (Bardon and Longeron 1980;Avraam and Payatakes 1995;Whitson et al. 2003;Ramstad et al. 2012;Schechter and Haynes 1992), i.e., To show how straight-line relative permeabilities are compatible with the results from Sinha et al. (2019b), we introduce (8), (9), (12) and (13) and into (6) and express the flow velocity in the high-Ca limit as This result is a special case of (10) with = −1.
Some of the experimental studies carried out that obtain straight-line relative permeabilities consider two phases near the mixture's critical point, see, e.g., Bardon and Longeron (1980); Schechter and Haynes (1992). In these cases, the interfacial tension approaches zero, and the capillary number infinity, as the phases are brought closer to the critical point. At the same time, the physical properties (e.g., composition, density and viscosity) of the two phases converge and they become chemically miscible. The two-phase flow then converges toward a literal single-phase flow, where the two phases have the same average velocity. The relative permeabilities approach (12) and (13) as a result. Furthermore, as n → w in this kind of miscible high-Ca limit, (14) reduces to the single-phase Darcy equation and the -exponent becomes meaningless. This paper is concerned with immiscible two-phase flow and we therefore do not seek to approach the high-Ca limit through small interfacial tensions.
Another way one might obtain the straight-line relative permeabilities, in the event that the high-Ca limit is approached while the fluids retain their immiscibility, is if the two fluids occupy similar, but separate and decoupled parts of the porous medium. In this case, the porous medium can essentially be described by the capillary tube model that gave = −1 above. The relative permeability of each fluid is then proportional to the cross-sectional area of the porous medium available to it, i.e., proportional to the saturation. Such a flow configuration seems compatible with one of the basic assumptions in the relative permeability framework, namely the fluids flow through connected pathways and the fluid-fluid interfaces behave as rigid partitions between them (Armstrong et al. 2016).
While the view that relative permeabilities are straight lines at high capillary numbers seems to be the prevailing one, there are also studies which indicate that this may not always be the case (Delshad 1981;Fulcher et al. 1985;Armstrong et al. 2016). A particularly illustrative example is the lattice-Boltzmann simulations by Armstrong et al. (2016), where the wetting and non-wetting relative permeabilities seem to approach curves that are concave up and concave down, respectively. High capillary numbers were simulated by varying interfacial tension, viscosity and driving force. Armstrong et al. (2016) show conclusively that ganglion motion and intermixing of fluids occur in a manner that is clearly inconsistent with the idea that the two fluids flow through static pathways. The latter point is part of a discussion that goes back to the flow regime micro-model studies of flow regimes by Avraam and Payatakes (1995).

Pore Network Simulations
In this section, we describe in some detail the pore network model used in this study, the numerical methods used to solve it and the procedure for calculating the time-averaged quantities described in the previous section. The pore network model is a dynamic model that takes both viscous and capillary forces into account. It keeps track of the fluid-fluid interface locations and evolves these in time according to the calculated pressures and flow rates in each pore. This provides a time series of the fluctuating total quantities Q w , Q n and Δp for the entire network, in steady state. Subsequently, these time series are averaged to get time-averaged steady-state flow properties.
Validation of the pore network modeling approach is provided through earlier works. These include Aker et al. (1998aAker et al. ( , 2000, Erpelding et al. (2013), Sinha et al. (2017) and the benchmark study by Zhao et al. (2019).

Pore Network Model
The model describes flow of two incompressible and immiscible fluids ( w and n ) in a porous medium. The porous medium is represented by a network consisting of N nodes that are connected by M links. The nodes are each given an index i ∈ [0, … , N − 1] , and the links are identified by the two nodes ij that they connect. An example pore network is shown in Fig. 2. The nodes have no volume, and the pore space volume is thus assigned to the links. Furthermore, it is assumed that each fluid fills the entire link cross sections. The location of a fluid-fluid interface can then be described by a single number which gives its position in the link. For each link, the vector ij contains the positions of the fluid interfaces in that link. The flow in each link is treated in a one-dimensional fashion, averaged over the link cross sections. We consider flows in relatively small cross sections only and therefore neglect any effects of fluid inertia. The volumetric flow rate from node j to node i through the link connecting the two nodes is then given by, Herein, p i is the pressure in node i, ij is the link's mobility and c ij is the net pressure difference across the link due to its fluid interfaces. Both ij and c ij depend on the interface positions ij . For two nodes i and j not connected by a link, ij = 0 . Applying mass conservation at each node i yields, The cross-sectional area of link ij is a ij . The interface positions ij therefore evolve in time according to the advection equation, when sufficiently far away from the nodes. Close to the nodes, they are subject to additional models that account for interface interactions in the nodes. This is further described in Gjennestad et al. (2018).

Link Mobility Model
The link mobility depends on link geometry and fluid viscosities. We assume cylindrical links when computing the mobilities and thus Here, L ij is the link length, r ij is the link radius and ij ij is the volume-weighted average of the fluid viscosities w and n .

Interfacial Pressure Discontinuity Model
There may be zero, one or more interfaces in each link. Their positions along the link are contained in ij . Each element in ij is thus between 0 and L ij . The symbol c ij denotes the sum of the interfacial pressure discontinuities in link ij. We assume that the links are much wider near the ends than in the middle and that the pressure discontinuities become negligibly small for interfaces near the ends. The pressure discontinuities are therefore modeled by Herein, wn is the interfacial tension and The effect of the -function is to introduce zones of length r ij at each end of the links where the pressure discontinuity of any interface is zero.

Numerical Solution Method
Solving the pore network model numerically involves stepping the fluid-fluid interface locations forward in time, from one discrete point in time to the next. This is be accomplished by applying a Runge-Kutta method to the M ODEs, one for each link, given by (17). Each evaluation of the right-hand sides during the integration requires simultaneously solving the mass conservation equations (16) and the constitutive equations (15) for the flow rate in each link. The time at the discrete point n is where the time step length Δt (i) is the difference between t (i+1) and t (i) . For ease of notation, we let quantities evaluated at the discrete time points be denoted by their associated time point index in superscript, e.g., Mobilities and pressure discontinuities with superscripts are evaluated with the interface positions at the time point indicated, In this work, we use the forward Euler method to integrate the ODEs (17). However, other and more sophisticated solution methods can also be used . Applying the forward Euler method to (17), we get The link flow rates herein are calculated by introducing the constitutive equation (15), evaluated at the current interface positions, into the mass conservation equations (16). This yields a system of linear equations, with one equation for each node i with an unknown pressure. One forward Euler time step is then taken by first obtaining the pressure in each node by solving the linear system, using the current interface positions. Subsequently, the link flow rates are calculated using (26) and the interface positions updated according to (25).

Time Step Restrictions
Since we use the explicit forward Euler method to integrate the interface positions in time, short enough time steps must be chosen to ensure numerical stability. The length of time step n is set according to the criteria in Gjennestad et al. (2018), where and the parameters C a and C c are set to 0.1 and 0.9, respectively.

Boundary Conditions
We carry out simulations in a network that can be laid out in two dimensions, as illustrated in Fig. 2b. The network is periodic both in the flow direction and in the transverse direction. Two different boundary conditions are explored: (1) A constant pressure difference of Δp is applied across the periodic boundary in the flow direction, and (2) a constant total flow rate Q is prescribed. The length of the network in the flow direction is denoted Δx and the average pressure gradient in the network is thus Δp∕Δx.
The boundary conditions corresponding to a constant applied Δp can be incorporated directly into the linear system (27) when this is solved during each time step. This is done by adding to or subtracting from each equation i in the linear system a term (n) ij Δp for each link ij that extends across the periodic boundary.
The boundary conditions corresponding to a prescribed flow rate Q are applied as follows. Since the model is linear, the total flow rate is related to the current pressure drop by where the constants C 1 and C 2 depend on the current interface positions in a non-trivial way. By selecting two different, but otherwise arbitrary, pressure drops Δp 1 and Δp 2 and solving the linear system (27) for both of them, we obtain the corresponding total flow rates Q 1 and Q 2 . The quantities C 1 and C 2 can then be calculated from, Subsequently, (31) can be solved for the pressure difference Δp necessary to get the prescribed total flow rate Q, and a forward Euler time step can be taken with this pressure difference incorporated into the linear system (27).

Computation of Steady-State Time Averages from Network Simulations
The porous medium we consider is modeled by a network of links, and the total volume of the links is the pore volume V p . The network is embedded in a three-dimensional block of solid material with thickness Δx in the flow direction and cross-sectional area A. The volume V of the porous block and its porosity are then easily calculated by (1) and (2), respectively.
The saturation S w may be computed at any time during the simulation, by adding up the fluid volumes for all links. However, since we use periodic boundary conditions, S w is a constant in each simulation. So is S n = 1 − S w .
In the case of constant applied pressure gradient Δp∕Δx , the quantities that we need to compute from the actual simulations are Q, Q w and Q n . These are time averages of the fluctuating quantities Q , Q w and Q n . The model is stepped forward in time as described in the previous section. We approximate the time-average Q by summing over the total flow rates Q (n) at each time step n (after steady state has been reached), The time-averaged quantities Q w and Q n are calculated from Q (n) w and Q (n) w in an analogous manner.
The instantaneous flow rate Q (n) can be computed by constructing a plane cutting through the network, transverse to the flow direction, and adding together the flow rates q (n) ij of all links intersecting the plane. We denote the set of intersecting links by B and add up, Since the fluids are incompressible, it does not matter where this cut is made. The instantaneous flow rate Q (n) w is computed by making several cuts, denote the set of cuts by C, and computing the sum Herein, |C| denotes the number of elements in C, i.e., the number of cuts, and s (n) ij is the volume fraction of wetting fluid in the volume of fluid that flowed past the middle of link ij during time step n. Q (n) n is computed in an analogous manner. Having computed the time averages Q, Q w and Q n we may the obtain the time-averaged flow velocity, mobility, fractional flow and relative permeabilities using (6), (7), (5), (8) and (9).
If Q is fixed instead of Δp∕Δx , the time-averaged value of the pressure gradient is computed by where Δp (n) is the pressure difference across the network during time step n.
Using steady-state averages calculated as described above, the capillary number is computed according to, where the average viscosity is defined as,

Dimensional Analysis
As can be surmised from the description above, the network and five numbers are given as input to steady-state simulations. In the case of constant pressure-difference boundary conditions, the five numbers are the fluid viscosities w and n , the fluid-fluid interfacial tension wn , the pressure gradient Δp∕Δx and the saturation S w . Any change in the steady-state averages is the response of the model to variations in these inputs. If we consider the network topology and aspect ratios fixed, and only allow for a linear scaling of the network size, any variations in the network can be described by a single length scale. We choose the average pore radius r.
By the Buckingham theorem (Rayleigh 1892), the total of six dimensional input variables can be reduced to three dimensionless variables. This means that any combination of the six inputs that give the same three dimensionless variables is similar and differs only in scale. Any dimensionless output from the model is therefore the same for the same values of the dimensionless input variables. One choice of dimensionless variables is where M is the viscosity ratio. The variable Π is a dimensionless pressure gradient. It represents the ratio of the average pressure drop over a length r to the Young-Laplace pressure difference over an interface in a pore of radius r . In particular, when Π = 1 , we have and the average pressure drop over the length r is equal to the typical Young-Laplace pressure difference.
Since it relates the average pressure drop to the capillary forces, Π may be expected to play a similar role as the capillary number. This should be true at least when capillary numbers are high and the average pressure drop is dominated by viscous contributions. However, Π is perhaps even more closely related to the ganglion mobilization number. This was defined by Avraam and Payatakes (1995) as the ratio between the driving force exerted on a ganglion and its resistance to motion resulting from capillary forces.

Simulations
Steady-state simulations were performed using the pore network model described in this section. All simulations were run on 72 × 48 hexagonal networks, similar to that shown in Fig. 2b. These networks consisted of 3456 nodes and 5184 links. All links had the same length L, and link radii were uniformly distributed between 0.1L and 0.4L. For each of the 288 combinations of the input parameters listed in Table 1, 21 values of S w were used, evenly spaced on the interval [0, 1] . In total, 288 × 21 = 6048 simulations were run. As this paper is concerned with immiscible two-phase flow, we do not use vanishingly small values of the interfacial tension. Time-averaged quantities were calculated from simulation results as described in Sect. 3.3. The averaging time corresponded to 10 pore volumes of flow.

Relative Permeabilities
Relative permeabilities r w and r n are perhaps the most extensively studied properties in two-phase flow in porous media, and the most obvious dimensionless steady-state quantities to calculate from the pore network model. In Fig. 3, computed relative permeabilities are plotted against saturation S w and the non-dimensional pressure gradient Π for three subsets of the simulations. Specifically, Fig. 3a, b, c shows relative permeabilities for the wetting phase and viscosity ratios M of 1, 0.25 and 4, respectively. Relative permeabilities for the non-wetting phase and a viscosity ratio of 4 are shown in Fig. 3d.
For each value of M, the relative permeability data fall on a single well-defined surface. This shows that the relative permeabilities produced by the model are indeed determined by the three dimensionless variables S w , M and Π , in agreement with the dimensional analysis in Sect. 3.4. Bardon and Longeron (1980) mention that gravity (Bond number), wettability (contact angle) and inertia (Reynolds number) could also affect the relative permeabilities. These effects are not considered in the present study, though gravity could be included in the model with relative ease. When measuring relative permeabilities (Oak et al. 1990;Bennion and Bachu 2005) and when using relative permeability models for continuum-scale calculations, it is often only their dependence on S w which is considered. Here, the calculated relative permeabilities are strongly dependent on Π , and they increase with increasing Π . There is a very strong correlation between Π and the capillary number, where high values of Π are also associated with high values of Ca . Thus, the results are consistent with those of Bardon and Longeron (1980), Avraam and Payatakes (1995), Armstrong et al. (2016) and Datta et al. (2014), who find that relative permeabilities increase with capillary number. This dependence seems to disappear, however, as Π → ∞ . For the viscosity ratios considered here, it disappears at Π ∼ 1 , where the contours in Fig. 3 approach vertical lines. At this Π-value, the average pressure drop over the length r is equal to the typical Young-Laplace interfacial pressure difference, as discussed in Sect. 3.4. The fact that the relative permeabilities become independent of the pressure gradient as capillary numbers increase is consistent with the high-Ca limit.
The general consensus in the literature (Bardon and Longeron 1980;Fulcher et al. 1985;Avraam and Payatakes 1995;Whitson et al. 2003;Ramstad et al. 2012;Schechter and Haynes 1992) seems to be that relative permeabilities approach straight lines (12) and (13), at high capillary numbers. In the equal-viscosity pore network simulations by , this was found to be the case. Here, however, we find straight lines, i.e., equidistant contour lines in Fig. 3 for Π ≳ 1 , only for M ∼ 1 . When M is different from unity, relative permeabilities converge to nonlinear functions of S w (and M) in the high-Ca limit.
As discussed in Sect. 2.1, relative permeabilities form straight lines in the high-Ca limit if this is realized by approaching a mixture critical point. In such a case Ca → ∞ would necessarily happen at the same time as M → 1 and our results would seemingly be in agreement with the literature consensus. However, our model assumes immiscible flow and it is questionable whether it can be relied on when the fluids become miscible.
Straight-line relative permeabilities could be obtained also in the immiscible high-Ca limit if the fluids flow in similar, decoupled flow channels. By visual inspection of fluid configurations, however, such decoupled flow channels are not found here (see, e.g., Fig. 10 in "Appendix" section). Instead, the fluids exhibit a large degree of intermixing and move as droplets or small ganglia at high capillary numbers. This was observed also by Sinha et al. (2019b), both in pore network model and lattice-Boltzmann simulations. Disconnected non-wetting droplets were also observed at high capillary numbers in the experiments by Avraam and Payatakes (1995), and were found to contribute significantly to the total flow rate, although connected pathways were also present.
We therefore conclude that our relative permeabilities deviate from straight lines at high capillary numbers because the fluids are not in the highly miscible near-critical region. Instead they have a viscosity disparity and intermix instead of forming decoupled flow channels. The effect of this on total mobility and fractional flow is discussed further in Sects. 4.3 and 4.4, respectively.
One example from the literature where relative permeabilities do not seem to form straight lines is found in Armstrong et al. (2016). Instead, they found a concave down r n and a concave up r w , similar to our results with M > 1 (see Fig. 3c, d). Their viscosity ratio was not given. From Fig. 3, it is evident that the relative permeabilities follow a nonlinear curve not unlike those produced by the classical Corey-type correlations for the lowest values of Π , i.e., the lowest capillary numbers. When working with such correlations, it is typically assumed that there exists a low capillary number limit below which relative permeabilities become independent of the flow rate, and the correlations are therefore valid (for the fluids used in the measurements). Ramstad et al. (2012) mention that viscous forces start to influence the fluid transport at capillary numbers around 10 −5 . This is consistent with the findings here, which are that relative permeabilities have a dependence on Π down to the lowest capillary numbers considered of approximately 10 −4 . We emphasize that the definition of capillary number used here differs from that used in Ramstad et al. (2012), since it includes the porosity. Adoption of the definition from Ramstad et al. (2012) would reduce all capillary numbers reported here by approximately half an order of magnitude.
Qualitatively, the changes in the relative permeability curves with Ca in Armstrong et al. (2016) are similar to those found here for M > 1 . In particular, their r n -curve shifts from a concave up shape at low Ca to a concave down shape at high Ca . We find the same here when M > 1 (Fig. 3d). Avraam and Payatakes (1995) find from their experiments that both relative permeabilities increase with M. This is not the case here, at least not at high capillary numbers. They attribute this to the existence of films of the wetting fluid, which are not included in our model.

Residual Saturations
For both the wetting and non-wetting fluids, there are regions in Fig. 3 where the relative permeabilities are zero. This was seen by  also, for M = 1 . These regions correspond to irreducible/residual saturations, two other ubiquitous dimensionless quantities in two-phase porous media flow. The residual saturations are often defined as the saturation of one fluid that remains after flooding with the other. This property, defined in this way, is somewhat difficult to measure using the type of steady-state pore network model simulations performed here. Therefore, we have chosen to define the residual saturation of the wetting fluid as the saturation on the fractional flow curve ( F w plotted vs. S w , all other input parameters kept constant) where the wetting fluid fractional flow falls below Computed residual wetting and non-wetting saturations are shown in Fig. 4. Residual saturations increase as capillary numbers are reduced, in accordance with findings of Ramstad et al. (2012), Bardon and Longeron (1980), Datta et al. (2014) and Fulcher et al. (1985). Furthermore, they reach zero at capillary numbers around 0.1. This means that it is possible to flush out all of one fluid from the network model through flooding with the other, provided that the flow rate is high enough. Bardon and Longeron (1980) observed that residual saturations were insensitive to changes in M, which is in agreement with our findings. However, wetting residual saturations are somewhat higher when the wetting fluid is more viscous and non-wetting residual saturations are a little higher when the non-wetting fluid is more viscous. It would seem that mobility of the minority fluid is impeded when it is more viscous than the majority fluid.

Average Flow Velocity and Mobility
The average mobility m and the average flow velocity v are other important quantities in two-phase flow. They are related through (7) and are discussed together here for reasons that will become apparent below. Sinha et al. (2019b) studied the high-Ca limit of two-phase porous media flow. They found that, at high capillary numbers, the average flow velocity followed a Darcy-type equation (10), with an effective viscosity (11). The existence of this high-Ca limit motivates the study of the average flow velocity and the average mobility, relative to their limit values. Dividing (7) by (10) gives where m D is the mobility in the high-Ca limit. The two quantities v∕v D and m∕m D are thus identical. Moreover, they are dimensionless and may be expected to vary, roughly, between 0 and 1. In particular, they should be 1 in the two single-phase cases, S w = 1 and S n = 1 , and in the high-Ca limit. Figure 5 shows v∕v D for M = 1 , M = 0.25 and M = 4 , plotted against S w and Π . As expected, all data points collapse to 1 in both single-phase cases. Furthermore, each value of M corresponds to a single well-defined v∕v D -surface, in accordance with the dimensional analysis. Each constant-M surface reaches values close to 1 at the highest values of Π , in agreement with the findings of Sinha et al. (2019b) for the high-Ca limit.
Interestingly, there are some values of v∕v D that are larger than 1. This means that, at a given saturation, the average mobility is higher at some lower capillary number than in the high-Ca limit. This somewhat counter-intuitive effect occurs for the more disparate viscosity ratios, at saturations where the more viscous fluid is in minority. Figure 6a shows v∕v D plotted against Π , for S w = 0.15 and three different viscosity ratios, 0.25, 1 and 4. The data points with M = 1 converge to 1, the limit value, from below and relatively fast as Π increases. The data points with M = 4 also approach 1 from below, but slower. For the lower viscosity ratio M = 0.25 , on the other hand, v∕v D increases fast, overshoots and then approaches 1 from above.
In Fig. 6b is shown the fractional flow for the same data points as in Fig. 6a. For the data points with M = 1 , convergence of v∕v D to the limit value occurs as the fractional In terms of mobility ratios m∕m D these observations may be understood as follows. At low pressure gradients, all wetting fluid is stuck, in the sense that F w = 0 , and the non-wetting fluid flows around it (see Fig. 6b). As the pressure gradient is increased, some of the wetting fluid is mobilized and F w increases above zero. This results in more active flow paths for both fluids and a sharp increase in the average mobility for all three viscosity ratios. An example of such mobilization can be seen by comparing Figs. 8 and 9 in "Appendix" section.
For M = 0.25 , average mobility reaches a maximum before all wetting fluid is mobilized, i.e., before F w converges to its value in the high-Ca limit. This maximum is caused by the competition between two different effects. First, an increase in pressure gradient makes more flow paths available, increasing mobility. Second, F w increases and the more viscous wetting fluid makes up a larger fraction of the flowing fluid. Thus, the average viscosity of the flowing fluid increases, reducing the average mobility. Eventually, a point is reached where the latter effect becomes more important and a further increase in the pressure gradient reduces the average mobility. For M = 1 , there is no such competition to generate a maximum, as the wetting and non-wetting fluids are equally viscous and mobilization of the wetting fluid does not affect the average viscosity.
For M = 4 , the two effects are again present. However, since the wetting fluid is now less viscous, they both lead to an increase in mobility with an increase in pressure gradient and we see no maximum.

Fractional Flow
The fractional flow for three subsets of the performed simulations is shown in Fig. 7, for viscosity ratios 1, 4, and 0.25. The data in Fig. 7a for M = 1 are in qualitative agreement with those from . They find F w ∼ S w at high capillary numbers and that deviation from the diagonal line representing F w = S w increases as the capillary number is reduced. Furthermore, curves for a specific capillary number are asymmetric w.r.t. S w = 0.5 and cross the diagonal line at S w > 0.5 , meaning that more of the curve lies below the diagonal than above it. This observation was explained by  by the propensity for the wetting fluid to occupy narrower pores where the flow rate is lower.
By comparing Fig. 7b, c we may deduce some of the impact of the viscosity ratio on the fractional flow. At high capillary numbers, F w > S w when M > 1 , i.e., when the wetting fluid is less viscous. Conversely, F w < S w when M < 1 and the wetting fluid is more viscous. The latter was observed also by Avraam and Payatakes (1995), at low viscosity ratios and high capillary number, the fractional flow curves tended to be concave up.
The dotted lines in Fig. 7b, c represent the fractional flows obtained if the relative permeabilities were r w = S w and r n = 1 − S w , i.e., what would be expected at high capillary numbers if the fluids followed separate, similar flow channels. The fractional flows from the network model in this limit curve the same way as the dotted lines, but lie much closer to the diagonal. As noted before, our pore network model gives instead a large degree of fluid intermixing in the high-Ca regime (see, e.g., Fig. 10). In such an intermixed state, fluid velocities will be more tightly coupled. We therefore propose that, in the event that fluids intermix rather than form decoupled flow channels, the tighter coupling causes the fluids to flow with velocities that are more similar. Thus, the fractional flow curves lie closer to the diagonal.
At lower capillary numbers, the fractional flow curves obtain the classical S-shape, as in the case for M = 1 . Also, as is intuitive and was observed by Avraam and Payatakes (1995), fractional flow for a given saturation and capillary number increases with viscosity ratio.

Conclusion
We have performed more than 6000 steady-state simulations with a dynamic pore network model of the Aker type, corresponding to a large span in viscosity ratios and capillary numbers. From these simulations, dimensionless time-averaged steady-state quantities such as relative permeabilities, residual saturations, mobility ratios and fractional flows were computed and discussed. By a dimensional analysis of the model, all dimensionless output was found to be functions of the saturation S w , the viscosity ratio M and the dimensionless pressure gradient Π . Effects of wettability, gravity and inertia were not considered. These effects may add additional dimensionless variables whose impact could be studied in future work.
Calculated relative permeabilities and residual saturations showed many of the same qualitative features observed in other experimental and modeling studies. In particular, the relative permeabilities increased with capillary numbers and converged to a limit, dependent on M and S w , at high capillary numbers. However, while consensus in the literature seems to be that relative permeabilities converge to straight lines at high capillary numbers, this was not the case for the results of the network model when M ≠ 1.
Our conclusion is that our relative permeabilities deviate from straight lines at high capillary numbers because the fluids are not in the highly miscible near-critical region. Instead they have a viscosity disparity and intermix rather than forming decoupled, similar flow channels. Such intermixing behavior has been observed previously in pore network and lattice-Boltzmann simulations (Sinha et al. 2019b) and, to some extent, in experiments (Avraam and Payatakes 1995). However, it would be very interesting to see whether experimental studies specifically designed to induce intermixing and measure steady-state properties at high capillary numbers would produce relative permeability curves that are nonlinear in S w .
Ratios of average mobility to their high capillary number limit values were also considered. These ratios varied, roughly, between 0 and 1, but values larger than 1 were also observed. For a given saturation, the mobilities were not always monotonically increasing with the pressure gradient. While increasing the pressure gradient mobilizes more fluid and activates more flow paths, when the mobilized fluid is more viscous, a reduction in average mobility may occur instead.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix: Example Fluid Distributions
In this section, we show example fluid distributions from steady-state simulations at different sets of parameters. We show two examples for each set of parameters, with a time difference between them corresponding to one pore volume of flow. For ease of viewing and discussion, the networks used in these illustrations are 12 × 8 , much smaller than those used in the steady-state simulations described in Sect. 3.5. Figures 8, 9 and 10 show example fluid configurations obtained during simulations with S w = 0.85 and M = 1 and different dimensionless pressure gradients. Comparing them illustrates how more and more of the non-wetting fluid is mobilized as the magnitude of Π is increased.
In Fig. 8, log 10 (Π) = −1.0 and all interfaces remain stationary. Only the wetting fluid is moving ( F w = 1.0 ) and necessarily does so through connected pathways.

3
In Fig. 9, log 10 (Π) = −0.7 and some interfaces remain stationary while others are moving. The non-wetting fluid is partially mobilized ( F w = 0.94 ) and moves in droplets or small ganglia. In Fig. 10, log 10 (Π) = 1.2 and the system starts to approach the high-Ca limit. Almost no interfaces remain stationary, and most of the non-wetting fluid is moving through droplets or ganglia. The fractional flow is F w = 0.88.  Figure 11 shows example fluid configurations when more non-wetting fluid is introduced, S w = 0.75 , while M = 1 and the pressure gradient log 10 (Π) = −1.0 is the same as in Fig. 8. As more non-wetting fluid is added, some of it is mobilized, while some remains stationary ( F w = 0.84 ). The non-wetting fluid moves as droplets or small ganglia.