Effective Rheology of Bi-Viscous Non-Newtonian Fluids in Porous Media

TL;DR

This paper investigates the flow behavior of bi-viscous non-Newtonian fluids in porous media, identifying a critical transition point characterized by specific critical exponents that depend on model parameters.

Contribution

It introduces a numerical model for bi-viscous non-Newtonian fluid flow in porous media and characterizes the critical transition with associated exponents.

Findings

Identification of a critical flow transition point.

Measurement of critical exponents related to the transition.

Dependence of exponents on model parameters.

Abstract

We model the flow of a bi-viscous non-Newtonian fluid in a porous medium by a square lattice where the links obey a piece-wise linear constitutive equation. We find numerically that the flow regime where the network transitions from all links behaving according to the first linear part of the constitutive equation to all links behaving according to the second linear part of the constitutive equation, is characterized by a critical point. We measure two critical exponents associated with this critical point, one of the being the correlation length exponent. We find that both critical exponents depend on the parameters of the model.

Equations

$$q(\nabla p)=\left\{\begin{array}[]{ll}-\alpha\nabla p&:|q|\leq q_{c}\;,\\ -\beta\nabla p+{\rm sgn}(q)\left[1-\frac{\beta}{\alpha}\right]q_{c}&:q_{c}\leq|q|\;,\\ \end{array}\right.$$

$$Q(\Delta P,\alpha ,\beta ,\{q_c\},N_x,{\lambda }_y{N}_y)$$

$$Q(\Delta P,\alpha ,\beta ,\{q_c\},N_x,N_y)$$

$$Q(\Delta P,{\lambda }_q\alpha ,{\lambda }_q\beta ,\{{\lambda }_q{q}_c\},N_x,N_y)$$

$$Q\left(\Delta P,\alpha ,\beta ,\{q_c\},N_x,N_y\right)$$

$$Q(\Delta P,\alpha ,\beta ,\{q_c\},N_x,N_y)$$

$$⟨q⟩(\nabla p,\alpha ,\beta ,\{q_c\})=Q(\nabla p,\alpha ,\beta ,\{q_c\},1,1).$$

$$Q(\Delta P,\alpha ,\beta ,\{q_c\},N_x,N_y)$$

$$⟨q⟩(\nabla P,\alpha ,\beta ,\{q_c\})$$

$$A=-\left(\frac{\partial ⟨q⟩}{\partial \alpha }\right),$$

$$B=-\left(\frac{\partial ⟨q⟩}{\partial \beta }\right),$$

$$\{c\}=\left\{\left(\frac{\partial ⟨q⟩}{\partial q_c}\right)\right\}.$$

$$A=A\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right),$$

$$B=B\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right),$$

$$\{c\}=\left\{c\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)\right\}.$$

$$⟨q⟩(\nabla P,\alpha ,\beta ,\{q_c\})$$

$$-A\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)\alpha$$

$$\left(\frac{\partial Q}{\partial \Delta P}\right)\Delta P=\left(\frac{\partial Q}{\partial \alpha }\right)\alpha +\left(\frac{\partial Q}{\partial \beta }\right)\beta ,$$

$$m\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)\nabla p$$

$$m=-\left(\frac{\partial ⟨q⟩}{\partial \nabla p}\right).$$

$$⟨q⟩=\left(\frac{\partial ⟨q⟩}{\partial \nabla p}\right)\nabla p+\sum _{\mathrm{links}}c{q}_c=-m\nabla p+\sum _{\mathrm{links}}c{q}_c,$$

$$⟨q⟩=$$

$$a\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)\nabla p=A\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right),$$

$$b\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)\nabla p=B\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right).$$

$$\frac{⟨q⟩}{\alpha }=\overline{q}\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right),$$

$$\overline{q}\left(\nabla p,\frac{\beta }{\alpha },\frac{\{q_c\}}{\alpha }\right)=$$

$$a\left(\nabla p,\frac{\beta }{\alpha }\right)=\Theta (q_c-|q|),$$

$$b\left(\nabla p,\frac{\beta }{\alpha }\right)=\Theta (|q|-q_c),$$

$$c\left(\nabla p,\frac{\beta }{\alpha }\right)=\Theta (|q|-q_c)\mathrm{sign}(q)\left(1-\frac{\beta }{\alpha }\right).$$

Full text

I Introduction

The behavior of complex fluids when being inside a porous medium may be very different from that when they are not. This is a problem encountered in many biological or industrial applications ranging from impregnation of fibrous materials to immiscible multi-phase flow in porous media. Among the different types of non-Newtonian fluids, many undergo behavioral changes depending on the stress or strain applied. One can mention the Carreau rheology which is Newtonian at low shear rate but behaves as a power law fluid above a certain shear rate c72 . Other examples are yield stress fluid that responds like a solid below a critical yield threshold. Above, the materials behave like a power law fluid hb26 . At the mesoscopic level, this rheological approach can also be extended to other situations. For example, inertial effects can be described as a rheological change from a Newtonian fluid to a power law (quadratic or cubic) for a certain large Reynolds number w96 . Another possible extension is the displacement immiscible fluids in porous media. In this case, the fluids may each be Newtonian. However, the interfacial tension between them, makes them effectively behave in a non-Newtonian way inside the porous medium sh12 . Indeed, a minimum amount of stress is then required for a non-wetting phase to invade small pore throats.

Non-Newtonian fluids are notoriously difficult to treat analytically and computationally. When in addition the flow is constrained by the very complex boundary conditions of the porous medium, the effective rheology of the fluid flow is not well understood. This might for example be seen in the fact that the leading theory for describing immiscible multi-phase flow in porous media is still the relative permeability theory dating from 1936 wb36 a theory which has evident faults.

The purpose of this manuscript is to investigate the coupling between the heterogeneities of the medium and a rheology with a change of behavior. We will study a very simple model called a bi-viscous fluid, where the fluid is Newtonian but with a change of viscosity at one particular shear rate (or shear stress) rhg87 ; hrh90 . The second viscosity might be lower (shear thinning) or higher (shear thickening). As we will see, the coupling between the disorder and such a simple rheological model is enough to generate a rich problem.

We also choose a simple porous medium, a square lattice oriented at 45∘ with respect to the average flow direction, see figure 1, consisting of $N_{x}$ links in the flow direction and $N_{y}$ links in the direction orthogonal to the flow direction.

The constitutive equation for the fluid in a link in the lattice is given by

[TABLE]

where $q$ is the volumetric flow rate in the link, $\nabla p$ is the pressure drop across the link. There are three parameters, $\alpha$, $\beta$ and $q_{c}$ The two first parameters, $\alpha$ and $\beta$ are mobilities when the fluid is either in the “$\alpha$-mode” or in the “$\beta$-mode.” The third parameter, $q_{c}$ is the flow rate at which the fluid changes from being $\alpha$-mode to $\beta$-mode. We illustrate the constitutive equation in figure 2. To simplify the problem as much as possible, we let the two mobilities $\alpha$ and $\beta$ be the same for all links in the lattice. However, each link has its own flow rate threshold $q_{c}$ drawn from a probability distribution $p(q_{c})$.

We will in the following study this system for arbitrary values of $\alpha$ and $\beta$ and for two threshold distributions; a uniform distribution and an exponential distribution.

In section II, we consider the symmetries inherent in the system. There are two types of symmetries. The first type is related to what happens to the volumetric flow rate through the system, $Q$ when we scale the parameters. Using the Euler theorem for homogeneous functions we are able to write down the most general form of the volumetric flow rate. If we define $\langle q\rangle$ as $Q/N_{y}$, we find that $\langle q\rangle=\alpha\ \overline{q}(\nabla p,\beta/\alpha,\{q_{c}\}/\alpha)$, where $\{q_{c}\}$ refers to the set of thresholds, one for each link. The second type of symmetry is the self-duality of the square lattice leading to a mapping between the behavior of the system for a given ratio $\beta/\alpha$ and its inverse, $\alpha/\beta$. Hence, we only need to discuss $\beta/\alpha\geq 1$.

We study in section III the lattice with $N_{x}=1$, i.e., there is only one layer. The model then becomes the capillary fiber bundle model which is analytically tractable. We find that for the uniform threshold distribution, the flow rate behaves as $\langle q\rangle-\langle q_{c}\rangle\sim(\nabla p-\nabla p_{c})^{2}$ where $(\langle q_{c}\rangle,\nabla p_{c})$ is a point only dependent on the value of the ratio $\beta/\alpha$ and the limits of the uniform distribution $q_{\min}$ and $q_{\max}$. This is reminiscent of a critical point. However, it is not a critical point. There are no correlations developing in the system as $\nabla p$ approaches $\nabla p_{c}$. Furthermore, the power law behavior is not seen when the threshold distribution is exponential.

Section IV is devoted to the numerical algorithm we use to solve the flow patterns. Our algorithm is based on the augmented Lagrangian algorithm, which we describe in this section.

We present our results in section V. First we note that the two limits $\beta/\alpha\to 1$ and $\beta/\alpha\to\infty$, or equivalently, $\beta/\alpha\to 0$ correspond to the directed percolation h00 and the directed polymer problems respectively hhr91 . This points us in the direction of there being a critical point in the problem in spite of the conclusion drawn for the capillary fiber bundle model in section III. Indeed, this is what we find: We find that $\langle q\rangle-\langle q_{c}\rangle\sim(\nabla p-\nabla p_{c})^{\mu}$ where $\mu$ depends on the ratio $\beta/\alpha$ for the same type of treshold distribution that gave a power law dependence in the capillary fiber bundle model studied in section III. We define and measure a correlation length $L_{\max}\sim(\nabla p-\nabla p_{c})^{-\nu}$. The correlation length exponent $\nu$ also depends on the ratio $\beta/\alpha$. In the limit $\beta/\alpha\to 1$, the longitudinal directed percolation correlation length exponent $\nu_{\parallel}=1.733847(6)$ j99 is expected and our numerical results are consistent with this. In the directed polymer limit $\beta/\alpha\to\infty$, however, the corresponding correlation length exponent is not the usual one, $\nu_{\parallel}=3/2$ rhh91 , but rather one that describes a correlated directed percolation problem.

The last section VI contains our summary and conclusions.

II Symmetries

In this section, we discuss the symmetries that lie hidden in the system we study, a square lattice of links obeying the constitutive equation (1). We consider two types of symmetry: one is based on scaling of the size and parameters of the model. Through the Euler theorem for homogeneous functions, we are able to write down the most general functional form the volumetric flow rate through the network takes. We then go on to exploring the geometrical symmetry inherent in the square lattice due to self duality in the same way as first done by Straley s77 . This symmetry demonstrates that we only need to explore the part of parameter space for which $\beta/\alpha\geq 1$.

II.1 Scaling symmetry

The volumetric flow rate $Q$ shows a number of scaling symmetries. We now combine these with the Euler theorem for homogeneous functions to deduce the functional form of $Q=Q(\Delta P,\alpha,\beta,\{q_{c}\},N_{x},N_{y})$ hsbkgv18 . Here $\{q_{c}\}$ is the set of thresholds, one for each link in the network. The volumetric flow rate is extensive in the width of the network, $N_{y}$. Hence,

[TABLE]

With respect to the length of the system, we find the symmetry

[TABLE]

A more subtle scaling symmetry is

[TABLE]

We also have the scaling symmetry

[TABLE]

The length $N_{x}$ and width $N_{y}$ of the network are discrete variables. By setting $\lambda_{y}=1/N_{y}$ we find from Equation (2) that

[TABLE]

The second scaling relation, Equation (3) gives when setting $\lambda_{x}=1/N_{x}$,

[TABLE]

where we have used the definition $\nabla p=\Delta P/N_{x}$. We now combine Equations (6) and (7) to get

[TABLE]

Hence, we define the average flow rate in the links as

[TABLE]

This is thus an intensive variable with respect to the width and the length of the network.

The two remaining scaling relations (4) and (5) involve continuous variables and we may thus make use of Euler’s theorem for homogeneous functions. The Euler theorem is easy to implement for each of these four scaling symmetries: we take the derivative with respect to the scaling variable $\lambda$ in each expression and set the variable equal to one.

The scaling relation (4) gives

[TABLE]

or in terms of the intensive variable

[TABLE]

We define the functions

[TABLE]

[TABLE]

and

[TABLE]

There is one function $c$ for each link in the network.

Whereas $\langle q\rangle$ is homogeneous of order one111A homogeneous function $f(x,y)$ of order $n$ in the variables $x$ and $y$ fulfills the scaling relation $\lambda^{n}f(x,y)=f(\lambda x,\lambda y)$. in the variables $\alpha$, $\beta$ and $\{q_{c}\}$, the functions $A$, $B$ and $\{c\}$ are homogeneous of order zero in these variables. This means that the parameters $\alpha$, $\beta$ and $\{q_{c}\}$ only appear as ratios in these functions,

[TABLE]

[TABLE]

and

[TABLE]

Equation (10) may thus be written

[TABLE]

Scaling equation (5) combined with the Euler theorem gives

[TABLE]

In terms of $\langle q\rangle$ and equation (17), we may rewrite this equation

[TABLE]

where we have defined the mobility

[TABLE]

From equations (10) and (19), we deduce that

[TABLE]

and with the help of equation (20) we find

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

We may take equation (23) one step further by dividing out the parameter $\alpha$,

[TABLE]

where

[TABLE]

We may as a check, compare equation (23) — our main result in this section — with the constitutive equation (1) in the case when there is no disorder, i.e., when all $q_{c}$ are equal. In this case, $\langle q\rangle$ should be equal to the constitutive equation. Hence, in this case we find,

[TABLE]

[TABLE]

and

[TABLE]

Here $\Theta$ is the Heaviside step function which is one for positive arguments and zero for negative arguments. We note that if $|q|<q_{\min}$, then equations (28) to (30) are correct as the disorder is not “noticeable” in this flow regime.

II.2 Self-duality of the square lattice

We define a dual network as sketched in Fig. 4. A node is located at the center of each cell and there is a link connecting each adjacent cell. On each link, a ”dual” current is defined from the pressure difference between pressure by the crossed link (from the original network),

[TABLE]

The current in the dual lattice satisfies the conservation of mass at each node (e.g. Kirchhoff condition) since: $j_{A\rightarrow B}+j_{A\rightarrow D}-j_{F\rightarrow A}-j_{E\rightarrow A}=0$

Moreover, one can define a pressure field $W$ on the dual lattice defined from this gradient,

[TABLE]

The definition is consistent once W is defined at a single point since the sum over a closed loop (and thus any) is equal to zero: $(W_{A}-W_{B})+(W_{B}-W_{C})+(W_{C}-W_{D})+(W_{D}-W_{A})=q_{1\rightarrow 4}+q_{1\rightarrow 5}+q_{1\rightarrow 6}+q_{1\rightarrow 2}=0.$

The “dual” pressure gradient and current follows thus the constitutive equation:

[TABLE]

Thus,

[TABLE]

The dual pressure and flow rate field satisfy thus the same kind of equation but with a local law which is inverted. It is important to note that the mean flow in the dual lattice is perpendicular to the original one.

III Capillary fiber bundle model

We now consider an analytically solvable model for the flow. Let us assume that the network consists of a set of parallel links placed between two fluid reservoirs kept at pressure $p=0$ and $p=\nabla p<0$, i.e., we are describing the capillary fiber bundle model s53 ; s74 ; rhs19 . The constitutive equation for the fiber bundle is then given by

[TABLE]

where we have labeled the links according to their position, $i=1,N_{y}$ and $q_{i}$ is the threshold of the $i$th link.

Let us now relabel the links in ascending order with respect to their thresholds: $q_{(1)}\leq q_{(2)}\leq\cdots\leq q_{(N_{y})}$. Equation (35) then becomes

[TABLE]

The thresholds are distributed according to the probability distribution $p(q_{c})$, with a corresponding cumulative probability given by

[TABLE]

According to order statistics, the mean value of $k$th largest threshold — mean value in the sense of averaging over an ensemble of networks — is given by

[TABLE]

Thus, the ensemble averages of the three types of sums in equation (36) are then

[TABLE]

[TABLE]

and

[TABLE]

Inserted into equation (36), these averages give

[TABLE]

where $\langle q\rangle=Q/N_{y}$.

III.1 Uniform threshold distribution

We now consider the concrete threshold distribution we will also employ in our numerical simulations on the square lattice: a uniform distribution on the interval $(q_{\min},q_{\max})$. Hence,

[TABLE]

We define

[TABLE]

and

[TABLE]

We also define

[TABLE]

Inserting these expressions into equation (42) gives

[TABLE]

We have here defined

[TABLE]

If we now define

[TABLE]

we may cast the middle regime where $|\nabla p_{\min}|<|\nabla p|<|\nabla p_{\max}|$ in the form

[TABLE]

It straight forward but somewhat tedious to rewrite the average flow rate $\langle q\rangle$, equation (47) in the general form (26) and (II.1) resulting from the scaling relations (2) to (5).

III.2 Exponential threshold distribution

Let us now consider the exponential threshold distribution

[TABLE]

for $0\leq q_{c}<\infty$. The corresponding cumulative distribution is

[TABLE]

Inserted into equation (42), this gives

[TABLE]

where we are still assuming $\nabla p<0$. Let us set $q_{0}=-\alpha\nabla p$. We then have the limits

[TABLE]

In contrast to the uniform distribution discussed in section III.1, there is not a transitional regime between the two limits of equation (54) which is on the form (50).

Hence, the uniform distribution on an interval, (43) results in $\langle q\rangle$ following a power law in $\langle q\rangle-\langle q_{c}\rangle$ vs. $\nabla p-\nabla p_{c}$, equation (50), whereas the exponential distribution (51) does not. From the simple capillary fiber bundle model we may conclude that the power law behavior seen in equation (50) is incidental and due to the uniform threshold distribution, which in itself is a power law.

We study a two-dimensional network mode in section V. Surprisingly, we find that also in this case, only the uniform distribution leads to a flow dependency on the pressure drop of the form

[TABLE]

In this case, however, the exponent $\mu$ depends on the parameter ratio $\beta/\alpha$.

IV Numerical method: Augmented Lagrangian

For completeness, this section describes the numerical method used to solve the non-linear Kirchhoff equation. This section is not required to understand the results that follow.

The method used is based on the Augmented Lagrangian method commonly used to solve the Stokes equation for yield stress fluids dl76 ; gl89 . It is based on a variational method of the problem. Indeed, if we rewrite the local equation (1) and introduce the function $f(q)$ as :

[TABLE]

We define the function $\phi(q)=\int_{0}^{q}f(q^{\prime})dq^{\prime}$. The flow field $\{q_{l}\}$ solution of equation (1), with the constraints of imposed inlet and outlet pressures at the boundaries $p_{in}$ and $p_{out}$, can be written as the saddle point of the functional

[TABLE]

where $\mathcal{L}$ represents the ensemble of links, $\mathcal{N}$ the ensemble of nodes and $\mathcal{V}(n)$ the ensemble of links connected to the node $n$. The symbol $\delta_{l,in}$ (resp. $\delta_{l,out}$) is equal to $1$ if the link is connected to the inlet (resp. outlet) node and to [math] otherwise. The $\{\lambda_{n}\}$ field is a set of Lagrangians which impose the conservation of mass at each node (and it can thus be associated to a pressure field).

The main idea of the Augmented Lagrangian method is to introduce a secondary set of velocities $\{j_{l}\}$ to decouple the nonlinear rheology to the Kirchhoff equation. Another constrain is then added $\{j_{l}\}=\{q_{l}\}$ via the Lagrangian method.

The velocity field is thus the solution of

[TABLE]

where $\{\mu_{l}\}$ is a Lagrangian set. The quadratic term is an additional penalty term which characterizes the augmented Lagrangian approach, where $\epsilon$ is a parameter.

The methods consists now in implementing an iterative algorithm to reach the saddle point starting from an initial guess $\{q^{0}_{l}\}$, $\{j^{0}_{l}\}$, $\{\lambda^{0}_{n}\}$ and $\{\mu^{0}_{l}\}$.

Knowing $\{q^{n}_{l}\}$, $\{j^{n}_{l}\}$, $\{\lambda^{n}_{n}\}$ and $\{\mu^{n}_{n}\}$, the algorithm is decomposed in the following steps:

Determination of $\{j^{n+1}_{l}\}$ and $\{\lambda^{n+1}_{n}\}$:

For this, we should solve:

[TABLE]

which reads

[TABLE]

where $\lambda^{n+1}_{l+}$ and $\lambda^{n+1}_{l-}$ are the lagrangian of the two nodes adjacent to the link $l$. For nodes adjacent to the outlet (resp. inlet), $\lambda_{+}$ (resp. $\lambda_{-}$) has to be replaced with $p_{out}$ (resp. $p_{in}$).

The most important point of this set of equations is that it is equivalent to solving the standard linear Kirchhoff equation with a constant permeability $1/\epsilon$ but with an additional source term $\mu_{l}^{n}-\epsilon q_{l}^{n}$. It can thus be solved by standard linear methods (Cholesky, LU decomposition, etc.).

Determination of $q_{l}^{n+1}$:

We solve

[TABLE]

yielding to the local, but implicit equation:

[TABLE]

Noting $y=\mu+\epsilon j^{n+1}_{l}$, the solution is:

[TABLE]

Determination of $\mu_{l}^{n+1}:$

For this we update in the direction of the gradient (Newton method):

[TABLE]

where $\gamma$ is a parameter set to $\gamma=\epsilon$ for simplicity.

In practice, this algorithm iterated until the relative variation of the total flow rate between two step is below $10^{-5}\%$. The computational time and the number of steps is strongly varying depending on $\beta$ but also on the applied pressure.

V Results

We now our numerical model based on the network show in figure 1 and the algorithm described in section IV. We use the link threshold distribution (43) with $q_{\min}=7.5$ and $q_{\max}=12.5$ in the following.

V.1 Criticality

As noted above, due to the distribution of thresholds, the links will reach their threshold at different macroscopic pressures. A link $l$ will be defined as being in $\beta$-mode if $q_{l}>q_{c}$ and in $\alpha$-mode otherwise. Similar to the percolation problem, a macroscopic change in flow regime is expected once the existence of percolation pathways of $\beta$-mode links. However, it is important to note a major difference with the percolation problem which lies in the fact that the mode of a link influences the neighboring links. Indeed, in the case of $\beta>\alpha$, once a link switches to $\beta$-mode, the flow will be easier through it. This will tend to concentrate the flow towards it. It will therefore increase the flow in the upstream and downstream neighboring links and therefore pushes these links towards the $\beta$-mode. In the opposite case, for $\beta<\alpha$, the $\beta$-mode have a lower conductivity once entering this mode compared to what it would have in $\alpha$-mode. Flow will therefore tend to go around it, increasing the flow in the lateral links. Consequently, $\beta$-mode links will tend to correlate in the streamwise (or lateral) direction for $\beta>\alpha$ and orthogonally to the streamwise direction for $\beta<\alpha$) wbhh96 .

The intermediate case $\beta=\alpha$ is interesting as the mode of a link has no influence on its neighbors. Since the mobility are the same for every link, the flow rate and the pressure gradient become homogeneous and equal to the mean flow rate and mean gradient. The problem is therefore identical to the directed percolation problem h00 .

The other limit $\beta/\alpha\gg 1$, the problem becomes identical to a yield stress fluid hhr91 ; grhbtc90 ; tb13 . The critical path is then related to the directed polymer problem hhr91 ; kz87 ; taph13 , as it corresponds to the path that minimize the sum of local pressure threshold $\Delta P_{c}=\min\sum(q_{c}/\alpha)$.

V.2 Pathscape method

To quantify this phenomenon and to determine the percolation pressure, we want to determine the longest directed path of the $\beta$-mode links. This quantity is essentially the longitudinal correlation length in directed percolation j99 . We map the length of all paths by invoking a pathscape approach as described in taph13 for yield-stress fluids.

We introduce the node field $L_{n}$ representing the longest upstream directed path ending at $n$. $L_{n}$ can be determined from a transfer matrix algorithm propagating from left to right (stream direction). If we note, at a given node $n$, $l_{1}$ and $l_{2}$ the two upstream neighbor links and $n_{1}$ and $n_{2}$ the corresponding nodes. We associate binary variables $m_{1}$ and $m_{2}$ with the two links $l_{1}$ and $l_{2}$. If link $l_{1}$ is in $\beta$-mode, then $m_{1}=1$, otherwise $m_{1}=0$ — and likewise for link $l_{2}$. We then have that

[TABLE]

We then construct the node field $R_{n}$ containing the longest directed path ending at $n$ but propagating in the downstream direction. The algorithm is identical to the previous one but it propagates in the upstream direction from the rightmost column.

Once both fields have been determined, we sum the two to obtain the pathscape $T_{n}=L_{n}+R_{n}$, which contains the length of longest directed percolating path passing by the node $n$. From this pathscape, we can then identify the longest directed path $L_{max}=\max(T_{n})$ In figure 5, we present two examples of such a pathscape at two different imposed pressure. In this figure, we can see the longest cluster path in dark blue. At low applied pressure, the longest cluster is quite low $L_{\max}=7$, whereas at higher pressure, $L_{\max}$ is closer to the system size.

It is important to note that the pathscape we have defined here is not the landscape of minimal paths taph13 . In the limit $\beta\to\alpha$ the pathscape reflects the clusters in directed percolation as noted in section V.1. However, when $\beta\neq\alpha$, the paths we identify correspond to directed percolation clusters. However, the directed percolation is now correlated.

V.3 Evolution of the correlation length $L_{\max}$

In figure 6, we investigate the evolution of $L_{\max}$ as function of the applied pressure. As it can be seen, the correlation length increases with pressure until it reaches the system length $N_{x}$. Similarly to percolation, one can see in figure 6(b) that the correlation length diverges as a power law close to a critical pressure gradient $\nabla p_{c}$,

[TABLE]

We note in this figure that the exponent $\nu$ seems to vary with $\beta$. In figure 7, we display the evolution of $\nu$ and the critical pressure gradient $\nabla p_{c}$ against the parameter $\beta$. As we can see, $\nu$ and $\nabla p_{c}$ decrease significantly with $\beta$. Where the limit $\beta\rightarrow 1$ is consistent with the results found in the literature on directed percolation, $\nu=\nu_{\parallel}=1.733847(6)$ j99 . Our best estimate of the threshold pressure is $\nabla p_{c}\approx 10.72$.

At the end of section V.2 we noted that the pathscape we have identified is not related to the pathscape spanned by minimal paths in the limit $\beta/\alpha\to\infty$. If that were the case, we would have expect $\nu$ to approach the value $\nu_{\parallel}=3/2$ rhh91 . Rather, we are identifying directed percolation clusters in a correlated landscape, and this directed percolation $\nu$ is approaching the value 1 in this limit.

V.4 Flow curve

We now investigate the flow curve. Figure 8 displays the evolution of the mean flow rate as function of the pressure gradient and for different $\beta$. In the lower figure, we show that, close to the critical pressure, the flow rate follows also a power-law which can be written in the form:

[TABLE]

where $q_{c}$ is a constant obtained by interpolating the data at the critical pressure. Here also we remark that the exponent $\mu$ varies with the the coefficient $\beta$. In figure 9, we report the evolution of this exponent as a function of $1/\log(\beta)$. For $\beta=\alpha=1$ we have the obvious limiting value $\mu=1$. As $\beta$ increases, so does the value of $\mu$. By plotting $\mu$ against $1/\log(\beta)$ we may estimate the limiting value for $\beta\to\infty$, which is not inconsistent with the value $\mu=2$; the value suggested by Roux and Herrmann in 1987 rh87 .

We note that the functional form $\langle q\rangle$, equation (68), based on the uniform threshold distribution (43), gives a behavior closely related to the one found for the capillary fiber bundle model with the same type of threshold distribution, see equation (50), but with $\mu=2$. The correlation length exponent $\nu$ cannot be defined in the capillary fiber bundle model.

In section III.2, we studied the capillary fiber bundle model with an exponential threshold distribution (51). We have used the same distribution for the network model considered here. As in the capillary fiber bundle model, we do not find a power law of the type (68) in this case, nor do we find a power law for the correlation length, (67).

VI Summary and Conclusions

We have in this paper explored the behavior of a bi-viscous fluid moving in a square lattice subject to the constitutive equation (1) for each link. This system contains a critical point which leads to the behavior $\langle q\rangle-\langle q_{c}\rangle\sim(\nabla p-\nabla p_{c})^{\mu}$ for the volumetric flow rate and $L_{\max}\sim(\nabla p-\nabla p_{c})^{-\nu}$ for the correlation length when a uniform threshold distribution is used. However, the two limits of the ratio between the two parameters representing the mobilities, $\beta/\alpha\to 1$ and $\beta/\alpha\to\infty$, or equivalently, $\beta/\alpha\to 0$ correspond to the percolation and the directed polymer problems respectively. These are problems containing critical points.

There are still a number of open questions concerning this system. We list them as follows:

• •

We have only considered $\nabla p\geq\nabla p_{c}$. What happens on the other side of the critical point?

• •

The critical exponents $\mu$ and $\nu$ are functions of the parameter ratio $\beta/\alpha$. Is this a crossover or are we dealing with a non-universal exponent?

• •

We have only dealt with $\beta\geq 0$. What happens for $\beta<0$? The limit $\beta\to-\infty$ turns the model into the fuse model. What happens when $\beta$ is barely negative? Our numerical algorithms is not capable of handling this problem.

• •

It would be more realistic, but also more challenging to consider a power-law type characteristic for the constitutive equation for $q\geq q_{c}$. How will this change our conclusions?

• •

Why do we not see critical behavior for the exponential threshold distribution in the network model?

We thank the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. AH thanks the Université de Paris-Sud for funding through a visiting professorship.