Rheology of immiscible two-phase flow in mixed wet porous media: Dynamic pore network model and capillary fiber bundle model results

Immiscible two-phase flow in porous media with mixed wet conditions was examined using a capillary fiber bundle model, which is analytically solvable, and a dynamic pore network model. The mixed wettability was implemented in the models by allowing each tube or link to have a different wetting angle chosen randomly from a given distribution. Both models showed that mixed wettability can have significant influence on the rheology in terms of the dependence of the global volumetric flow rate on the global pressure drop. In the capillary fiber bundle model, for small pressure drops when only a small fraction of the tubes were open, it was found that the volumetric flow rate depended on the excess pressure drop as a power law with an exponent equal to 3/2 or 2 depending on the minimum pressure drop necessary for flow. When all the tubes were open due to a high pressure drop, the volumetric flow rate depended linearly on the pressure drop, independent of the wettability. In the transition region in between where most of the tubes opened, the volumetric flow depended more sensitively on the wetting angle distribution function and was in general not a simple power law. The dynamic pore network model results also showed a linear dependence of the flow rate on the pressure drop when the pressure drop is large. However, out of this limit the dynamic pore network model demonstrated a more complicated behaviour that depended on the mixed wettability condition and the saturation. In particular, the exponent relating volumetric flow rate to the excess pressure drop could take on values anywhere between 1.0 and 1.8. The values of the exponent were highest for saturations approaching 0.5, also, the exponent generally increased when the difference in wettability of the two fluids were larger and when this difference was present for a larger fraction of the porous network.


Introduction
The study of rheology of two-phase flow in porous media is pivotal for many disciplines, and the wettability conditions of the system is an important factor that directly affects the rheology. Examples for relevant disciplines include drug delivery in biology (Vafai 2010), studies of human skin behavior relevant for cosmetic and medical sectors (Elkhyat et al. 2001), creation of self-cleaning and fluid repelling materials relevant for textile industry (Li et al. 2017) and oil recovery (Kovscek et al. 1993) and carbon dioxide sequestration (Krevor et al. 2015) in geophysics (Blunt 2017;Marle 1981). All of these examples, dealing with different kinds of porous media, will benefit from a better understanding of two-phase flow under different wetting conditions. Two-phase flow means simultaneous flow of two fluids in the same space. When an immiscible fluid is injected into a porous medium filled with another fluid, different transient flow mechanisms occur depending on the flow conditions, such as capillary fingering (Lenormand and Zarcone 1989), viscous fingering (Toussaint et al. 2005;Måløy et al. 1985;Løvoll et al. 2004) and stable displacement (Frette et al. 1997;Méheust et al. 2002). After the transient flow mechanisms have surpassed, steady state sets in, which is the regime in which the rheology of two-phase flow under different wetting conditions is examined in this work.
Darcy's law is widely used to describe the flow of fluids through a porous medium which states that the volume of fluid flowing per unit area per unit time depends linearly on the applied pressure drop across a representative elementary volume in that porous medium (Blunt 2017). That is indeed the case for large applied pressures, however the linearity gets modified into a power law at the low pressure limit. For the flow to start, the applied pressure has to overcome the disordered capillary barriers (Sinha and Hansen 2012). When the applied pressure is so small that it exceeds the capillary barriers in only parts of the porous medium, the capillary forces will be comparable to the viscous forces. In this case, the volumetric flow rate scales nonlinearly with the pressure drop due to the fact that increasing the pressure drop by a small amount creates new connecting paths in addition to increase the flow in the previously connected paths. Earlier works (Roy et al. 2019;Sinha et al. 2021;Tallakstad et al. 2009a;Rassi et al. 2011;Tallakstad et al. 2009b;Aursjø et al. 2014;Gao et al. 2020a;Zhang et al. 2021) have provided experimental, theoretical and numerical evidences for this phenomena in porous media under uniform wetting conditions. Instead of assuming uniform wetting conditions, we here investigate the same phenomena using non-uniform wetting conditions, theoretically and numerically.
The wetting condition of a porous medium is a major factor controlling the location, flow and distribution of fluids (Anderson 1986), and is a result of the interplay between the attractive forces on the surface of the adjoining materials. When two immiscible fluids flow in a porous medium, the relative values of the surface tensions between each pair of the three phases, namely the fluids and the solid, determine the wetting angle and hence the equilibrium configuration of the fluids. In nature, the wettability of a porous medium tends to alter along the system and results in a range of different wetting angles. For instance, the internal surface of reservoir rocks is composed of many minerals with different surface chemistry and absorption properties, which can cause wettability variations (Anderson 1986). There are different types of non-uniform wetting conditions depending on the de-gree of non-uniformity as well as the geometrical and topological distribution of regions with different wettability. The examples include fractional wettability where grains with same type of wettability are packed together in different proportions or mixed wettability where there are continuous paths with one type of wettability (Anderson 1986;Salathiel 1973). It is often useful to make these distinctions because the physical processes which create non-uniform wetting conditions can result in different forms of connectedness. In this work, we want to study how the deviation from uniform wetting conditions affect the rheology. Hence, it is desirable to isolate the effect of non-uniform wettability in terms of the mean wetting angle and the spread of the wetting angles. To this end, we use mathematical models with wetting angles determined from various distributions. We use the term mixed wet to denote the resulting non-uniform wetting conditions, but note that this term can also imply geometrical effects mentioned above which are not considered here. We leave for future work the problem of how others types of non-uniform wetting conditions can affect the rheology further. A mechanism for a correlated wettability distribution for pore-network modeling, where the wettability depends on the connected oil paths, was demonstrated previously by some of the authors of this manuscript (Flovik et al. 2015) and may be adopted in future.
Several works in the past have investigated multiphase flow in mixed wet porous media, and discovered clear discrepancies in the fluid behavior in uniform wet systems and in mixed wet systems. Experimental studies have found that the main determinant of the filling sequence in a porous medium is the wettability rather than the pore size (Scanziani et al. 2020;Gao et al. 2020b). There were also findings from experimental studies indicating that the processes where it is necessary to allow the flow of both fluids favor mixed wetting conditions (AlRatrout et al. 2018;Alhammadi et al. 2017), such as oil recovery or fluid transport through membranes or in biological tissue. These experimental findings show the importance of understanding the effect of wettability even further, which is easier to do through analytical and numerical studies where large range of wetting conditions can be examined in short time. In the papers by Sinha et al. (2011) andFlovik et al. (2015), pore network models similar to the one used in the present article were used to investigate the effect of wettability alteration due to changes in salinity in oil-brine mixtures. The wettability alterations were done by changing between either complete wetting and complete non-wetting conditions in the first article (Sinha et al. 2011), and by changing the wetting angles continuously between two limits depending on the cumulative flow of the wetting phase in the second article (Flovik et al. 2015). The results from both show that local alterations of the wettability introduce qualitative changes in the flow patterns by destabilizing the trapped clusters. While such past numerical studies provide important insight into the behaviors of mixed wet porous media and support the usefulness of mixed wettability, they consider limited cases of the wetting angle conditions and do not consider the effect of the applied pressure on the flow. In the present work, we conduct a systematic analysis of the effect of mixed wetting conditions, both in terms of a wide range of different mean wetting angles as well as different spread of the wetting angles. In doing so we manage to perform a direct study of the relation between the total volumetric flow rate and the pressure drop across the system as influenced by the mixed wettability.
Stated more in detail, the investigations in this work have been carried out by, firstly, calculating the total volumetric flow rate in a model consisting of a bundle of capillary tubes with mixed wet properties (Roy et al. 2019;Sinha et al. 2013). Thereafter, case studies with various specific wetting angle distribution have been performed through numerical calculations which confirmed the analytical results in addition to providing a holistic picture. Secondly, mixed wetting conditions have been implemented into a dynamic network model (Sinha et al. 2021) where the motion of the fluid interfaces are followed through the porous medium. The results confirm that the volumetric flow rate Q indeed depends on the applied pressure drop ∆P as where P t is the minimum pressure drop necessary for flow. The exponent β > 1 in the low pressure limit and β = 1 in the high pressure limit. More specifically in the low pressure limit, the capillary fiber bundle model considering a simple system gives β = 2 and β = 3/2, while the dynamic pore network model considering a more sophisticated system gives values varying anywhere between β ∈ [1.0, 1.8] depending on the system wettability configuration. The models and the wetting condition implementing methods, as well as previously existing relevant theories, are explained in Sec. 2. The theoretical and numerical results are presented and discussed in Sec. 3 and a conclusion summarizing the findings is given in Sec. 4.

The Capillary Fiber Bundle Model Description
The first model that is used to investigate immiscible two phase flow in mixed wet porous media is a capillary fiber bundle (CFB) model (Roy et al. 2019;Sinha et al. 2013). This model consists of a bundle of parallel capillary tubes, disconnected from each other, each carrying the two immiscible fluids. A typical porous medium normally has a varying radius for the links in the system, which is a factor that contributes to the capillary pressure being position dependent. To emulate this effect, sinusoidal shaped tubes have been used in this model. A sketch of the model is shown in Fig. 1.
As the main goal of this work is to examine the effect of mixed wettability, each one of the tubes in the CFB model has been given a wetting angle θ chosen randomly from a certain predefined distribution ρ(θ). This means that each tube has the same assigned wetting angle over its entire length. The flow is driven by applying a global pressure drop ∆P over the system. The total global volumetric flow rate Q of the bundle of tubes is then calculated by considering the contributions to the flow given by each tube. This calculation has been carried out both analytically and numerically.

The Dynamic Pore Network Model Description
The second model that is used for the investigations in this work is a dynamic pore network (DPN) model, a complex numerical model which is not analytically solvable (Sinha et al. 2021;Aker et al. 1998;Knudsen et al. 2002;Tørå et al. 2012;Gjennestad et al. 2018). A sketch of the network used in the model is given in   The mixed wettability has been implemented into the DPN model by randomly choosing a wetting angle θ for each one of the links in the network from a certain predefined wetting angle distribution ρ(θ). As in the case with the first model, the flow in the DPN model is driven by a pressure drop ∆P over the system. When using periodic boundary conditions, ∆P is defined across a period of the system. The total volumetric flow rate Q is constant over all the cross sections normal to the direction of the overall flow, as the one illustrated with a horizontal line in Fig. 2.

Commonalities
There are several commonalities in the two models. The smallest computational unit, which will be denoted SCU for ease of reference, in the CFB model is a single tube and that in the DPN model is a single link. Even though each SCU in the two models has uniform wettability, the entire system consisting of various such entities together describes a mixed wet porous media. In both models, the radius of each SCU i, with cylindrical symmetry around the center axis x has the form where r 0,i is a constant and a is the amplitude of the periodic variation. In the DPN model, l in Eq. (2) is the length of the links, since the shape of the links covers only one period of oscillation, giving an hourglass form as shown in Fig. 2. In the CFB model, l in Eq.
(2) is the wavelength of the shape of the tubes, as shown in Fig. 1. The flow within SCUs are governed by the following equations. In SCU i, the capillary pressure p c,i (x) across an interface between the two immiscible fluids with wetting angle θ i can be derived from the Young-Laplace equation to be (Blunt 2017) where σ is the surface tension. For each SCU with length l experiencing a pressure drop ∆p across its body, the fluid within it is forced to move due to the force exerted by the total effective pressure. Total effective pressure is the difference between ∆p and the total capillary pressure k p c (x k ) due to all the interfaces with positions x k ∈ [0, l ]. Assuming that the radius does not deviate too much from its average valuer i , the volumetric flow rate q i in SCU i is given by (Sinha et al. 2021;Washburn 1921) where µ i is the saturation weighted viscosity of the fluids given by Here, s A,i = l A,i /l and s B,i = l B,i /l are saturations of the two fluids A and B with viscosities µ A and µ B and lengths l A,i and l B,i . In the DPN model, l is the same as l from Eq. (2). In the CFB model, l is the length of the whole tube. The capillary number Ca, which is the ratio of viscous to capillary forces, is related to q i through Ca = µQ/(σα) where Q is the sum of all q i through a cross sectional area α (Sinha et al. 2021). Note that due to the incompressible nature of the fluids examined in this work, q given by Eq. (4) is the same for any position along a single SCU. Also note that all θ in this work are defined through fluid A, as shown in Fig. 3, which means the wetting angles of fluid B are 180 • − θ. The fluid that makes the smallest angle with the solid wall is the wetting fluid in that region of the pore space and the other fluid is the non-wetting fluid. 3 Results

The Capillary Fiber Bundle Model Results
The analysis of flow in a single tube is presented in Subsec.

A Single Tube
In the paper by Sinha et al. (2013), they calculate the flow properties in a capillary tube with cos(θ) = 1. Here, for a single tube, we will follow their calculations while keeping θ as a variable as it is needed for the rest of the work. The parameters r 0 , µ, l, σ, a, l and ∆p, as given in Eqs.
(2) to (5), are kept constant for all the tubes. All the tubes have the same length l = L and global applied pressure ∆p = ∆P . We start by considering a capillary tube with N bubbles. A "bubble" is one type of fluid restricted on two sides by the other fluid. Each bubble j has the center of mass position x j and a width ∆x j . From Eqs.
(2) to (5) we find that the volumetric flow rate through one tube is where δx j = x j − x 0 . Due to the incompressible nature of the fluids, the velocity of the bubbles is approximately constant along the axis of flow and equal to dx 0 / dt ≈ dx j / dt ≈ q/(πr 2 ). In addition, the effect of the variation in ∆x j and δx j can be assumed to be small. With this, Eq. (6) can be rewritten as where and With algebraic manipulations, Eq. (7) can be rewritten as We wish to calculate the average velocity of the bubbles as they travel from one end to the other end of a tube segment with length l using a time T , T can be calculated by using the equation of motion in Eq. (13), Inserting the result in Eq. (15) into Eq. (14) and using the relation dx/ dt ≈ q/(πr 2 ) gives the average volumetric flux equation where Θ is the Heaviside step function. From Eq. (16) we see that, on average, the direction of flow is opposite to the pressure drop, as expected. Additionally, we see that for a nonzero flow, ∆P needs to exceed a certain threshold γ that is specific for the tube.

A Bundle of Tubes
In the CFB model, the global volumetric flow rate Q of a bundle of tubes is the sum of the time-averaged individual volumetric flow rates q of all the tubes that carry flow. As remarked at the end of Subsec. 3.1.1, the tubes that carry flow are those that have a threshold γ that satisfies the requirement |∆P | > |γ|. We will define a quantity P t which is the minimum possible γ a tube can have. This means that the first active path across the entire system occurs once ∆P exceeds P t . Let us also define γ max as the maximum possible γ a tube can have, for later use. The factors that determine γ can be seen from Eq. (11). Among those, θ is the only variable that varies from tube to tube, while the other quantities are set to be universal. Under a constant ∆P , what determines which tubes in the bundle will conduct flow, while others do not is therefore their θ. Using Eq. (11) and that θ ∈ [0 • , 180 • ], the requirement for flow to happen in a tube can be rewritten as with θ a = arccos min 1, Note that θ + P t = θ − P t when |P t | = 0. Since the flow requirement in Eq. (17) indicates the range of θ that give the system a nonzero flow, Q can be expressed as a function of the probability distribution ρ(θ) of the wetting angles through Inserting Eq. (16) into Eq. (21) gives Case studies with several different forms of ρ(θ) will be done numerically in Subsec. 3.1.3. Here, we will solve Eq. (22) for a general ρ(θ) and show that the exponent Here, |γ max | is the maximum possible threshold pressure a tube can have.
Case 1 : |∆P | − |P t | |γ max | Using that |∆P | is large in this case, we can write Inserting Eq. (25) into Eq. (22) and using that the distribution ρ(θ) is normalized to 1, gives In the last step, we have used that |P t | which is the minimum |∆P | needed to achieve Q > 0, is a much smaller number than |∆P |. For the equations derived for all three cases in Eq. (24), we wish to express Q in terms of |∆P | − |P t | for ease of comparison with Eq. (23). Comparing Eq. (26) with Eq. (23), one gets β = 1 for case 1.
Common for Case 2 and Case 3: Eq. (24) states that the effective pressure obeys |∆P | − |P t | |γ max | in case 2, while it obeys |∆P |−|P t | |P t | in case 3. From Eq. (11), the threshold pressure γ is so that |γ| ≤ k γ . This criterion should also be followed by the maximum possible |γ| which is |γ max | and the minimum possible |γ| which is |P t |. Combining these information, a common requirement for cases 2 and 3 should be Eq. (18) then becomes where Eq. (27) was used in the last step. Next, based on Eq. (27), Taylor expanding the integrands of Eq. (22) with respect to θ around θ ± P t gives Performing this integrations using the integration limits given in Eqs. (20) and (28) and that θ b = 180 • − θ a (Eq. (19)) gives Eq. (30) holds true for both cases 2 and 3.
In the transition region between case 1 and cases 2 and 3, where |∆P | − |P t | ≈ k γ , the volumetric flow depends more sensitively on the wetting angle distribution function ρ(θ), and is in general not a simple power law. Nevertheless, we can use the analysis presented here to compute the height of that transition region. Taking the logarithm of Q we find that where η is a constant that is the same for all the three cases. Evaluating this at |∆P | − |P t | = k γ we see that the height difference between case 2 and case 1 in a logarithmically scaled plot is and the height difference between case 3 and case 1 is This is shown in Fig. 4. It is assumed in the above analysis that ρ θ + P t +ρ θ − P t < 1 − (P t /k γ ) 2 , since ρ(θ) is a normalized probability distribution. Hence, h < 0, which means that the β in the transition region must be larger than outside the transition region.

log(k γ )
-h β = 1 β >2 β = 2 or 1.5 Fig. 4: Height h of the transition region is defined as the distance between the two lines that are extrapolations of the linear parts of the curve at ∆P − P t = k γ . Here, Q is the volumetric flow rate and ∆P − P t is the excess pressure.

Numerical Results
Here, we study numerically the volumetric flow rate Q's response to a wide range of an applied pressures ∆P . In addition to verify the analytical results from Subsec. 3.1.2, this numerical study allows us to probe the transition region where ∆P is the same order of magnitude as γ max . We present results with three different normal and uniform distributions of θ, and note that we find similar results also for many other distributions.
The CFB model has also been numerically tested with ρ(θ) being normal distribution with meanθ and standard deviation δθ given by The results from when θ is normally distributed with (θ = 40 • , δθ = 10 • ), (θ = 40 • , δθ = 30 • ) and (θ = 90 • , δθ = 30 • ) are shown in Fig. 5b. These results once again confirm the analytical calculations performed in Subsec. 3.1.2. Toward the right in Fig. 5b where |∆P |−|P t | |γ max |, case 1, all three examples follow β = 1. Toward the left in Fig. 5b where |P t | |∆P | − |P t | |γ max |, the results follow β = 2, case 2. Notice that, P t = 0 for all these three cases. Since a nonzero Q occurs only when |∆P | > |P t |, the requirement for β = 1.5, namely |∆P | − |P t | |P t |, case 3, is not satisfied here. In Fig. 5b, the transition region between when β = 1 and when β = 2 exhibits a gradient β > 2, in accordance with the analysis above. The same effect can also slightly be seen in Fig. 5a, but is much more apparent in Fig. 5b. One can also see from Fig. 5b that β is larger for smaller δθ and meansθ further away from 90 • . This can be understood from Eq. (34), since smaller δθ and larger θ − 90 • implies smaller ρ(90 • ). From Eq. (34), we see that smaller ρ(90 • ) implies a larger height difference h. Physically, this can be understood from the fact that a smaller ρ(90 • ) means that a smaller fraction of the total number of tubes are active in the low pressure regime. Here, ∆P is the applied pressure and P t is the minimum threshold pressure. The flow rates have been normalized with Q 0 = −(πr 4 sgn (∆P ) k γ )/(8µL) and the pressures have been normalized with k γ . Straight lines with gradients β has been added.

The Origin of β
We can write volumetric flow rate as Q = Nq, where N is the number of open tubes andq is the average flow per open tube. Put differently,q is the average of q for all the open tubes. We propose that β can be understood from how each of these factors change with applied pressure ∆P . Suppose that increasing the pressure difference from |∆P 0 | − |P t | to |∆P | − |P t | =x (|∆P 0 | − |P t |) transforms the number of open tubes and the flow per tube according to N 0 →xāN 0 and q 0 →xbq 0 , respectively. In this case we see that the volumetric flow change from so β =ā +b. Consider firstq. From Eq. (16), we know that the volumetric flow for a tube with threshold pressure γ is q ∝ γ. Hence, if most of the active tubes have threshold pressure just below the applied pressure (case 3), thenā = 1/2. On the other hand, if most of the active tubes have threshold pressure well below the applied pressure (cases 1 and 2), thenā = 1.
Next, consider the number of active tubes transporting fluid, N. This is given by where N max is the total number of tubes in the system. When the applied pressure is larger than the maximal threshold pressure γ max , then all the tubes are active and N = N max . Thus, for case 1, we haveb = 0 and consequently β =ā +b = 1. On the other hand, when |∆P | − |P t | |γ max |, then θ 1, as seen from Eq. (28). Thus, sob = 1. Combining this with the result forā we see that in case 2, we get β =ā +b = 2, and in case 3, we get β =ā +b = 1.5. This explains why β ∈ {1, 1.5, 2} when either all tubes are active (case 1) or only a small fraction is active (cases 2 and 3).

Data Collecting and Processing Procedures
Using the method described in Subsec. 2.2, numerical simulations of the DPN model have been performed. The following factors and parameters have been kept constant during all simulations. The 2D network used was made of 64 × 64 links and had periodic boundary conditions in all directions. All the links had length l = 1 mm, average radiir ∈ [0.1l, 0.4l] and amplitude of the periodic variation a = 1 mm, see Eq.
(2). The viscosities of the fluids A and B were µ A = µ B = 0.01 Pa·s, see Eq. (5). The surface tension between the fluids were 0.03 N/m, see Eq.
(3). In this closed network system, the control parameter was the saturation of one of the fluids in the whole system, which was tested for values S A ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}.
The distributions of the wetting angles, ρ(θ), have firstly been tested for a uniform distribution with θ ∈ [0 • , 180 • ]. Secondly, normal distributions with meansθ ∈ {0 • , 30 • , 60 • , 90 • } and, with each mean, standard deviations δθ ∈ {0 • , 30 • , 60 • } were also tested. The ρ(θ) mentioned until now were implemented into the network with all of the above-mentioned values of S A . In addition to this, at S A = 0.5, several more normally distributed θ were examined. They were with meansθ ∈ {0 • , 30 • , 60 • , 90 • } and, with each mean, standard deviations δθ ∈ {15 • , 45 • , 75 • }, as well as, with meansθ ∈ {15 • , 45 • , 75 • } and, with each mean, standard deviations δθ ∈ {0 • , 30 • , 60 • }. Note that the normally distributed θ could go outside the interval [0 • , 180 • ], which is equivalent to a slightly increased weight around 0 • or 180 • because the angle only comes in through cos θ, as seen in Eq. (3). Another thing to note is that only distributions withθ ≤ 90 • have been considered. This is because since the fluids have the same viscosity, a symmetry is in place where the case with meanθ and saturation S A is the same as the case with mean 180 • −θ and saturation 1 − S A .
For each S A , the system was driven by various different Q, and for each Q, 20 different realizations of the network were performed and averaged over. The global applied pressure ∆P in the direction of the flow was measured, and was calculated by averaging over the fluctuations after the system had reached a steady state. After obtaining a set of Q and ∆P for every ρ(θ) at every S A , data analysis had to be performed to determine the global threshold pressure P t below which there is no flow through the whole system, as well as the exponent β in Q ∝ (|∆P | − |P t |) β . Note that similar to the CFB model, the first active path across the entire system in DPN occurs once ∆P exceeds P t .
The process of determining P t and β started with deciding the indices of the data points that belonged to the linear and power law regimes through visual examination. As shown with an example in Fig. 6, this meant deciding the indices of the datapoints that lied between n start β≈1 to n end β≈1 and belonged to the region with β ≈ 1, as well as the indices between n start β>1 to n end β>1 that belonged to the region with β > 1. The error bars were calculated as the absolute values of the difference between the results and the results that would have been if the range of data points included from each region were reduced. The next step was to perform linear fitting of Q 1/β against ∆P on the data points from n start β>1 to n end β>1 . The linear fitting was of the form c 1 Q 1/β + c 2 with c 1 and c 2 real numbers. Due to the definition that P t = ∆P exactly when Q 1/β becomes nonzero, P t = c 2 . This procedure was repeated for a range of different βs, and the P t that gave the least root-mean-square error was chosen as the final candidate. Thereafter, linear fitting log(Q) versus log(∆P −P t ) separately for the data points from index n start β≈1 to n end β≈1 and data points from index n start β>1 to n end β>1 gave the values of β in those regions.  Fig. 6: Linear fitting the data points from index n start β≈1 to n end β≈1 and data points from index n start β>1 to n end β>1 separately gives β in those regions.

Simulation Results
All of the simulations performed, using the parameters and the different wetting angle distributions ρ(θ) described in Subsec. 3.2.1, resulted in a Darcy-like flow with β ≈ 1 in the high pressure limit where most of the links were active. The transition region with β > 2 flow, as in the case with the CFB model, was not observed with the DPN model. Note that the DPN model, compared to the CFB model, simulates a porous medium that has a more complex interplay of the fluids in the links which separate and rejoin. It could be speculated that this advanced behavior of the network eliminates the transition region originally observed in the CFB model, in other words, the transition region may be an artifact of the CFB model. In the low pressure limit result of the DPN model, the exponent β shows dependence on the saturation and the wettability properties of the network. This is also the case with the threshold pressure of the network P t . A closer exploration of the latter two factors will now be presented.
The results for β in the low pressure limit are shown in Fig. 7 and takes on various values in the range β ∈ [1.00 ± 0.05, 1.82 ± 0.05]. The phenomenon with β > 1 originates from that many links in the network are not yet opened in the low pressure regime, which means increasing ∆P increases the number of active links in addition to increasing the flow within each active link. The overall combined effect of these allows the volume of fluid transported to rise much more than if all the links were already open. This is the same as in the capillary fiber bundle model, but here β takes on a larger range of values depending on the saturation and the wetting conditions.
The results for P t are shown in Fig. 8. The exponents β in Fig. 7, as well as, the minimum threshold pressures P t in Fig. 8 have a tendency to be largest for saturations around 0.5 and decrease steadily with increasing saturation of either one of the fluids. The reason is that when one of the fluids dominates the system, S A → 0.1 or 0.9, it is easy for that dominating fluid to create an active flow-path through the system. This is because those connected links that contain the same fluid will not experience a interfacial capillary pressure barrier. This decreases the overall threshold P t of the system, which is the cumulative effect of the interfacial capillary barriers in the network. There will be few new links to become active as ∆P increases under these circumstances which will further make Q less reactive toward changes in ∆P , meaning a decreased β. In contrast, when there are comparable amounts of the fluids A and B in the system, S A → 0.5, ∆P has to overcome the cumulative capillary pressure barrier created by the large number of interfaces between A and B. This naturally has an increasing effect on P t . When ∆P is increased under such conditions, the requirement for nonzero flow for many links are satisfied at once, causing a drastic increase in Q as a response, which increases β. Lastly regarding the effect of saturation, in the three middle rows in Fig. 8, the maxima of the P t plots are skewed to the left of S A = 0.5. In those cases, θ is concentrated around aθ < 90 • , which makes fluid A is the most wetting fluid while B is the most non-wetting fluid. It is easier for the wetting fluid to get transported in a porous medium, meaning when the saturation of the wetting fluid is lower, S A < 0.5, the system will require a higher applied pressure to achieve a non-zero flow making the system's P t higher.
The variations in β for different ρ(θ) are more subtle than in the case with P t . To get a clear overview of the differences, the maximum values, β(S A = 0.5), have been plotted as functions ofθ and δθ in Figs. 9a and 9b, respectively. Both the exponents β in Figs. 9a and 9b and the minimum threshold pressure P t in Fig. 8 vary withθ and δθ of the normal distributions. The interfacial capillary pressures, given by Eq. (3), increase with the distance between the wetting angles and 90 • . This happens with increasing δθ for wetting angles with mean aroundθ ≈ 90 • , or with decreasing δθ forθ deviating from 90 • . Note that a larger δθ means that the wettability is allowed to deviate more fromθ. This reflects in the values of P t in Fig. 8 where the peaks of the plots in row 2 and 3 decrease from left to right while the peaks in row 5 increase from left to right, and the peaks decrease from top to bottom. The same effect also creates the trend of decreasing β asθ → 90 • in Fig. 9a. When links have a wetting angle θ close to 90 • , many links will open at very small pressure ∆P , which means increasing ∆P in the typical low pressure limit does not open many new links, hence raises Q with a small β. Fig. 9b also supports this phenomena. Fig. 9b, in addition, shows the expected result that as δθ increases, β for the various normal distributions approach a value close to that of the uniform distribution with θ ∈ [0 • , 180 • ]. Instead of having a small range of wetting angles around aθ, uniform distribution provides wetting angles anywhere between 0 • to 180 • with equal probability. Therefore, it makes sense that uniform distribution results in β and P t values most similar to those of normal distributions with largest δθ which have the most variation in the wetting angles.
From the results presented here we see that both P t and β depend on the wetting angle distribution ρ(θ). In particular, they are affected not only by the mean wetting angle, but also by the spread of wetting angles. Thus, in order to fully characterize the flow through a porous media it is in general insufficient to assume a uniform wettability.

Conclusion
We studied systematically the effect of mixed wetting conditions on the effective rheology of two-phase flow in porous media by using the capillary fiber bundle (CFB) model and the dynamic pore network (DPN) model. Although the two   role in the transition regime between low and high ∆P limit where most of the tubes opens. In the DPN model, the whole process leading up to opening of all the possible links produces a Q that depends on the wettability of the system. Hence, the studies carried out in this work show that the behavior of immiscible two-phase flow in porous media changes when we move from uniform to mixed wet conditions. The wettability distribution of the porous media is therefore an important factor that should be taken into account when studying the rheology in porous media. Future works may study spatial correlation in the wettability distributions, as well as, the effect of varying viscosities of the fluids, which have not been done in this study.
From the CFB model, we found that the exponent β in Q ∝ (|∆P | − |P t |) β is 1, 2 and 3/2 in the high, low and very low effective pressure limits, respectively. The numerical solutions in addition revealed the β > 2 behavior in the transition region between these extreme limits, which was due to the rapid opening of tubes in that pressure range. In this model, the functional form of the wetting angle distribution ρ(θ) has the largest effect in the transition region, since this is the region where most tubes become active.
In the DPN model, on the other hand, we found that ρ(θ) is influential for the behavior at all pressure drops below the Darcy regime. The transition region with β > 2 flow could not be observed with the DPN model, leading to the speculation that the transition region could have been an artifact of the CFB model. In the low pressure limit, β had values varying anywhere between β ∈ [1.00±0.05, 1.82±0.05]. Both β in this pressure limit and the threshold pressure P t showed the tendency to be largest for saturations around 0.5 and decrease steadily with increasing saturation of either one of the fluids. The reason is that when there is comparable amount of both fluids in the system, there will be large number of interfaces in the system giving large interfacial capillary pressure for the link. This works to increase P t as well as making Q more reactive toward changes in ∆P , which increases β. Finally, we found that β generally increases when the difference in wettability of the two fluids is larger, and when this difference is present for a larger fraction of the porous network. This is because a larger difference in wettability, meaning that the wetting angle is further away from 90 • , gives rise to a larger interfacial capillary pressure and the overall threshold pressure. This in turn makes the effect of opening new pathways more prominent.

Declaration
This version of the article has been accepted for publication in Transport in Porous Media, after peer review but is not the Version of Record and does not reflect postacceptance improvements, or any corrections. The Version of Record is available online at: http://dx.doi.org/10.1007/s11242-021-01674-3. Author Contributions: HF developed the theory, performed the analytical and numerical calculations of the CFB model, contributed to editing of the code of the DPN model, ran the simulations and analyzed the data of the DPN model and wrote the first draft. SS suggested the idea of the problem and developed the code for the DPN model. SR sketched the initial calculation related to the CFB model and helped in data analysis. All the authors contributed in developing the theory and writing the manuscript to its final form. Funding: This work was partly supported by the Research Council of Norway through its Center of Excellence funding scheme, project number 262644. SS was partially supported by the National Natural Science Foundation of China under grant number 11750110430.

Conflicts of interest:
The authors declare no conflict of interest.