Role of pore-size distribution on effective rheology of two-phase flow in porous media

The flow of immiscible fluids inside a porous medium shows non-linearity in the form of a power law in the rheological properties of the fluids under steady state flow conditions. However, different experimental and numerical studies have reported different values for the exponent related to this power law. Here we explore how the rheological properties of the two-phase flow in porous media depends on the distribution of the pore sizes and how it affects the power-law exponent. The pore-size distribution controls fluctuation in the pore radii and their density in a porous material. We present two approaches, analytical calculations using a capillary bundle model and numerical simulations using dynamic pore-network modeling. We observe crossover from a non-linear to linear rheology when increasing the flow rate where the non-linear part is highly affected by the pore-size distribution. We have also carried out the study for different saturations of the two fluids.

vary between 2.2 (≈ 1/0.45) and 3.3 (≈ 1/0.3) [14] depending on the saturation. It was later found to converge to a certain value ≈ 2.17 when a global yield pressure is considered in the system, under which there is no flow [15]. By using pore-network modeling with 2D and 3D pore networks, Sinha et al. found the exponent to be close to 2 in the non-linear regime [15,16]. They also have reported a crossover to linear Darcy type regime at high capillary number when capillary forces are insignificant. Yiotis et al. performed Lattice-Boltzmann simulations with stochastically reconstructed porous system and studied the dynamics of fluid blobs in the presence of gravity [17]. In the steady state, they found a non-linear regime with quadratic dependence with an exponent 2, which is bounded by two linear regimes at both the high and low capillary numbers. The blobs were then studied experimentally and the non-linear exponent was found as 1.54 (≈ 1/0.65) [18]. Very recently, Gao et al. performed experiments of two-phase flow in sandstone samples. They used x-ray micro-tomography measurements and for a fractional flow of 0.5 they found the exponent in the non-linear regime to be equal to 1.67 (≈ 1/0.6) [19]. They also reported a regime with linear Darcy type behavior at lower capillary numbers where the conductance does not change significantly. Further experiments by Zhang et al. [20] explore the dependence of the exponent on fractional flow and reported values in the range of 1.35 (≈ 1/0.74) to 2.27 (≈ 1/0.44). They presented a theory that can predict the boundary between the linear regime and the non-linear intermittent flow regime.
A simple explanation for the observed power law relation between flow rate and pressure drop may be found by following the arguments of Roux and Herrmann [21], concerning the conductivity of a disordered network of resistors, where each resistor has a threshold voltage to start conducting current. If we compare the voltage across a resistor in this system to the pressure drop over a link in a network of pores, and the threshold voltage to the pressure drop necessary to overcome capillary forces due to the presence of fluid interfaces [22], we may translate the Roux and Herrmann arguments into a language appropriate for porous media. When the pressure drop ∆P across a network of pores containing fluid interfaces is increased by an amount d∆P , an additional number of pores (dN ) will start contributing to the flow. This leads to an increase in the effective conductivity of the network as more links are participating in the flow. If correlations between the opened links are ignored, the increase in conductivity will be proportional to the increase in the number of opened links, dK ∝ dN . We integrate to find where the lower integration limit is determined by the threshold pressure necessary to overcome to induce flow across the network. The flow rate is then given by Tallakstad et al. [12,13] provided another explanation by considering the scaling of clusters that are trapped by capillary forces. They assumed that the flow occurs in channels in between trapped clusters. In a two-dimensional system of length L under a pressure drop ∆P , a cluster will be trapped if the capillary force p c > λ |∆P |/L, where the λ is the length of such a cluster. The maximum length of a such a trapped cluster is therefore given by λ m = Lp c /|∆P |. By assuming the distance between the flow channels equal to the typical cluster length, the total number of flow channels will be n c = L/λ m . The total flow through all the channels is therefore the number of channels multiplied by the flow rate in each channel, which leads to Q ∝ n c |∆P | ∝ |∆P | 2 . However, if this formalism is extended to three dimensions, it leads to a cubic relationship, Q ∝ |∆P | 3 , which is in contrary to what is observed in experiments and simulations.
Sinha and Hansen [16] developed a mean-field theory for a disordered network. By analytically calculating the average rheological behavior for such a pore [22] and using Kirkpatrick's self-consistent expression for the equivalent conductivity for a homogeneous network [23], they derived the relationship, Note that the above theoretical approaches find the exponent in the non-linear regime β to be equal to 2, thus hinting at universality.
Recently, Roy et al. have studied the effect of the threshold distribution on the effective rheology of two Newtonian fluids in a capillary bundle model [24]. The model consists of a bundle of parallel capillary tubes with variable diameters along their lengths which introduce thresholds for each tube [25,26]. For power-law type distributions of the thresholds they find analytically and numerically that the non-linear exponent β can be related to α, the exponent for the power-law distribution, by the relationships β = α + 1 or β = α + 1/2 depending on whether the distribution starts from zero or has lower cut off respectively. This means, for α = 1, the uniform threshold distribution, β will be equal to 2 and 3/2 respectively for the two cases. This study clearly hints at the non-linear exponent depends on the distributions related to the system properties. We note that the capillary fiber bundle In this article, we present a detailed study how the distribution of pore sizes controls the effective rheology of the two-phase flow in the steady state. We will first study analytically the capillary bundle model, as it is an analytically tractable model for two-phase flow and provides deeper understanding of the underlying physical mechanism. There we will show that there exists a transition point as the applied pressure is increased below which the relation between flow rate and pressure drop is non-linear. We will investigate how the degree of such non-linearity depends on the shape and the width of the distribution. Above the transition point we observe Darcy-like linear flow where the distribution of pore sizes do not have any affect on the flow equation. We will then move to numerical simulation with dynamic pore-network modeling, where the similar transition is observed. There the variation in exponent in the non-linear regime and the transition point is studied by varying three parameters: the saturation of the wetting fluid, and the span and shape related to the pore-size distribution. Finally, with a two-dimensional plane of the non-linear exponent vs the transition point, we show how the above three parameters control the effective rheology of two-phase flow.

Capillary Fiber Bundle Model
In this section, we will study the analytically solvable capillary fiber bundle model (CFBM) [25,26], which is a prototype for a one-dimensional porous medium. The present authors have recently explored [24] this model in order to explore a transition from non-linear to a Darcy-like behavior when the pressure gradient across the system is continuously increased. CFBM is a hydrodynamic analog of the fiber bundle model [27], which is a disordered system driven by threshold activated dynamics and often used as a model system to study mechanical failure under stress. volumes are therefore given by πr 2 L w and πr 2 L n , where r is the average radius of the capillary tube. The saturations in each tube are then S w = L w /L and S n = L n /L. Each tube making up the bundle contains the same amount of each fluid but with its own division of the fluids into bubbles. The total volumetric flow rate q in a capillary tube at any instant of time is given by, where |∆P | is the pressure drop across the capillary tube, p c is the instantaneous capillary pressure given by the sum of all the capillary forces along the capillary tube due to the interfaces and µ av is the effective viscosity given by µ av = S w µ w + S n µ n . Here Θ(|∆P | − p c ) is the Heaviside function which is 0 for negative arguments and 1 for positive arguments. When the pressure difference across the tube is kept fixed, the average volumetric flow rate q in the steady state can be obtained by averaging Equation 4 over a time interval, where sgn(∆P) is the sign of the argument. If the tube is assumed to have a sinusoidal variation in radius with amplitude a about an average radius r and a period of length l, then γ is given by, where and Here x j is the position of center of the jth bubble if the tube is filled with 2K + 1 bubbles. Its width is ∆x j . The surface tension times the average contact angle is σ.
We are interested in the effect due to the variation in the pore radii assuming that the average contribution due to the bubble sizes are the same for all tubes. In such a scenario, the γ can be expressed in terms of the link radius r as, where k is a proportionality constant. Equation 9 implies that the higher the radius of a tube, the lower the threshold pressure will be and the fluids will start flowing at a relatively lower value for |∆P |.
We now consider a bundle of N fibers with average radii drawn from a distribution ρ(r). We assume there to be a smallest radius r min and a largest radius r max , so that ρ(r) = 0 for r < r min and for r > r max . This means that there is a smallest threshold P m = k/r max and a largest threshold P M = k/r min among the N fibers as N → ∞.
Let us define a radius The total flow rate is then given by which combined with Equations 9, 5 and 11 give

Uniform distribution
We use a uniform distribution of r between r min > 0 and r max > r min as a first illustration. We have that Equation 12 then gives for |∆P | > k/r max We now have two possibilities, either |∆P | > k/r min , making r c (|∆P |)|∆P |/k = r min |∆P |/k. The other possibility is that |∆P | < k/r min making r c (|∆P |)|∆P |/k = 1.
We consider the |∆P | > k/r min case first. This is when there is flow in all the fibers. We get For large pressure drops |∆P | k/r min , this expression reduces to We now consider the opposite case, i.e. when |∆P | < k/r min . In this case, r c (|∆P |)|∆P |/k = r min |∆P |/k. The flow rate is then If we now assume that |∆P | − k/r max = |∆P | − P m P m , we may expand this expression in terms of |∆P | − P m , finding to lowest order

Power law distribution
We now consider a power law distribution where link radii are chosen with different probabilities depending on the slope of the distribution. The expression for ρ(r) we assume to be The uniform distribution (Equation 13) illustrated in the previous section is a special case of this distribution with α = 0. The global flow rate obtained from Equation 13 is As for the uniform distribution, we have two cases to consider: |∆P | > k/r min for which r c (|∆P |) = r min and We consider the case |∆P | > k/r min first. Then, there is flow in all the fibers and we have When the pressure drop |∆P | becomes very large, i.e., |∆P | k/r min , the integral in Equation 21 simplifies by having √ u 2 − 1 → u in the integrand, and we find The other case, |∆P | < k/r min leads to r c (|∆P |) = k/|∆P | and Equation 20 becomes We now assume that |∆P | − P m P m . We expand the integral in Equation 23 to find to lowest order in ∆. Hence, the total flow rate is to lowest order in |∆P | − P m We see that the exponent α does not affect the exponent β for |∆P | larger than but close to P m . Furthermore, for α = 0 we retrieve Equation 18.
Here P m = k/r max plays the role of a threshold pressure P c below which there is no flow, whereas P M = k/r min is the crossover pressure at which there is flow in all fibers. This pressure (P M ) signals the transition to the Darcy-type flow where Q is proportional to |∆P |.
In gist, the analytical calculations with the capillary bundle model show that as soon as the pressure drop over the fiber bundle is large enough for flow to start, it enters a non-linear regime where the total flow rate Q is proportional to (|∆P | − P m ) 3/2 irrespective of the exponent α, defined in Equation 19. When the threshold pressure P M is crossed, the non-linear behavior subsides and we find Darcy-type flow. Mobility is sensitive to the details of the radius distribution for low flow rates and insensitive at higher flow rates. We like to point out here that, unlike the radii distribution, when a distribution of pore thresholds is considered, a quadratic regime with β = 2 can also obtained for CFBM when the distribution has no lower cutoff [24]. With any non-zero lower cut-off in the thresholds, β is also 3/2 there. In the present study we considered the distribution of pore radii rather than the thresholds, and for any finite pore radii, there is always a lower cutoff in the thresholds.

Dynamic Pore Network Model (DPNM)
Pore-network modeling is a computational technique to simulate two-phase flow in porous media where the porous matrix is represented by a network of pores with simplified geometries. The fluid displacements inside the pores are governed by equations for fully developed flow. We use a dynamic pore-network model where the menisci positions between the fluids track the flow [29,30]. All the pore space in this model are assigned to links and the positions where the different links meet are indicated by nodes. The flow rate in such a link is given by [28], where ∆p j is the local pressure drop across the j th link. The terms l j , g j and µ j respectively define the length, mobility and effective fluid viscosity related to that link. If µ w and µ n are the wetting and non-wetting viscosities respectively, then µ j = s n,j µ n + s w,j µ w , where s n,j and s w,j are wetting and non-wetting saturations inside that link. In this study, we consider links with circular cross sections with radii r j , for which g j = a j r 2 j /8 where a j = πr 2 j , the cross-sectional area [31]. The interfacial pressure due to surface tension between the fluids is indicated by p c which is summed over all the interfaces inside the link j, taking into account the direction of the capillary forces.
The links here represent the total pore space that consists of a pore throat in between pore bodies and the variation in the link radii along its length is therefore modeled by a sinusoidal periodic shape. The interfacial pressure at a meniscus inside such a pore can be expressed by a modified Young-Laplace equation [22], where x ∈ [0, l j ], the position of a meniscus inside the j th link. Here γ and θ are the surface tension and the contact angle for the set of fluids and pores which are kept constant throughout the simulations. We considered σ cos(θ) = 0.03N/m here. We study in-compressible fluids and therefore at every time step ∆t we have for each node i from Kirchhoff law, q i = 0, where the sum is over all the links connected to i th link. This, together with Equations 26 and 27 construct set of linear equations for every node. We solve these equations by conjugate gradient method [32] and determine the flow rates q j in every link. All the menisci in each link are then advanced by an amount ∆x j = q j ∆t. Further technical details related to the menisci displacements can be found in [30]. We consider periodic boundary conditions that lead the system to evolve to a steady state.
The network we consider in this study consists of 64 × 64 links in two dimensions (2D) which form a diamond lattice. All the links are therefore at an angle 45 • with respect to the overall flow direction. The links have equal lengths, l j = 1mm, and the disorder appears in their radii r j . We choose the values of r j from a distribution ρ(r j ) with a power α and in the range r min to r max given by the same type of distribution as in Equation 19, We denote the width of the distribution by δ = r max − r min . The power α generates different distribution types, α = 0 corresponds to a uniform distribution whereas positive and negative values of α imply higher probability to find narrower and wider links respectively.

Numerical results from DPNM
We perform simulations with constant global pressure drop ∆P and evolve the systems to steady states, where the macroscopic quantities fluctuate around a steady average. In the steady state, we measure the total flow rate Q.
Results are averaged over 20 different samples of the pore network. By fitting the numerical data with Equation 3 for the low capillary number regime, we calculate the exponent β and the threshold pressure P c . As there are two parameters to determine from each data set, we considered a method of minimizing the error related to the least square fit [16]. We illustrate this procedure briefly here. First we choose a trial value of β and perform least square fitting with the data points. This will provide a value of P c and an error associated with it. We perform this for a set of β values and find the corresponding error values. We then plot the errors as a function of β. This is shown in figure 2 for five different saturations where we see non-monotonic behavior of the errors with a minimum. We then consider the value of β that corresponds to the minimum error and use the corresponding value of P c . When plotting log(∆P − P c ) with log Q, we find the crossover point between a non-linear and a linear regime from eye approximation. From this crossover, we then identify the crossover pressure P t . In the following we present the results showing how the steady-state rheology depends on the three parameters, We will focus on the non-linear exponent β and the pressure P t related to the cross-over from non-linear to linear regime. In a P t vs β plane we will then highlight how the two quantities vary when we change the values of S w , α and δ.

Effect of S w
The variation of the flow rate Q with the pressure drop ∆P is shown in Fig. 3  All the plots show a non-linear regime at lower pressure drops and then a crossover to a linear regime. However, the slope in the non-linear regime is much higher for S w = 0.5 than 0.1 or 0.9. Similarly, the threshold pressure P c is also higher for S w = 0.5. This indicates the fact that as S w → 0 or S w → 1, the two-phase flow essentially approach to the single phase flow and it should eventually follow the linear Darcy law without any threshold pressure. In order to see the nature of this variation towards the linear regime, we plot in Fig. 3 (d), (e) and (f), the values of β and P t and P c as a function of S w , where S w is varied from 0.1 to 0.9 in the interval of 0.1. All three plots show a peak around S w = 0.5 with β = 1.97 and then continuously decrease in both sides. The decrease in both β and P t suggests that non only the non-linearity gradually disappears as S w deviates from 0.5, but also the linear Darcy regime can be obtained at a relatively lower pressure drop. At the same time, the overall threshold pressure P c also decreases as the saturation approaches to 0 or 1.
If we consider that the flow occurs in channels with capillary barriers [21], that is, the increase in Q with ∆P is contributed from two factors, the increase in the number of conducting flow paths and the increase in the flow in each path, then that explains the reduction in both β and P t as S w → 0 or 1. As the saturation of a certain fluid decreases, either it reduces the number of fluid-fluid interfaces or produces smaller bubbles of one fluid. In both the cases the effective capillary barriers corresponding to any flow path decreases. This results in more flow paths with one fluid or with negligible capillary barriers which will start flowing as soon as any pressure drop applied, making the β to move towards 1. At the same time, the maximum capillary barrier that the model needs to overcome to make all possible paths flowing also decreases, which moves the non-linear to linear transition point to a lower value of ∆P . Experimental observation of the variation of β with saturation was first reported in [14]. No threshold pressure was considered in that study while analyzing the results in that study. A recent experimental study [20] explores the variation of β and the crossover point as a function of fractional flow, F w = Q w /Q. By balancing the surface energy to create fluid meniscus to the injection energy, they developed a theory that can predict the crossover point between the two regimes they have studied. However, they have studied the crossover from the linear regime  α < 0 is intuitive, since the wider pores are higher here which will cause less capillary barrier to start the flow. The decrease in P c for α > 0 is rather counter intuitive and may be related to the decrease in slope in the non-linear regime.
In Fig. 4 (d), we plot β as a function of α which shows the decreasing trends of β on both sides of α = 0 and then becomes constant after it reaches to ≈ 1.5 around |α| = 1. This indicates that any change in the fluctuations among the link radii than the uniform distribution causes slower increase in the conductive paths when increasing the ∆P . The crossover pressure P t is plotted in Fig. 4 (e) which shows that P t remains constant around 27KPa independent of the value or sign of α. Fig. 4 (f) shows that P c shows a maximum at S w = 0.5 and decreases on either side. The decrease of P c for negative α is understandable as there will more links with larger radius in this case. For positive α we have less links with large radius. P c is still observed to decrease, most probably because of the nature of the fitting since β decreases here.
Notice that, such a variation in the non-linear exponent β while varying α was not observed in case of CFBM where β has a value of 3/2 irrespective of the radii distribution (Equation 25). This indicates that the mixing of the fluids at the nodes in a pore network has more complex effect than the flow in individual channels in CFBM. The crossover point P t here is analogous to the maximum threshold P M in CFBM which do not depend on α in both the models as it only depends on the span of the distribution and not on the the shape. Therefore, as we will see in the next section, P t has a strong dependence on δ, the span of the distribution.

Effect of δ
We now perform simulations by varying the width of the radii distribution give by, δ = r max − r min = 0.1, 0.2, 0.3, 0.5 and 0.7, where r min was kept constant at 0.1. The wetting saturation and the distribution power are kept constant here at S w = 0.5 and α = 0 respectively. The results are illustrated in Fig. 5. Interestingly, unlike the variation of β with the distribution power α as seen before, here the slopes in the non-linear regime are almost independent of the distribution width δ and remains constant around 2.0 ( Fig. 5 (a)). The crossover pressure P t , on the other hand, decreases with increase in δ, which was almost constant when we varied the distribution power α.
We also found that the global threshold pressure P c gradually decreases with increasing δ, from P c ≈ 2.5KPa at δ = 0.1 to ≈ 1.5KPa at δ = 0.7. This is because as we increase δ, r max also increases and the system contains more pores with larger pore radii. This decreases the capillary barriers and hence reduces the global threshold pressure P c .

The P t vs β plane
To illustrate the all variations in the rheological behavior as a function of different parameters, we constructed a P t − β plane. This is shown in Fig. 6. We indicate all the data points there corresponding to the three sets. The variation in the β and P t as we increase S w , α and δ are indicated by arrows in the plot. They can be summarized as: Varying α: When α is increased, β decreases keeping P t constant, and the points moves to left along a horizontal line (green squares).
Varying δ: When δ is increased, P t decreases keeping β constant, and the points move down along a vertical line (red circles).
Varying S w : When S w deviates from 0.5, both β and P t decreases and the points move to lower values along a diagonal line (purple triangles).

Discussion
In this article we explored how the steady-state two-phase flow in porous media is affected by the underlying system disorder. We performed analytical calculations with a capillary bundle model and numerical simulations with dynamic pore-network model. By considering pore-radii distributions with a power α and a width δ we studied the transition from a non-linear flow regime to a linear Darcy flow. There we measured the exponent β related to the non-linear flow and cross-over pressure P t to the linear flow. We see that the exponent β is affected by α in the dynamic network model but not in the capillary fiber bundle model. The crossover point is affected by the width δ in both models. As |α| > 0 and we deviate from a uniform radii distribution, the slope in the non-linear regime decreases from 2 in the dynamic network model as it decreases the fluctuation in the pore-radii. On the other hand, when the width of the pore size distribution δ is increased, the linear region becomes achievable at relatively lower pressure gradient. When S w deviates from 0.5, the two-phase flow moves closer to a single phase flow. In such circumstances, both β and P t decreases simultaneously making the linear region more and more prominent. These numerical results can be explained qualitatively with the hypothesis that the total flow rate is contributed from the number of conducting paths and the flow in those paths [21], which makes it to increase faster than a linear behavior.
It is worth mentioning here that the single phase flow of Bingham fluid in porous media also shows similar non-linearity with a global yield threshold and a crossover to a linear regime [21,34,34,35]. Nash and Rees [36] has studied Binghum fluid flow in 1d and obtained different relations between applied pressure gradient and Darcy velocity depending on the distribution of channel widths. Talon et al. [37] has analytically explored the flow of Binghum fluid in 1d channels with aperture variation and observed different scaling between pressure gradient and flow rate depending on such variation in apertures. Above studies support the fact that, when up-scaled to a certain pore-network, both Binghum flow and two-phase flow is affected by the pore-size distribution and hence the topology of the of the network. The effect of pore size distribution has been explored extensively in case of a single phase flow where the distribution of local fluid velocity is observed to be affected by the distribution of pore-sizes [38][39][40][41][42]. In the present work we demonstrated how flow equations are affected by the pore-size distribution during a two-phase flow. It will be really interesting to extend the study of velocity distribution for the two phase flow and establish a link with the flow equations through the distribution of pores.

Acknowledgment
This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. SS was partially supported by the National Natural Science Foundation of China under grant number 11750110430.