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

TL;DR

This study investigates the effective rheology of immiscible two-phase flow in porous media with mixed wetting properties, revealing regimes with zero threshold pressure and power-law flow behaviors through a dynamic pore network model.

Contribution

It introduces a detailed analysis of flow regimes in porous media with mixed wettability, highlighting the existence of zero-threshold flow paths and quantifying flow exponents near percolation thresholds.

Findings

Percolation regime with zero capillary threshold pressure.

Flow rate transitions from linear to power-law and back to linear with pressure.

Mobility scales with occupation probability near critical point with exponent ~5.7.

Abstract

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.

Equations

$$Q=-\frac{AK}{\mu L}\Delta P,$$

$$Q_w=-\frac{AK{k}_{rw}}{{\mu }_{w}L}\Delta P,$$

$$Q_n=-\frac{AK{k}_{rn}}{{\mu }_{n}L}\Delta P,$$

$$Q=Q_w+Q_n,$$

$$Q=-\frac{AK}{L}\left[\frac{k_{rw}}{{\mu }_w}+\frac{k_{rn}}{{\mu }_n}\right]\Delta P,$$

$$Q\propto |\Delta P{|}^{\beta },$$

$$k_{rn}\sim (P_c-P_c^c{)}^t,$$

$$Q=\ \left\{\begin{array}[]{ll}0&\mbox{if $|\Delta P|\leq P_{t}$}\;,\\ M(|\Delta P|-P_{t})^{\beta}&\mbox{if $P_{t}<|\Delta P|<P_{\max}$}\;,\\ M_{D}(|\Delta P|-P_{t})&\mbox{if $P_{\max}\ll|\Delta P|$}\;,\end{array}\right.$$

$$\beta=\ \left\{\begin{array}[]{ll}1&\mbox{if $|\Delta P|-P_{t}\gg P_{\max}$}\;,\\ 2&\mbox{if $P_{t}\ll|\Delta P|-P_{t}\ll P_{\max}$}\;,\\ 3/2&\mbox{if $0<|\Delta P|-P_{t}\ll P_{t}$}\;.\end{array}\right.$$

$$q=-\frac{\pi {\bar{r}}^4}{8\mu l}\left(\Delta p-\hat{x}\cdot \left(\sum _k{p}_t(x_k){\hat{J}}_k\right)\right),$$

$$p_t(x)=\frac{2\sigma \cos \theta }{r(x)},$$

$$r(x)=\frac{r_0}{1-a\cos \left(\frac{2\pi x}{l}\right)}$$

$$Q=\ \left\{\begin{array}[]{ll}M|\Delta P|&\mbox{if $|\Delta P|<P_{l}$}\;,\\ M_{m}|\Delta P|^{\beta_{3}}&\mbox{if $P_{l}<|\Delta P|<P_{u}$}\;,\\ M_{D}|\Delta P|&\mbox{if $P_{u}<|\Delta P|$}\;,\\ \end{array}\right.$$

$$M\sim\ \left\{\begin{array}[]{ll}(p^{+}-(1-p_{c}))^{t^{\prime}}&\mbox{for $p^{+}\to(1-p_{c})^{+}$}\;,\\ (p_{c}-p^{+})^{t^{\prime}}&\mbox{for $p^{+}\to(p_{c})^{-}$}\;,\\ \end{array}\right.$$

$$M\sim L^{-t^{\prime }/\nu }$$

$$t^{\prime }=5.7\pm 1.3.$$

$$Q=\ \left\{\begin{array}[]{ll}0&\mbox{if $|\Delta P|\leq P_{t}$}\;,\\ M_{m}(|\Delta P|-P_{t})^{\beta}&\mbox{if $P_{t}<|\Delta P|<P_{u}$}\;,\\ M_{D}(|\Delta P|-P_{t})&\mbox{if $|\Delta P|<P_{u}$}\;,\\ \end{array}\right.$$

Full text

I Introduction

It was in 1827 that Ohm published his law stating that electrical current is proportional to the voltage drop across a conductor Ohm (1827), 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 Darcy (1856). 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 Whitaker (1986).

The Darcy law for single-phase flow through a porous sample is

[TABLE]

where $Q$ is the volumetric flow rate along the axis of the cylindrical sample, $\Delta 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, $\mu$ 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 Wyckoff and Botset (1936),

[TABLE]

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. Barenblatt et al. (2002) 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,

[TABLE]

and as a consequence, the generalized Darcy equations (2) and (3) predict

[TABLE]

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 $\nabla P_{c}=(dP_{c}/dS_{w})\nabla 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. Tallakstad et al. (2009a, b) reported a flow rate $Q$ depending on $\Delta P$ as

[TABLE]

with $\beta\approx 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. Rassi et al. (2011) finding an exponent $\beta$ varying between 2.2 and 3.3. Aursjø et al. Aursjø et al. (2014) using the same model porous medium as Tallakstad et al. Tallakstad et al. (2009a, b), but with two incompressible fluids, found $\beta\approx 1.5$ or $1.35$ depending on the fractional flow rates. Similar results in the sense that $\beta$ is considerably larger than one, have since been observed by a number of groups, see Sinha et al. (2017); Gao et al. (2020); Zhang et al. (2021, 2022). There has also been a considerable effort to understand these results theoretically and reproduce them numerically Tallakstad et al. (2009a, b); Grøva and Hansen (2011); Sinha and Hansen (2012); Sinha et al. (2013); Xu and Wang (2014); Yiotis et al. (2019); Roy et al. (2019a, b); Fyhn et al. (2021); Sales et al. (2022); Feder et al. (2022); Lanza et al. (2022); Cheon et al. (2023).

It should be pointed out that the power law behavior seen in Eq. (6) is different from the one described by Wilkinson in 1986 Wilkinson (1986). 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

[TABLE]

where $t$ is the percolation conduction exponent Stauffer and Aharony (2018). 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 Roux and Herrmann (1987). Why it should be a power law is less obvious. The best argument for why was perhaps already given by Tallakstad et al. Tallakstad et al. (2009a, b) through comparing pressure drop across fluid clusters with the capillary pressures holding them in place. Capillary fiber bundle models Scheidegger (1953, 2020) are porous media in the form of bundles of capillary fibers and they are typically simple enough to be mathematically solvable Sinha et al. (2013); Roy et al. (2019a, b); Fyhn et al. (2021); Lanza et al. (2022); Cheon et al. (2023). 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

[TABLE]

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 non-zero 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 Sinha and Hansen (2012); Sinha et al. (2017); Fyhn et al. (2021). The existence of a non-zero threshold pressure makes the measurement of $\beta$ much harder than when it is zero as this implies determining two parameters simultaneously, $(P_{t},\beta)$, rather than just one, $\beta$.

A central unanswered question is whether the exponent $\beta$ 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 $|\Delta P|\to P_{t}$ from above Tallakstad et al. (2009a, b). The experimental measurements of $\beta$ 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 $\beta$. Roy et al. Roy et al. (2019a) found using a capillary fiber bundle model that $\beta=2$ if the fiber-to-fiber probability distribution of thresholds includes $P_{t}=0$ with a finite probability, otherwise $\beta=3/2$. The fibers here had smoothly undulating radii along the flow direction. Lanza et al. Lanza et al. (2022) who studied non-Newtonian a mixture of immiscible Newtonian and non-Newtonian fluids in a capillary fiber bundle model found a different value of $\beta$ 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 $|\Delta 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. Fyhn et al. (2021) studied the exponent $\beta$ 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

[TABLE]

The network model studies showed a less clear picture, with $\beta$ 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 $\beta$ 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 $\beta$.

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. Irannezhad et al. (2022, 2023). 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^{+}={\rm constant}$ critical lines. Each line signifies a percolation transition Stauffer and Aharony (2018). 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 Knudsen and Hansen (2006).

In the region where $P_{t}=0$, while still having two-phase flow, we observe an exponent $\beta=\beta_{3}=2.56\pm 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}\approx 0.5927$ is the site percolation threshold for square lattices, i.e., $\beta=\beta_{2}=\beta_{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 $\beta=\beta_{1}=2.25\pm 0.1$. The large uncertainty in $\beta$ seen here stems from $P_{t}>0$.

It is surprising that $\beta_{2}=\beta_{3}$ within the precision we are able to obtain. The exponent $\beta_{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 Grossman and Aharony (1986). The reason for not seeing critical behavior reflected in $\beta$ 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^{\prime}$ of the same type as the conductivity exponent $t$ percolation in ordinary percolation in the vicinity of the critical $p^{+}$ Stauffer and Aharony (2018); Redner (2007). We find that $t^{\prime}\approx 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 Sinha et al. (2021); Aker et al. (1998); Knudsen et al. (2002); Tørå et al. (2012); Gjennestad et al. (2018) for this study. It has been used earlier in the context of modeling mixed wetting porous media, see (Sinha et al., 2011; Flovik et al., 2015; Fyhn et al., 2021). We describe the model in Section II including our use of the wetting model similar to the one introduced by Irannezhad et al. Irannezhad et al. (2022, 2023). 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 $\Delta 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^{\circ}$ angle from the flow direction. The $\Delta 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 $\theta$ to the links, where $\theta$ is consistently defined through one of the fluids. In contrast to earlier models Sinha et al. (2011); Flovik et al. (2015); Fyhn et al. (2021) 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. Irannezhad et al. (2022, 2023). We assume two types of grains, they being either fully non-wetting with $\theta=180^{\circ}$ or fully wetting with $\theta=0^{\circ}$. 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, $\theta$ 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 $\theta=180^{\circ}$. If both of the adjacent grains are assigned “$-$”, then $\theta=0^{\circ}$. 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 $\theta=90^{\circ}$.

The networks have periodic boundary conditions in both directions. Two fluids that flow through the network are immiscible and their movement is traced through the position of their interfaces at each instant in time. Whenever the fluids flowing in a link reach the crossing point with the three other links, namely a node, the fluids get distributed into the neighboring links in the same time step instead of being retained in the node itself Sinha et al. (2021).

The volumetric flow in each link with length $l$, pointing along the link’s center-axis $\boldsymbol{x}$, is given by

[TABLE]

where it has been assumed that the radius does not deviate too much from its average value $\bar{r}$ Sinha et al. (2021). Here, $\mu=s_{A}\mu_{A}+s_{B}\mu_{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 $\mu_{A}$ and $\mu_{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 $\hat{J}$ lies along the $\boldsymbol{x}$ axis, and points in the direction out of the fluid within which $\theta$ is defined. In Fig. 3, $\theta$ 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, $\hat{J}$ 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 $\hat{J}_{k}$ with positions $x_{k}\in[0,l]$ along $\boldsymbol{x}$. The dot product of this sum with the unit vector $\hat{x}$ in the positive $\boldsymbol{x}$ direction is taken afterward to obtain the total capillary pressure. The capillary pressure across one interface at position $x$ which has an angle $\theta$ with the solid through fluid A is modeled by using the Young-Laplace equation Blunt (2017),

[TABLE]

where $\sigma$ is the surface tension and

[TABLE]

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\cdot 10^{-5}$ N$/$mm and viscosity $0.1$ Pa$\cdot$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 Stauffer and Aharony (2018). The links altogether form a square lattice. The nodes of the dual lattice, form another square lattice Straley (1977) 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^{+}\geq p_{c}$, where $p_{c}$ is the site percolation threshold $0.5927\dots$. If we, on the other hand, focus on the “$-$” sites, there will be a cluster of such sites that percolate if $p^{-}=1-p^{+}\geq p_{c}$, or $p^{+}\leq 1-p_{c}\approx 0.4073\dots$ Xun et al. (2021). Hence, if $0\leq p^{+}\leq 1-p_{c}$, the “$-$” clusters percolate, if $1-p_{c}\leq p^{+}\leq p_{c}$, neither the “$-$” sites nor the “$+$” sites percolate, and if $p_{c}\leq 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.

We note that if neither the “$+$” sites nor the “$-$” sites percolate ($1-p_{c}\leq p^{+}\leq 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$ Grossman and Aharony (1986). 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}\leq p^{+}\leq 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\times L$ for each $p^{+}\in\{0.3000,0.3001,0.3002,\dots,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 $|\Delta P|$ can be written as

[TABLE]

in the region $1-p_{c}\leq p^{+}\leq 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 Redner (2007), we make the assumption that the mobility vanishes as

[TABLE]

where $t^{\prime}$ is a transport exponent of the same type as the conductivity exponent $t$ in ordinary percolation, which is $1.303(8)$ according to Cen et al. (2012). In Eq. 14, $p^{+}\to(1-p_{c})^{+}$ means $p^{+}$ approaches $1-p_{c}$ from above and $p^{+}\to(p_{c})^{-}$ means $p^{+}$ approaches $p_{c}$ from below. By using finite size scaling analysis, we have that

[TABLE]

where $\nu$ is the correlation length exponent in percolation, which is known to be $4/3$ Den Nijs (1979).

To investigate the relation given in Eq. 15, we set $p^{+}=0.5927\approx p_{c}$ and network-dimensions $L\times L$ for $L$ between $50$ to $90$ links. The lowest numerically feasible $|\Delta P|$ are used in order to stay in the lower linear regime in Eq. 13, specifically $2.8$ Pa/link $\lessapprox|\Delta P|/L\lessapprox 5.8$ Pa/link. When operating at low $|\Delta 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 $|\Delta P|$. For these simulations, the flow is driven for approximately $40$ pore-volumes of fluid through the network, where one pore-volume 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^{\prime}/\nu=4.3\pm 1.0$, giving

[TABLE]

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\times 100$ links2. For each $|\Delta 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.

V.1 Hysteresis

We pose here the question of whether there are any hysteretic effects from raising and lowering the pressure drop $|\Delta 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, $|\Delta P|$ applied across a network is raised and then lowered in steps. The $|\Delta 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 $|\Delta P|$ is returned to the same value, $Q$ also quickly stabilizes back to the previous value it had with the same $|\Delta P|$. This shows that the steady state results generated using the DPN model do not depend on long-term memory Erpelding et al. (2013).

V.2 Volumetric flow rate dependence on pressure drop

The results relating $Q$ and $|\Delta P|$ in systems with zero $P_{t}$ and different values of $p^{+}$ are shown in Fig. 9. We used $p^{+}\in\{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 $|\Delta P|$ have been plotted in Fig. 10. The velocity maps show the steady state averaged absolute velocities, in other words, they show the average speed of the fluid. The velocities are color coded so that those through neutral links with $\theta=90^{\circ}$ are in shades of red and the rest that are through links with $\theta\in\{0^{\circ},180^{\circ}\}$ are in shades of blue. The results in Figs. 9 and 10 show three regimes in terms of $\beta$ as indicated in Eq. (13).

The lowest regime in Eq. 13 seems to correspond to $\log|\Delta P|\lessapprox 2.8$ in Fig. 9, in other words $|\Delta P|/L\lessapprox 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 $|\Delta P|$, we find $\beta=1.00\pm 0.01$. The velocity maps of a network with $p^{+}=0.5$ at two different $|\Delta 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|\Delta 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 $\beta$ 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}\approx 0.5927$ and $1-p_{c}\approx 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 $|\Delta P|$ hence resulting in a larger $Q$, meaning a larger $M$ when $p^{+}$ moves towards 0.5. For instance, at $\log|\Delta P|\approx 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\lessapprox\log|\Delta P|\lessapprox 4.1$ in Fig. 9, in other words $20.0$ Pa/Link $\lessapprox|\Delta P|/L\lessapprox 125.9$ Pa/link. Here, the exponent in Eq. 13 is $\beta=\beta_{3}=2.56\pm 0.05$ and $M_{m}$ is the same for all $p^{+}$ examined. When $\log|\Delta P|$ increases from Fig. 10(c) to Fig. 10(d) in this regime, the velocity maps show that there is a significant increase in the number of flow carrying links, meaning $Q$ increases significantly also. The opening of new paths in addition to increased flow in the already active paths explains $\beta$ being large. At this level of $|\Delta P|$, $M_{m}$ and $\beta$ being the same for all $p^{+}$ examined makes sense as the connected paths that differentiate networks with different $p^{+}$ no longer are the main contributors to the flow.

The highest regime in Eq. 13 seems to corresponds to $\log|\Delta P|\gtrapprox 4.5$ in Fig. 9, in other words $|\Delta P|/L\gtrapprox 316.2$ Pa/link. Here, the exponent in Eq. 13 is $\beta=1.00\pm 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 $|\Delta 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 $|\Delta P|$. Increasing $|\Delta 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}\leq p^{+}\leq p_{c}$ examined in Fig. 9. If $p^{+}$ were very close to the range examined in Fig. 9, the behavior of $\beta$ 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

[TABLE]

where $P_{u}$ is the crossover pressure between non-linear and Darcy behavior. By varying $P_{t}$ from $0.00$ Pa to the lowest $|\Delta P|$ in the data sets with an increment of $0.01$ Pa, mathematical linear fits with slopes $\beta$ were calculated at the lowest pressures to find the candidate that gave the least room-mean-square error. This gave $\beta=2.23\pm 0.05$ and $P_{t}=\left(3.4\pm 0.5\right)$ kPa for $p^{+}=0.2$ and $\beta=2.29\pm 0.05$ and $P_{t}=\left(2.0\pm 0.5\right)$ kPa for $p^{+}=0.3$. The regime these $\beta$ 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. (2021) has observed $\beta>1$ behavior even in networks with the same wetting angle everywhere which would be the same as having $p^{+}\to 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 $\beta\approx 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 immiscible 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. Irannezhad et al. (2022, 2023). 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}\leq p^{+}\leq 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^{+}\ll 1-p_{c}$ or $p^{+}\gg p_{c}$, easy links do not percolate. We find two classes of constitutive equations for volumetric flow rate $Q$ vs. pressure drop $|\Delta P|$: For $1-p_{c}\leq p^{+}\leq p_{c}$ we observed a constitutive equation as in Eq. (13), see Fig. 9, whereas for $p^{+}\ll 1-p_{c}$ or $p^{+}\gg 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}\leq p^{+}\leq p_{c}$ we observed the following: At the regimes with lowest and highest $|\Delta P|$, it seems that $\beta=1.00\pm 0.01$, because there is no significant change in the paths fluids are flowing through and increasing $|\Delta P|$ only increases the flow in the already active paths. At the lowest $|\Delta P|$, the flow is mainly always through connected paths with zero resistance. When $p^{+}\to 0.5$ in this regime, there are more connected paths which means more fluid gets transported, making $Q$ hence $M$ higher. At the highest $|\Delta P|$, almost the entire network is always active. On the other hand, $\beta>1$ in the middle regime where an increase in $|\Delta 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 $\beta=\beta_{3}=2.56\pm 0.05$. We saw no systematic dependence of $\beta$ on $p^{+}$.

For $p^{+}=0.2$, however, $\beta=2.23\pm 0.05$ and $P_{t}=\left(3.4\pm 0.5\right)$ kPa, and for $p^{+}=0.3$, $\beta=2.29\pm 0.05$ and $P_{t}=\left(2.0\pm 0.5\right)$ kPa. Due to the necessity of determining $P_{t}$ and $\beta$ simultaneously at these $p^{+}$ values, there is more uncertainty associated with the measurements of $\beta$. It is not possible to verify or falsify whether there is a fixed $\beta=\beta_{2}$, or whether it depends on $p^{+}$ and $S_{w}$.

The existence of connecting paths is a percolation problem. They disappear when $p^{+}\to(1-p_{c})^{+}$ or $p^{+}\to(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^{\prime}\approx 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.

Acknowledgment: We thank S. B. Santra and D. Jnana for discussions on the percolation problem that this system poses.