Effective rheology of immiscible two-phase flow in porous media consisting of random mixtures of grains having two types of wetting properties

We consider the effective rheology of immiscible two-phase flow in porous media with random mixtures of two types of grains with different wetting properties using a dynamic pore network model under steady-state. Two immiscible fluids A and B flow through the pores between these two types of grains denoted"+"and"-". Fluid A is fully wetting and B is fully non-wetting with respect to"+"grains and opposite with"-"grains. The direction of the capillary forces in the links between two"+"grains is therefore opposite compared to that between two"-"grains, whereas the capillary forces in the links between two opposite types of grains average to zero. For a window of grain occupation probabilities, a percolating regime appears where there is a high probability of having connected paths with zero capillary forces. Due to these paths, no minimum threshold pressure is required to start a flow in this regime. While varying the pressure drop across the porous medium from low to high in this regime, the relation between the volumetric flow rate and the pressure drop goes from being linear to a power law with exponent 2.56 to linear again. Outside the percolation regime, there is a threshold pressure. No linear regime is observed for low pressure drops. When the pressure drop is high enough for there to be flow, we find that the flow rate depends on the excess pressure drop to an exponents around 2.2-2.3. At even higher pressure drops, the relation is linear. We see no change in exponent for the intermediate regime at the percolation critical points where the zero-capillary force paths disappear. We measure the mobility at the percolation threshold at low pressure drops so that the flow rate versus pressure drop is linear. Assuming a power law, the mobility is proportional to the difference between the occupation probability and the critical occupation probability to a power of around 5.7.


I. INTRODUCTION
It was in 1827 that Ohm published his law stating that electrical current is proportional to the voltage drop across a conductor [1], meeting fierce resistance from the physics community in the beginning. Darcy arrived in 1856 at a similar law for single-phase flow in porous media, i.e., the volumetric flow rate is proportional to the pressure drop across the porous medium [2]. Both of these fundamental laws are examples of there being a linear relationship between current and driving force. In the case of the Darcy law, the derivation based on pore scale physics has been a challenge, see e.g., Whitaker's derivation based on momentum transfer [3].
The Darcy law for single-phase flow through a porous sample is where Q is the volumetric flow rate along the axis of the cylindrical sample, ∆P the pressure drop along it in the flow direction, A the area of the sample orthogonal to the flow direction, K the permeability of the sample, µ the viscosity of the liquid and L is the system length.
In 1936 the Darcy law (1) was generalized to the simultaneous flow of two immiscible liquids by Wyckoff and Botset by essentially splitting it into two [4], where the subscripts w and n refer to the wetting properties of the two fluids with respect to the matrix; w refers to the more wetting fluid and n to the less wetting fluid. The idea behind this split is simple. The wetting fluid will see a pore space reduced by the presence of the other fluid, leading to a reduction in effective permeability for the wetting fluid. The reduction parameter is the wetting relative permeability k rw . Completely analogously the non-wetting fluid sees an effective reduction of the permeability by a factor k rn , the non-wetting relative permeability. The split was given physical contents when Wyckoff and Botset assumed that the two relative permeabilities were functions of the wetting saturation S w alone; the wetting saturation being the pore volume occupied by the wetting fluid divided by the total pore volume and assuming the fluids are incompressible. Barenblatt et al. [5] have later shown that this assumption is valid if there exists a local phase equilibrium between the fluids, a condition that is fulfilled only for slow flows. A further assumption built into equations (2) and (3) is that there are no macroscopic saturation gradients present. The total volumetric flow rate is given by the sum of the volumetric flow rates of each fluid, and as a consequence, the generalized Darcy equations (2) and (3) predict that is, a total volumetric flow rate being proportional to the pressure drop. Equations (2) and (3) assume that there are no macroscopic saturation gradients. If this is not the case, the pressure is split into one associated with the non-wetting fluid, P n , and one associated with the wetting fluid, P w . Their difference is equal to the capillary pressure function, P n − P w = P c (S w ), which is also assumed only to depend on the saturation S w . Equations (2) and (3) will then contain terms of the type ∇P c = (dP c /dS w )∇S w , thus setting up the pressure gradient and the saturation gradient as driving forces. When these equations are combined with mass conservation, the result is a closed set of equations that determine how the saturation develops within the porous medium.
When the saturation changes inhomogeneously in the porous medium with time, one implicitly assumes that fluid interfaces move within the porous medium. It was then a surprise when Tallakstad et al. [6,7] reported a flow rate Q depending on ∆P as with β ≈ 1.85 for a two-dimensional glass-bead-filled Hele-Shaw cell filled with a water-glycerol mixture and air in the flow regime where the generalized Darcy equations (2) and (3) are supposed to be valid. This study was followed up by an NMR study of the three-dimensional glass bead packings by Rassi et al. [8] finding an exponent β varying between 2.2 and 3.3. Aursjø et al. [9] using the same model porous medium as Tallakstad et al. [6,7], but with two incompressible fluids, found β ≈ 1.5 or 1.35 depending on the fractional flow rates. Similar results in the sense that β is considerably larger than one, have since been observed by a number of groups, see [10][11][12][13]. There has also been a considerable effort to understand these results theoretically and reproduce them numerically [6,7,[14][15][16][17][18][19][20][21][22][23][24][25].
It should be pointed out that the power law behavior seen in Eq. (6) is different from the one described by Wilkinson in 1986 [26]. In that work, Wilkinson used the invasion percolation model to work out the dependence of the relative permeabilities on the capillary pressure, which could be linked to the saturation. He would find that the non-wetting relative permeability k rn would depend on the difference between the capillary pressure P c and a critical capillary pressure P c c related to the percolation critical point as where t is the percolation conduction exponent [27]. This is, however, a very different problem from the one giving rise to Eq. (6). The power law in (7) is a direct reflection of the geometry of the clusters of the non-wetting fluid in the system after the invasion process. Hence, it is a static problem. The power law in (6) is, as we shall see, the result of a dynamic process caused by the motion of the fluid interfaces. The power law behavior in equation (6) is due to a competition between the capillary and the viscous forces. It is straightforward to understand why the flow rate should increase faster than linear when these forces are in competition. When the pressure difference across the porous medium is increased, more interfaces are beginning to move leading to a higher effective permeability [28]. Why it should be a power law is less obvious. The best argument for why was perhaps already given by Tallakstad et al. [6,7] through comparing pressure drop across fluid clusters with the capillary pressures holding them in place. Capillary fiber bundle models [29,30] are porous media in the form of bundles of capillary fibers and they are typically simple enough to be mathematically solvable [16, 19-21, 24, 25]. When the fibers have undulating radii along the long axis, they show non-linear volumetric flow rate vs. pressure drop, but not quite of the form (6), but rather where P t is a threshold pressure necessary for flow to occur, P max is the maximum threshold pressure found in any capillary fiber, and M and M D are mobilities. A nonzero threshold pressure is in general necessary in porous media when neither of the two immiscible fluids percolates when dealing with porous media and not just the capillary fiber bundle model [10,15,21]. The existence of a non-zero threshold pressure makes the measurement of β much harder than when it is zero as this implies determining two parameters simultaneously, (P t , β), rather than just one, β.
A central unanswered question is whether the exponent β is universal in the sense that there are classes of systems that all have the same value, i.e., can one define universality classes? Intuitively this is a very appealing idea as one has a diverging length scale as in equilibrium critical phenomena as |∆P | → P t from above [6,7]. The experimental measurements of β have so far not given any indication of the existence of universality classes and neither have the computational efforts due to the difficulties in dealing with two unknown parameters, P t and β.
Roy et al. [19] found using a capillary fiber bundle model that β = 2 if the fiber-to-fiber probability distribution of thresholds includes P t = 0 with a finite probability, otherwise β = 3/2. The fibers here had smoothly undulating radii along the flow direction. Lanza et al. [24] who studied non-Newtonian a mixture of immiscible Newtonian and non-Newtonian fluids in a capillary fiber bundle model found a different value of β when the radius distribution is jagged from when it is smooth.
Equation (8) which was derived for the capillary fiber bundle model, predicts there being a pressure drop |∆P | = P t below which the flow rate Q is zero. This threshold may be zero. Within the capillary fiber bundle model, this means that some capillary tubes belonging to the bundle have interfaces that move as soon as there is a pressure difference across them. There is, however, one important mechanism missing in the capillary fiber bundle models: In the porous medium, the immiscible fluids may be percolating. That is, there are pathways through the porous medium along which there are no interfaces. In this case, there will be a linear regime when the pressure drop is low enough so that the interfaces surrounding the percolating paths do not move. When the pressure drop is increased sufficiently for them to do so, the non-linear power law regime sets in.
Recently Fyhn et al. [21] studied the exponent β and the threshold pressure P t in a capillary fiber bundle model and a dynamic pore network model under mixed wet conditions. In the dynamic pore network model, each link was given a wetting angle -in the sense that if there is an interface in the link, this is the angle it will make with the walls of the tube -drawn from a given probability distribution. In the capillary fiber bundle model, each undulating tube is given a wetting angle from a given probability distribution. In both models, a constitutive law of the form (8) was found. The capillary fiber bundle model could be solved analytically, yielding The network model studies showed a less clear picture, with β varying between 1 and 1.8 depending on the saturation and the wetting angle distribution. It was not possible to resolve whether there were regions of fixed β or whether it varied continuously with the parameters of the model. This was due to the non-zero threshold pressure P t which needed to be determined together with β.
We study here a model for immiscible two-phase flow in a porous medium made from two types of grains that have different wetting properties with respect to the fluids. The model treats the interfacial tension between the two fluids similarly to a model introduced by Irannezhad et al. [31,32]. We imagine a packing of two types of grains, say type "+" and type "−". Two immiscible fluids, denoted by "A" and "B" flow through the pores between the grains denoted by "+" or "−". Fluid "A" is fully wetting and "B" is fully non-wetting with respect to "+" grains whereas it is the opposite with respect to "−" grains. The direction of the capillary forces in the links between two "+" grains is therefore opposite compared to the direction in the links between two "−" grains, whereas the capillary forces in the links between two opposite types of grains are zero.
The probability that a grain is of "+" type, is p + . A second parameter is the wetting saturation S w . There is a rich phase diagram when plotting the threshold pressure P t as a function of the two control variables p + and S w which is illustrated in Fig. 1. Note in particular in this phase diagram that there is a region in the middle where the threshold pressure P t = 0. This region is limited by two p + = constant critical lines. Each line signifies a percolation transition [27]. The two curved grey lines signify a possible shift of the two blue transition lines due to the dynamics of the model. There are also two other lines, one green line marked "hysteretic transition" and one red line marked "non-hysteretic" transition. Crossing such a line, one of the two fluids stops moving and we are essentially dealing with a single-phase flow problem. When the wetting fluid stops moving, there is no hysteresis. On the other hand, when the non-wetting fluid stops, there is hysteresis in the sense that the wetting saturation has to be lowered more to get the non-wetting fluids moving again [33].
In the region where P t = 0, while still having twophase flow, we observe an exponent β = β 3 = 2.56 ± 0.05 for saturation S w = 0.5. This value is also seen when setting p + = p c or p + = 1 − p c , where p c ≈ 0.5927 is the site percolation threshold for square lattices, i.e., β = β 2 = β 3 . For p + much lower than 1 − p c or p + much higher than p c , still for saturation S w = 0.5, we see β = β 1 = 2.25 ± 0.1. The large uncertainty in β seen here stems from P t > 0.
It is surprising that β 2 = β 3 within the precision we are able to obtain. The exponent β 2 is obtained at the percolation threshold where the paths where fluid surfaces meet no resistance are fractal with fractal dimension 4/3 as they are the external perimeters of percolation clusters [34]. The reason for not seeing critical behavior reflected in β comes from there also being other links that have no interfacial tension in them as they contain no interfaces, thus driving the system away from criticality. In order to investigate whether there are any traces at all in the transport properties of the percolation critical point, we have studied the mobility M at low pressure drops as p + approaches a critical value, where we expect it to vanish with an exponent t ′ of the same type as the conductivity exponent t percolation in ordinary percolation in the vicinity of the critical p + [27,35]. We find that t ′ ≈ 5.7, indicating that the system is not critical and M falls off faster than algebraic. We therefore expect there to be two extra transition lines (marked in grey in Fig. 1) that distinguish between P t = 0 and P t > 0. The nature of these lines is unknown.
We use a dynamic pore network model [36][37][38][39][40] for this study. It has been used earlier in the context of modeling mixed wetting porous media, see [21,41,42]. We describe the model in Section II including our use of the wetting model similar to the one introduced by Irannezhad et al. [31,32]. Section III explains how we identify the paths through the network that have no capillary forces associated with them and relate them to a site percolation problem. Section IV presents the analysis of the low pressure drop mobility at the percolation critical points. Section V constitutes our investigation of volumetric flow rate Q vs. pressure drop ∆P . We fix the saturation S w = 0.5 and scan through this line in the phase diagram in Fig. 1 for different values of p + . We also tested whether there would be hysteresis with respect to increasing or decreasing the pressure drop, finding none. Section VI contains a summary and our conclusions.

II. DYNAMIC PORE NETWORK MODEL
A sketch of the dynamic pore network (DPN) model used in this work is given in Fig. 2, showing a square two-dimensional network with links with the same length tilted 45 • angle from the flow direction. The ∆P across the network drives the flow leading to a Q which is measured over a cross-section of the system normal to the direction of the overall flow. The zoomed-in sketch to the right in Fig. 2 illustrates the rules for using the wetting properties of the grains to assign wetting angles θ to the links, where θ is consistently defined through one of the fluids. In contrast to earlier models [21,41,42] that assign the wetting angles to the pores or links directly, the physical basis for this model is a mixture of grains and the wettability of the pore space in-between depends on the wettability of the surrounding grains, similar to the system introduced by Irannezhad et al. [31,32]. We assume two types of grains, they being either fully nonwetting with θ = 180 • or fully wetting with θ = 0 • . Having fully non-wetting or fully-wetting grains maximizes the difference between the two types of grains in terms of their wettability and hence maximizes any impact on the rheology that comes as a result of this difference. The grains are denoted fully non-wetting and assigned a notation "+" with an occupation probability p + and the rest of the grains are then fully wetting with a notation "−". For each link, θ is determined based on the link's adjacent grains. Each grain in the network is connected to four links which means each link has two adjacent grains, as shown in Fig. 2. If both of the adjacent grains are assigned "+", the link will have θ = 180 • . If both of the adjacent grains are assigned "−", then θ = 0 • . Lastly, if one of the adjacent grains is "+" and the other one is "−", the link in the middle should be easy to pass through for both fluids and the wettability should be neutral with θ = 90 • . The dynamic pore network model implemented on a square lattice consists of links oriented 45 • from the overall flow direction. The flow is driven by a global applied pressure ∆P and the total volumetric flow rate Q is measured over a cross-section normal to the direction of the overall flow. The wetting angle θ of each link is based on its adjacent grains. The grains are assigned "+" with an occupation probability p + and the rest of the grains are assigned "−". If both of the adjacent grains are assigned "+", θ = 180 • (marked pink). If both of the adjacent grains are assigned "−", θ = 0 • (marked blue). Lastly, if one of the adjacent grains is "+" and the other one is "−", then θ = 90 • (marked black) and hence there are no capillary forces associated with interfaces in the link. along the link's center-axis x, is given by where it has been assumed that the radius does not deviate too much from its average valuer [36]. Here, µ = s A µ A + s B µ B is the saturation weighted viscosity of the fluids where s A and s B are saturations of the two fluids A and B respectively with viscosities µ A and µ B in the links (in contrast to S w which is the average saturation over the whole network). Figure 3 can be used to further explain the variables in Eq. (10). The unit vector J lies along the x axis, and points in the direction out of the fluid within which θ is defined. In Fig. 3, θ is consistently measured through fluid A in both examples (a) and (b) regardless of if fluid A is more or less wetting with respect to the solid. Hence,Ĵ also consistently points across the interface starting from fluid A towards fluid B. The sum in Eq. (10) is taken over the interfaces numbered k with varyingĴ k with positions x k ∈ [0, l] along x. The dot product of this sum with the unit vectorx in the positive x direction is taken afterward to obtain the total capillary pressure. The capillary pressure across one interface at position x which has an angle θ with the solid through fluid A is modeled by using the Young-Laplace equation [43], where σ is the surface tension and is the radius where a is the amplitude of the periodic variation and r 0 /a is randomly chosen from the interval [0.1l, 0.4l]. This way, p t varies with both the position along a link and from link to link.
For all simulations in the following, the two immiscible fluids have been given surface tension 3.0 · 10 −5 N/mm and viscosity 0.1 Pa·s for both. The overall network saturation is kept constant at 0.5, meaning there is equal amounts of the two fluids. The links in the network have length l = 1 mm. In all the figures, the logarithms are in base 10.

III. EASY LINKS AND CONNECTED PATHS
There are three types of links in the model: those that are of the "+ +" type, those that are of the "− −" type, and the "+ −"="− +" type. We will in the following refer to the latter type as "easy links" since they offer no capillary resistance to interfaces that happen to be in them. Paths of connected easy links may percolate, i.e., stretch across the network forming loops as we are implementing bi-periodic boundary conditions. We will refer to such percolating paths of easy links as "connected paths", see Fig. 4. The geometry of the easy links and connected paths may be mapped onto an ordinary site percolation problem [27]. The links altogether form a square lattice. The nodes of the dual lattice, form another square lattice [44] and are assigned "+" or "−". These values are placed at random. The distribution of neighboring "+" sites in this dual lattice form an ordinary site percolation problem. In an infinitely large lattice, there will be a percolating "+" cluster when p + ≥ p c , where p c is the site percolation threshold 0.5927 . . . . If we, on the other hand, focus on the "−" sites, there will be a cluster of such sites that percolate if p − = 1−p + ≥ p c , or p + ≤ 1−p c ≈ 0.4073 . . . [45]. Hence, if 0 ≤ p + ≤ 1 − p c , the "−" clusters percolate, if 1 − p c ≤ p + ≤ p c , neither the "−" sites nor the "+" sites percolate, and if p c ≤ p + , the "+" sites percolate. We show in Fig. 5 a map of the wetting angles associated with different values of p + . The easy links are shown in black.  Fig. 2. This means that the pore space between these grains, namely the links in DPN, all have θ = 0 • . Oppositely at p + = 1.0, shown in (e), there are only links that have θ = 180 • . In these two extreme cases, there is no easy link in the network with neutral-wettability θ = 90 • . Moving away from these extremes, when p + = 0.3 in (b) or when p + = 0.7 in (d), links with θ = 90 • are present but not enough to create a connected path that crosses the entire system. At the middle point of p + = 0.5,(c), DPN has half of each type of grain creating the highest possible probability for having connected paths with only θ = 90 • links. In these examples, p + = 0.5 (c) is the only one that lies within the limit 1 − pc < p + < pc, and it is only here we find connected paths.
We note that if neither the "+" sites nor the "−" sites percolate (1 − p c ≤ p + ≤ p c ), there must be connected paths. We furthermore note that if either of the two site types percolates, there cannot be any connected paths. At the two thresholds, p + = 1 − p c and p + = p c , the connected paths appear together with the appearance of a percolating cluster of either "−" or "+" type, as the perimeter of the incipient percolating cluster is a connected path. At the percolation thresholds, we know that the fractal dimension of the perimeter, and hence the corresponding connected path, is 4/3 [34]. For values away from the critical points, the connected paths are not fractal. Hence, the structure of the easy link clusters and the connected path is very different away from the critical points, while still being in the interval 1 − p c ≤ p + ≤ p c .
The probability of finding a connected path as a function of p + is investigated by testing 1000 randomly generated networks with size L × L for each p + ∈ {0.3000, 0.3001, 0.3002, . . . , 0.7000}. The results are shown in Fig. 6 for L = 50 links and L = 100 links. We see that the two curves cross very close to 1 − p c and p c .

IV. MOBILITY
As we will show with the results presented in the next section, the constitutive law between the volumetric flow rate Q and the pressure drop |∆P | can be written as in the region 1 − p c ≤ p + ≤ p c . Here, P l and P u are two crossover pressures. There are three regimes: (1) a linear regime at low pressure drop, (2) a non-linear regime for intermediate pressure drops and (3) a linear regime for high pressure drops. Each regime is characterized by a mobility, M (p + , S w ), M m (p + , S w ), and M D (p + , S w ), respectively.
If we move to values of p + where P t > 0, regime (1) disappears. Hence, we have that M (p + , S w ) tends to zero as p + reaches the boundary between the P t = 0 region and the P t > 0 region. We hypothesize in the following that the boundaries of this region are given by the percolation thresholds 1 − p c and p c . Expecting that M (p + , S w ) shows similar behavior to the conductivity in percolation [35], we make the assumption that the mobility vanishes as where t ′ is a transport exponent of the same type as the conductivity exponent t in ordinary percolation, which is 1.303(8) according to [46]. In Eq. (14), p + → (1 − p c ) + means p + approaches 1 − p c from above and p + → (p c ) − means p + approaches p c from below. By using finite size scaling analysis, we have that where ν is the correlation length exponent in percolation, which is known to be 4/3 [47].
To investigate the relation given in Eq. (15), we set p + = 0.5927 ≈ p c and network-dimensions L × L for L between 50 to 90 links. The lowest numerically feasible |∆P | are used in order to stay in the lower linear regime in Eq. (13), specifically 2.8 Pa/link ⪅ |∆P |/L ⪅ 5.8 Pa/link. When operating at low |∆P |, the flow, which is mainly through the connected paths, stabilizes quickly and retains approximately a constant value compared to the fluctuating flow at higher |∆P |. For these simulations, the flow is driven for approximately 40 porevolumes of fluid through the network, where one porevolume is equal to the total volume of pore space in the network. The values of Q are calculated by averaging over the last 20 pore-volumes simulated. Variation in the connected paths a network can have is covered by averaging the results over 50 network realizations. The results are shown in Fig. 7, where we get t ′ /ν = 4.3 ± 1.0, giving This is a huge value. A possible explanation for the observed value is that the system is not at a critical point in spite of the geometry of the easy links and the connected paths indicating this. In our argumentation, we have not taken into account the empty links, i.e., those links that do not contain any interfaces. They will be indistinguishable from the easy links with respect to the dynamics. These empty links drive the system away from the percolation critical point, and Fig. 7 is in reality indicative of non-algebraic behavior. We have indicated this possible shift in transition in the phase diagram shown in Fig.  1.

V. NON-DARCY BEHAVIOR
In the simulations done for this section, networks have dimensions 100 × 100 links 2 . For each |∆P |, the flow is driven for approximately 100 pore-volumes of fluid through the network. This ensures the steady-state flow and the value of Q in steady state is calculated by averaging over the total flow rate during approximately the last 25 pore-volumes simulated.

A. Hysteresis
We pose here the question of whether there are any hysteretic effects from raising and lowering the pressure drop |∆P | on the volumetric flow rate Q. The result is shown in Fig. 8. With the passing of time, measuring in terms of injected pore-volumes, |∆P | applied across a network is raised and then lowered in steps. The |∆P | values used, 200 Pa, 266 Pa, 355 Pa, 473 Pa and 631 Pa, are from the lowest numerically feasible range. It can be observed from Fig. 8 that whenever |∆P | is returned to the same value, Q also quickly stabilizes back to the previous value it had with the same |∆P |. This shows that the steady state results generated using the DPN model do not depend on long-term memory [48].

B. Volumetric flow rate dependence on pressure drop
The results relating Q and |∆P | in systems with zero P t and different values of p + are shown in Fig. 9. We used p + ∈ {0.42, 0.46, 0.50, 0.54, 0.58, 0.5927} for the simulations. For each of these p + , the results were averaged over ten randomly chosen networks that have connected paths, meaning ten networks were randomly chosen from a subset of networks with zero threshold pressure P t . To assist the understanding of Fig. 9, velocity maps of a network with p + = 0.5 at various |∆P | have been plotted in Fig. 10.  to log |∆P | ⪅ 2.8 in Fig. 9, in other words |∆P |/L ⪅ 6.3 Pa/link. The transition between this regime to the next is more gradual for p + away from 0.5 in Fig. 9. In this regime with very low |∆P |, we find β = 1.00 ± 0.01. The velocity maps of a network with p + = 0.5 at two different |∆P | in this regime are shown in Fig. 10(a) and (b) and they indicate that the flow is mainly through the neutrally wet (red) links. When increasing log |∆P | from Fig. 10(a) to Fig. 10(b), the impact mainly manifests in the increase of the speed of the fluids rather than the creation of new paths. Therefore, it makes sense that the flow remains Darcy-like with β approximately equal to 1. In the lowest regime in Fig. 9, it is apparent that the mobility M in Eq. (13) decreases when p + moves away from 0.5 towards p c ≈ 0.5927 and 1−p c ≈ 0.4073. In this regime, the flow is mainly through the connected paths. The network has more connected-path links, and transports more fluid for the same |∆P | hence resulting in a larger Q, meaning a larger M when p + moves towards 0.5. For instance, at log |∆P | ≈ 2.3, the number of active connected-path links is high at p + = 0.5, as can be seen in Fig. 10(a), slightly lower at p + = 0.54, as can be seen in Fig. 11(a), and significantly lower at p + = 0.58, as can be seen in Fig. 11(b), making Q at p + = 0.58 significantly less than in the other two cases.
The middle regime in Eq. (13) seems to corresponds to 3.3 ⪅ log |∆P | ⪅ 4.1 in Fig. 9, in other words 20.0 Pa/Link ⪅ |∆P |/L ⪅ 125.9 Pa/link. Here, the exponent in Eq. (13)  The highest regime in Eq. (13) seems to corresponds to log |∆P | ⪆ 4.5 in Fig. 9, in other words |∆P |/L ⪆ 316.2 Pa/link. Here, the exponent in Eq. (13) is β = 1.00 ± 0.01 and M D is the same for all p + examined. The velocity maps taken from two different points in this regime are shown in Fig. 10(e) and (f). In both cases, almost all the links in the network are carrying flow regardless of their wettability, hence increasing |∆P | does not create new paths. The effect of capillary barriers in the links becomes insignificant in comparison to the enormous pressure drop across the links, making all p + produce the same Q at the same |∆P |. Increasing |∆P | in this regime, increases Q linearly, indicative of Darcy flow.
As the results in Section III show, there are very few to zero connected paths outside of the range 1 − p c ≤ p + ≤ p c examined in Fig. 9. If p + were very close to the range examined in Fig. 9, the behavior of β and M would have been expected to be the same as in Fig. 9 since the flow will similarly be carried by the connected paths. To test p + further away, simulations have been performed with p + = 0.2 and 0.3 and the results are shown in Fig. 12. Here, P t is not zero, unlike the systems used for Fig. 9 and corresponding constitutive Eq. (13). In this case, we find a constitutive equation where P u is the crossover pressure between non-linear and Darcy behavior. By varying P t from 0.00 Pa to the lowest |∆P | in the data sets with an increment of 0.01 Pa, mathematical linear fits with slopes β were calculated at the lowest pressures to find the candidate that gave the least room-mean-square error. This gave β = 2.23 ± 0.05 and P t = (3.4 ± 0.5) kPa for p + = 0.2 and β = 2.29±0.05 and P t = (2.0 ± 0.5) kPa for p + = 0.3. The regime these β correspond to is the middle regime discussed in Fig. 9 where the behavior was also non-linear due to the capillary barriers created by the interfaces between the two fluids. Fyhn et al. [21] has observed β > 1 behavior even in networks with the same wetting angle everywhere which would be the same as having p + → 0.0 or 1.0 here. Due to the lack of connected paths in systems with p + = 0.2 and 0.3, the lowest regime in Fig. 9 does not appear for the results in Fig. 12. Lastly, the highest regime where β ≈ 1 should occur for all where flow pushes through almost the entire network and there is almost no influence of p + . That is indeed what we see in Fig. 12 also.

VI. CONCLUSION
We studied the effect of having porous media consisting of randomly mixed dual-wettability grains on the im-miscible two-phase flow using a dynamic pore network model. The model treats the interfacial tension between the two fluids similarly to a model introduced by Irannezhad et al. [31,32]. The model has two parameters, the saturation S w and the probability p + to have a grain of "+" type. The model, which is explained in Fig. 2, contains links (pores) of three types when filled with two immiscible fluids A and B: Links that are wetting with respect to fluid type A, links that are wetting with respect to fluid type B and easy links where there are no capillary forces associated with interfaces. The parameter p + controls the number of links generating capillary forces and easy links. The model has a rich phase diagram, sketched in Fig. 1. There is a region 1 − p c ≤ p + ≤ p c , where p c is the site percolation threshold, where the easy links form connected paths across the network. Outside this region, i.e., for p + ≪ 1−p c or p + ≫ p c , easy links do not percolate. We find two classes of constitutive equations for volumetric flow rate Q vs. pressure drop |∆P |: For 1 − p c ≤ p + ≤ p c we observed a constitutive equation as in Eq. (13), see Fig. 9, whereas for p + ≪ 1 − p c or p + ≫ p c , we observed a constitutive equation (17), see Fig. 12. The crucial point that distinguishes these is whether there is a non-zero threshold pressure P t .
When 1 − p c ≤ p + ≤ p c we observed the following: At the regimes with lowest and highest |∆P |, it seems that β = 1.00 ± 0.01, because there is no significant change in the paths fluids are flowing through and increasing |∆P | only increases the flow in the already active paths. At the lowest |∆P |, the flow is mainly always through connected paths with zero resistance. When p + → 0.5 in this regime, there are more connected paths which means more fluid gets transported, making Q hence M higher. At the highest |∆P |, almost the entire network is always active. On the other hand, β > 1 in the middle regime where an increase in |∆P | increases the flow in the active paths and in addition opens new conducting paths. In the middle and the highest regimes, the flow is no longer mainly through the connected paths and the differences between the pressures across the links and the capillary barriers in the links are large. With the diminished role of the connected paths and capillary barriers at higher pressure drops, M m and M D do not depend strongly on p + . The exponent in the middle regime was found to be β = β 3 = 2.56 ± 0.05. We saw no systematic dependence of β on p + .
For p + = 0.2, however, β = 2.23 ± 0.05 and P t = (3.4 ± 0.5) kPa, and for p + = 0.3, β = 2.29 ± 0.05 and P t = (2.0 ± 0.5) kPa. Due to the necessity of determining P t and β simultaneously at these p + values, there is more uncertainty associated with the measurements of β. It is not possible to verify or falsify whether there is a fixed β = β 2 , or whether it depends on p + and S w . The existence of connecting paths is a percolation problem. They disappear when p + → (1 − p c ) + or p + → (p c ) − . It would therefore be expected that the mobility M defined in Eq. (13) would exhibit a critical behavior similar to the conductance near a critical point. By making the hypothesis that M behaves as in Eq. (14) and using finite size scaling, we determined t ′ ≈ 5.7, see Fig. 7. This is a huge value and raises the suspicion that the system is not critical where percolation theory dictates that it should be. Possible suspects for causing this push away from criticality are the links that do not contain interfaces. They are not easy links, but they have precisely the same effect on the dynamics of the flow as the easy links. If this is so, the transition lines would then be shifted as shown in Fig. 1.
We have only explored a small part of the phase diagram of this rich model in this first study. The phase diagram should be investigated in more detail and over a wider range of parameters. The nature of the transition lines is as of now unknown, and should also be further investigated. There are percolation transitions in the model. Where are they and what are their properties as the transport is not through percolation clusters?

VII. DECLARATION
Author Contributions: HF performed the numerical simulations and the data analysis and wrote the first draft of the manuscript. HF wrote the code specific for this project based on algorithms written by SS. AH and SS suggested the idea of the problem. AH worked out the relation to percolation theory. Funding: This work was supported by the Research Council of Norway through its Center of Excellence funding scheme, project number 262644.