# Non-Newtonian flow in macroscopic heterogeneous porous media: from power-law
fluids to transitional rheology

Laurent Talon

###### Abstract

Many natural or industrial fluids exhibit non-linear behavior. Non-Newtonian
fluids can therefore be found in different applications related to porous or
fractured media such as mud flows, oil recovery, hydraulic fracturing or foam
injection clean-up. This work focuses on the understanding of the flow of
these fluids in large scale heterogeneous porous media, and particularly in
the coupling between flow heterogeneity and rheology. This work consider
fluids that exhibit a change in behavior, where the viscosity changes above a
certain velocity threshold. This change of viscosity affects drastically the
inhomogeneity of velocity field. In this paper, we analyse some general
mathematical properties for non-linear rheologies at the microscopic and
macroscopic scale. Most non-linear rheologies imply the existance of two
asymptotic limit, at low and high flow rate. The velocity field of both limits
are obtained by solving a power-law rheology in the medium, which is analysed
using a perturbation expansion. As an example, we will consider fluids which
are Newtonian for low shear-rate but becomes shear-thinning or shear-
thickening at higher shear rate. Three macroscopic flowing regimes are then
observed. It will be shown that the limits at low and high flow-rate are
governed by the properties of power-law fluids. The transition between these
two regimes is the analysed. The range depends on the heterogeneity
permeability field but also on the rheological parameters. The statistical
properties of the flow field displays interesting statistical and geometrical
properties such as multi-scale (fractal), similarly to other critical systems
(percolation, avalanches, etc.)

##  1 Introduction

Many natural or industrial fluids exhibit non-Newtonian behaviors, (Bird et
al., 1987; Barnes et al., 1989; Coussot, 2005), they are thus found in many
applications related to porous or fractured media. An important industrial
application is, for example, enhanced oil recovery (EOR) (see Sorbie et al.
(1989); Sorbie (2013)). Indeed, oil is usually recovered by displacing it with
another fluid (e.g. water). The main problem lies in the fact that the
displacing fluid is often less viscous, resulting in viscous fingering. The
fingering tends to create preferential flow paths, leaving many areas of oil
inplace. The idea is then to inject a non-Newtonian fluid in order to prevent
fingering. Another interesting application is the description of blood in the
capillary network, which should be regarded as a suspension and thus with a
non-Newtonian viscosity: shear-thinning or yield stress (Boyd et al., 2007;
Bessonov et al., 2016). Non-Newtonian fluid is also present in the proppants
used for hydraulic fracturing or the mud produced by drilling wells. Non-
Newtonian fluids (cements, polymers, etc.) are also commonly used for fracture
sealing (Tongwa et al., 2013).

A very recurring problem when dealing with porous media is that of upscaling.
If the equations of motion are generally well known at the pore scale (_e.g_
typically $\sim 10^{-3}\;\rm{m}$), a particular interest is in understanding
the flow at much larger scales ($\sim 1-10^{3}\;\rm{m}$). This is usually done
by deriving constitutive equations for average quantities at an intermediate
scale. This is illustrated by the famous Darcy’s law for Newtonian fluids,
which relates linearly the mean gradient of pressure to the mean flow rate. At
the microscopic level, fluids obey the Stokes equation (if inertia is
neglected) :

$\vec{0}=-\vec{\nabla}P+\mu\Delta\vec{u},$ (1)

where $\vec{u}$ is the fluid velocity, $P$ is the pressure and $\mu$ is the
viscosity. After averaging the velocity and pressure field over a large number
of pores, it results Darcy’s law (Darcy, 1856) :

$\langle\vec{u}\rangle=-\frac{\kappa}{\mu}\langle\vec{\nabla}P\rangle,$ (2)

where $\kappa$ is the permeability of the porous medium. $\kappa$ depends on
its topology and is homogeneous to a squared length. In the following, for the
sake of simplicity, the $\langle.\rangle$ is usually omitted when the scale of
interest is explicit. This scale, containing enough pores for Darcy’s law to
be applicable, is called mesoscopic.

At the gelogical scale the type of material may however spatially vary leading
to a macroscopic heterogeneous permeability field. To understand the flow at
large scale then requires solving the heterogeneous Darcy’s law :

$\vec{u}=-\frac{\kappa(\vec{r})}{\mu}\vec{\nabla}P\;\;\;{\rm
and}\;\;\;\vec{\nabla}.\vec{u}=0.$ (3)

It is also important to recall that a very similar equation is used for
solving flow in rough fractures, usualy refered as the Reynolds equation
(Reynolds, 1886; Zimmerman et al., 1991; Mourzenko et al., 1995). Indeed, in a
fracture with varying opening and under the lubrication approximation (small
thickness and small variation of the opening), flow obeys :

$\vec{q}(x,y)=-\frac{\kappa(x,y)}{\mu}\vec{\nabla}_{2D}P\;\;{\rm
et}\;\vec{\nabla}_{2D}.\vec{q}=0,$ (4)

where the fracture is in the $(x,y)$-plane, $b(x,y)$ represents the local
opening, $\vec{q}(x,y)=\int_{0}^{b(x,y)}\vec{u}\;dz$ and
$\kappa(x,y)=\frac{b^{3}(x,y)}{12}$. The fracture can thus be treated as a 2D
porous medium. There is, however, a small caveat because the dimensions are a
slightly different. In a porous medium, Darcy’s law implies an average
velocity ($m.s^{-1}$) and permeability has the dimension of $m^{2}$, whereas
in a fracture, $\vec{q}$ is a flow per unit length ($m^{2}.s^{-1}$) and the
”permeability” has the dimension of $m^{3}$. We will see later, this
difference has some consequences in the model equations.

All the above mentioned equations apply to Newtonian fluids. A question that
naturally arises is: how should this approach be modified when considering
non-Newtonian fluids? While there is a very large variety of non-Newtonian
fluids (Bird, 1976; Bird et al., 1987; Coussot, 2005), we will limit our study
to fluids whose viscosity varies with shear rate, which are the most commonly
non-Newtonian fluids. For these fluids, the rheology is characterized by the
relationship between shear rate $\dot{\gamma}$ and viscosity
$\mu(\dot{\gamma})$, or equivalently between $\dot{\gamma}$ and shear stress
$\tau(\dot{\gamma})$, with $\tau(\dot{\gamma})=\mu(\dot{\gamma})\dot{\gamma}$.

Due to the importance of the subject, a large number of experimental,
numerical and theoretical studies have been carried out on the determination
of a Darcy’s law for non-Newtonian fluids. The most classical approach (see
for instance (Christopher & Middleman, 1965; Sadowski & Bird, 1965; Slattery,
1967; Hirasaki et al., 1974; Chauveteau, 1982; Sorbie et al., 1989; Pearson &
Tardy, 2002; Lopez et al., 2003; Sochi & Blunt, 2008; Balhoff & Thompson,
2004; Balhoff, 2005)) consists in determining an effective shear rate
$\dot{\gamma}_{pm}$ in order to derive an effective viscosity. Another
approach is to determine an effective stress (McKinley et al., 1966), or an
average viscosity (Eberhard et al., 2019). A common feature of these
approaches is that they are based on the determination of effective
quantities. They can be synthesized using scaling arguments. Indeed, by
defining a typical length scale $\lambda$ (pore size, grain diameter,
$\sqrt{\kappa}$, etc.) and using the average flow rate $u$, a typical shear
rate $\dot{\gamma}_{eff}=\alpha u/\lambda$ can be defined. Similarly, using
the mean pressure gradient, a typical shear stress
$\tau_{eff}=\beta\lambda\nabla P$ is also defined. The idea is then to use
these quantities in the rheological function $\dot{\gamma}=f(\tau)$ to derive
a generalization of Darcy’s law in the form :

$u\propto f(\beta\nabla P),$ (5)

where the pre-factors must be determined (experimentally, numerically or
theoretically). It is therefore expected that the flow/pressure curve will
keep the overall shape of the rheological curve. We will refer to these
approaches as the “mean field” approach.

Despite its apparent simplicity, it has successfuly been used for different
types of rheologies: see for example Christopher & Middleman (1965); Shah &
Yortsos (1995, 1996); Auriault et al. (2002) for the power laws, Chauveteau
(1982); Sorbie et al. (1989); Lopez et al. (2003); P. et al. (2006); Eberhard
et al. (2019); Zami-Pierre et al. (2016) for the truncated power laws/tile
rheologies or Al-Fariss & Pinder (1987); Chevalier et al. (2014); Shahsavari &
McKinley (2016) for yield stress. We note that most of the studies require a
certain number of fitting parameters. Indeed, it seems quite natural that the
prefactors of the effective quantities depends on the local topology and
rheological parameters.

This approach has, however, some limitations which has been illustrated in the
case of the yield stress fluid (Talon & Bauer, 2013; Chevalier & Talon, 2015;
Liu et al., 2019). Indeed, these models assume that the whole fluid in the
medium yields at the same mean gradient of pressure. The description is
therefore particularly inaccurate for heterogeneous media as the shear rate
(or stress) is heterogeneously distributed. Some regions are in fact easier to
yield than others. This induces a transient flow regime where the material
progressively yields and which cannot be taken into account by the mean field
approach. A similar effect should be expected for any fluid with a change of
rheological behavior: a transitional flow regime should occur as the material
progressively and inhomogeneously changes its local rheology. However, as we
will see later, the effective ”mean field” approach can be considered
approximately valid in the asymptotic limits (high or low flow rate).

The problem of non-Newtonian flow in porous media at the macroscopic scale has
not been the subject of many studies. Entov (1967) studied homogeneous porous
media with different geometries and we studied (Kostenko & Talon, 2019) the
case of the yield stress in heterogeneous porous media. This work has shown
that the influence of the disorder is similar to that at the micro scale:
there is a progressive creation of preferential flowing paths which induce a
transient mean flow regime. The problem has been more studied for fractures
considering different rheologies: Auradou et al. (2008); Boschan et al.
(2008); Lavrov (2013b, a, 2015); Hewitt et al. (2016); Roustaei et al. (2016);
Rodríguez de Castro & Radilla (2017). This studies showed a modification of
the structure of the flow field as a function of the non-linear rheology
(yield stress or shear-thinning) and of the mean aperture. There are also the
works of Di Federico (1998); Talon et al. (2014); Felisa et al. (2018) which
have analytically studied the particular cases where the opening is correlated
in layers either perpendicular or transverse to the direction of the flow. In
these cases, transitional flow regimes are observed depending on the
distribution of the opening. There are also studies on foam flows in porous
media or fracture Rossen (1990); Chen et al. (2005); Géraud et al. (2016)
which can be, in some extend, be considered as non-Newtonian fluids.

It is important to note here that, although we are mainly dealing with non-
Newtonian rheology, much of this work can also be used to understand the role
of inertial effect in porous media or fracture. Indeed, since Forchheimer’s
pioneering work, it is usually accepted that, in order to include the effect
of inertia, a modified Darcy’s law of the following form can be used:
(Forchheimer, 1901; Mei & Auriault, 1991; Andrade Jr et al., 1999):

$\vec{\nabla}P\propto au^{3}+bu^{2}+cu.$ (6)

The latter model can also be considered as a shear thickening rheology with a
smooth transition from Newtonian at low flow rate to shear-thickening power-
law at high flow rate.

The main objective of the article is to investigate transitional rheologies in
macroscopic heterogeneous porous media. In particular, we are interested in
understanding the interplay between the hetergogenities and the non-linearity
of the viscosity. One motivation of this work is for instance to assess the
use of non-Newtonian fluids to characterize the degree of heterogeneity of a
field. Indeed, as each location modifies its behavior at a different applied
pressure difference, recording the evolution of the average velocity can
potentialy give indications on the amplitude of field heterogeneities.

This paper is structured as follow. Section 2 is dedicated to the problem
presentation. In section 3 will be presented some mathematical properties of
the microscopic and macroscopic heterogeneous Darcy’s equation. First, a brief
discussion will be made on the ”mean field” approach at the microscale to
derive the non-Newtonian Darcy’s law. Then will be discussed the problem of
macroscopic heterogenous permeability field. It will be shown for instant
that, if the rheological curve presents power law asymptotes at low and high
flow rate, it is possible to deduce the asymptotes of the relationship between
the mean flow rate and the mean pressure gradient, by considering a power-law
rheology. A perturbation expansion method will be used to predict the
asymptotic mean flow rate and the spatial fluctuations of the velocity field.
Section 4 contains the numerical results using the so-called ”truncated”
rheology, which behaves as a Newtonian fluid at low shear rate and presents a
power law rheology at high shear rate, with a sharp transition. In particular,
the transition regime presents interesting statistical properties, similar to
those of displaying criticality (_e.g_ percolation). Section 5 will be
dedicated to a discussion and conclusion.

##  2 Problem description

We will first recall the basic equations governing the motion of non-
Newtoninan fluids both at the micro and meso (Darcy) scale.

###  2.1 Governing equations

####  2.1.1 Microscale: Stokes equation

At the micro scale, assuming a steady state and negligible inertia, the flow
is governed by the non-Newtonian Stokes equation :

$\nabla\cdot\bm{\Pi}-\nabla P=\bm{0}$ (7)

where $P$ is the pressure and $\bm{\Pi}$ is the deviatoric stress tensor. The
rheology of the fluid is then characterize by the relationship between the
deviatoric stress tensor and the strain rate tensor defined as:

$\bm{\Delta}_{ij}=\frac{1}{2}\left(\frac{\partial{u}_{i}}{\partial{{x}}_{j}}+\frac{\partial{u}_{j}}{\partial{{x}}_{i}}\right),$
(8)

where $\vec{u}$ is the velocity field. For isotropic and incompressible
fluids, it is sufficient to relate the second invariant of these tensors given
by:

$\dot{\gamma}=\sqrt{2\bm{\Delta}_{ij}\bm{\Delta}_{ij}}\;\;{\rm
and}\;\;\tau=\sqrt{\frac{1}{2}\bm{\Pi}_{ij}\bm{\Pi}_{ij}},$ (9)

where $\dot{\gamma}$ and $\tau$ are the shear rate and shear stress
respectively. The rheological response is thus defined by either the function
$\tau(\dot{\gamma})$ or, equivalently, the effective viscosity function:

${{\mu}_{eff}}(\dot{\gamma})=\frac{\tau(\dot{\gamma})}{\dot{\gamma}}.$ (10)

And we have:

$\bm{\Pi}=2\mu(\dot{\gamma})\bm{\Delta}=2\frac{\tau(\dot{\gamma})}{\dot{\gamma}}\bm{\Delta}.$
(11)

The simplest rheology is the Newtonian one, where viscosity is independent of
shear rate. However, there is a wide variety of non-Newtonian rheologies
(Bird, 1976; Coussot, 2005).

One common rheology is the power-law rheology which is defined as
${{\mu}_{eff}}(\dot{\gamma})=K\dot{\gamma}^{n-1}$, where $K$ is the
consistency and $n$ the flow index. If $n<1$, the viscosity decreases with the
shear rate and the fluid is called shear-thinning. On contrary if $n>1$, the
fluid is shear-thickening.

Some fluids might display a change of behavior depending on the shear rate or
stress. For instance, one can cite yield stress fluids often described by the
Herschel-Bulkley rheology:

$\left\\{\begin{array}[]{lll}\tau(\dot{\gamma})<\tau_{0}&\;\;\;{\rm
if}&\dot{\gamma}=0\\\ \tau(\dot{\gamma})=\tau_{0}+K\dot{\gamma}^{n}&\;\;\;{\rm
if}&\dot{\gamma}>0\\\ \end{array}\right.\;\;\;{\rm
or}\;\;\;\left\\{\begin{array}[]{lll}\dot{\gamma}(\tau)=0&\;\;\;{\rm
if}&\tau<\tau_{0}\\\
\dot{\gamma}(\tau)=\left[\frac{1}{K}(\tau-\tau_{0})\right]^{1/n}&\;\;\;{\rm
if}&\tau>\tau_{0}\\\ \end{array}\right.,$ (12)

where $\tau_{0}$ is the yield stress, $n$ the flow index and $K$ the
consistency. If the applied stress is below the yielding threshold, the
material responds like a solid. Above it, the fluid is able to be sheared and
thus flowing.

Many fluids also diplays a change of rheological behavior, they have a
constant viscosity at low shear rate, but as the shear rate increases they
become shear-thinning of thickening with a power law viscosity. There are many
different models to describe those fluids, such as the Carreau-Yasuda (Bird et
al., 1987) model:

$\frac{{{\mu}_{eff}}(\dot{\gamma})-\mu_{\infty}}{\mu-\mu_{\infty}}=[1+(\frac{\dot{\gamma}}{\dot{\gamma}_{0}})^{a}]^{\frac{n-1}{a}},$
(13)

where $\mu$ and $\mu_{\infty}$ are the viscosity at low and high shear rate
respectively, and $\dot{\gamma}_{0}$ a shear rate characterizing the
transition. If $\mu_{\infty}=0$, the fluid is Newtonian at low shear rate and
power-law at higher shear rate. Similarly, there is the Cross model:

${{\mu}_{eff}}(\dot{\gamma})=\frac{\mu}{1+(\frac{\dot{\gamma}}{\dot{\gamma}_{0}})^{n}}.$
(14)

In this article, we will use a simplified model to describe such rheology
which is the “truncated” rheology:

$\left\\{\begin{array}[]{lll}\tau(\dot{\gamma})={{\mu}}\dot{\gamma}&\;\;\;{\rm
if}&\dot{\gamma}<\dot{\gamma}_{0}\\\
\tau(\dot{\gamma})=K\dot{\gamma}^{n}&\;\;\;{\rm
if}&\dot{\gamma}>\dot{\gamma}_{0}\\\ \end{array}\right.\;\;\;{\rm
or}\;\;\;\left\\{\begin{array}[]{lll}\dot{\gamma}(\tau)=\frac{1}{{{\mu}}}\tau&\;\;\;{\rm
if}&\tau<\tau_{0}=\mu\dot{\gamma}_{0}\\\
\dot{\gamma}(\tau)=\left[\frac{1}{K}\tau\right]^{1/n}&\;\;\;{\rm
if}&\tau>\tau_{0}\\\ \end{array}\right.,$ (15)

where $\mu$ is the viscosity at low shear rate and $\dot{\gamma}_{0}$ a
threshold shear rate at which the fluid changes its behavior as sketched in
Fig. 1 (left).

Figure 1:  Left: Relationship between shear rate and shear stress for a
truncated rheology. For a shear rate lower than $\gamma_{0}$ (or a stress
lower than $\tau_{0}$), the fluid is Newtonian. Beyond this thresholds the
fluid is a power law fluid. Right: Extended Darcy’s law, mean flow rate as
function of the mean gradient of pressure in a porous medium, for the
truncated rheology using the “mean field” rheology approach.

####  2.1.2 Mesoscale: Darcy’s law

These rheological laws are applicable to solve Stokes equation and thus to
determine the flow at the microscopic level, inside each pore. As described in
the introduction, in porous media, it is important to determine a mean
constitutive equation, which relates the mean flow rate, $u$ to the mean
pressure gradient, $\nabla P$ at the mesoscale, viz. a volume containing a
large amount of pores. By applying the effective shear rate approach, it
results the following generalization of Darcy’s law:

$\nabla P=-g(u)\;\;\;{\rm or}\;\;\;u=-f(\nabla P),$ (16)

where $g$ and $f$ are assumed to be monotonically increasing (and thus
invertible) functions.

To solve the flow field, a vector formulation of the non-Newtonian Darcy’s law
is required. Assuming that the medium is locally isotropic, so that the mean
flow is collinear and opposite to the mean pressure gradient, it follows:

$\displaystyle\vec{\nabla}P=-g(\vec{r};|\\!|\vec{u}|\\!|)\frac{\vec{u}}{|\\!|u|\\!|}$
$\displaystyle\;\;\;{\rm or}\;\;\;$
$\displaystyle\vec{u}=-f(\vec{r};|\\!|\vec{\nabla}P|\\!|)\frac{\vec{\nabla}P}{|\\!|\vec{\nabla}P|\\!|},$
(17) $\displaystyle\rm{and}$ (18) $\displaystyle\vec{\nabla}.\vec{u}=0$

This form is very generic and applies for instance to Carreau, Herschel-
Bulkley or the truncated rheology. It is also very important to note that the
functions $g$ and $f$ may vary spatially (_e.g._ heterogeneous permeability
field). For better readability, the $\vec{r}$ will be omitted in the
following.

For the truncated rheology, the equation then reads:

$\left\\{\begin{array}[]{lll}\vec{\nabla}P=-\frac{{{\mu}}}{\kappa(\vec{r})}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|<{u_{c}}(\vec{r})\\\
\vec{\nabla}P=-\frac{{{\mu}}}{\kappa(\vec{r})}\left[\frac{|\\!|\vec{u}|\\!|}{{u_{c}}(\vec{r})}\right]^{n-1}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|>{u_{c}}(\vec{r})\\\ \end{array}\right.,$ (19)

where $\kappa(\vec{r})$ is the local permeability and the prefactor in the
second relation has been determined by continuity. Equivalently, by defining
${u_{c}}=\frac{\kappa{g_{c}}}{{{\mu}}}$, one may also write:

$\left\\{\begin{array}[]{lll}\vec{u}=-\frac{\kappa(\vec{r})}{{{\mu}}}\vec{\nabla}P&\;\;\;{\rm
if}&|\\!|\nabla P|\\!|<{g_{c}}(\vec{r})\\\
\vec{u}=-\frac{\kappa(\vec{r})}{{{\mu}}}\left[\frac{|\\!|\vec{\nabla}P|\\!|}{{g_{c}}(\vec{r})}\right]^{1/n-1}\vec{\nabla}P&\;\;\;{\rm
if}&|\\!|\nabla P|\\!|>{g_{c}}(\vec{r})\\\ \end{array}\right.,$ (20)

As described in the introduction, this approach is based on the assumption
that an effective shear rate $\gamma_{eff}$ can be related to an effective
shear stress $\tau_{eff}$. Despite its relative simplicity, this approach has
been validated by my numerous experimental and numerical works (see
introduction for references). In section 3.1, we will examine in more detail
the validity of this approach. For example, we will show that these
expressions actually represent the asymptotic limit at low and high flow rate.

In this article, where we intend to study the flow in a heterogeneous
macroscopic field, eq. (19) is assume to be valid but with a hetergenous
permeability field $\kappa(\vec{r})$ and threshold velocity field
${u_{c}}(\vec{r})$ (or ${g_{c}}(\vec{r})$) to be heterogenous.

###  2.2 Porous media

The permeability field was chosen to be distributed according to a log-normal
distribution which is a common assumption (Gelhar & Axness, 1983). This is
obtained by generating a Gaussian field $\delta f$ with zero mean and a given
standard deviation $\sigma_{f}$. $\delta f$ is correlated with a correlation
length ${\lambda}$ (see (Kostenko & Talon, 2019) for more details). The
permeability field is then given by:

$\kappa=\exp{(f_{0}+\delta f)}=\kappa_{0}\exp{(\delta f)},$ (21)

where $\kappa_{0}=\exp{(f_{0})}$, is a parameter characterizing the average
permeability of the medium.

As for the threshold parameters ${u_{c}}$ and ${g_{c}}$, intuitively they are
expected to be somehow related to the permeability $\kappa$. A relationship
can be established using phenomenological arguments. For porous media,
defining the typical pore size $d$, a scaling analysis leads to $\kappa\sim
d^{2}$ and ${g_{c}}\sim\tau_{0}/d$, thus: ${g_{c}}\sim\kappa^{-1/2}$ and
${u_{c}}\sim\kappa^{1/2}$.

For fracture, this reasoning is slighlty different, defining $h$ is a typical
opening leads to $\kappa\sim h^{3}$ and ${g_{c}}\sim\tau_{0}/h$. This leads
to: ${g_{c}}\sim\kappa^{-1/3}$ and ${u_{c}}\sim\kappa^{2/3}$. As these
arguments are purely phenomenological, we have studied a more general relation
in the form

${u_{c}}=A\kappa^{\gamma},$ (22)

where $A$ is a prefactor and $\gamma$ is a parameter with $\gamma\in[0;1]$. We
also have:

${g_{c}}=\mu A\kappa^{\gamma-1}.$ (23)

###  2.3 Simulations

Macroscopic simulations were performed using a second order finite difference
method combined with an augmented Lagragian approach described in the
appendix. In practice, the flow rate is driven by imposing a pressure
difference $\Delta P$ between inlet and outlet, but it is more convenient to
use the mean pressure gradient $\langle\nabla P\rangle=\frac{\Delta P}{L}$ as
control parameter, where $L$ is the length of the system. Lateral boundary
conditions are periodic. Althugh there are many different parameters, the
focus of this study will be on the amplitude of the disorder $\sigma_{f}$, the
flow index $n$ and the exponent $\gamma$. Parameters $A$, $\mu$ and
$\kappa_{0}$ will thus be kept constant with $A=1$, $\mu=0.1$ and
$\kappa_{0}=1$. Based on these parameters a characteristic velocity

${u_{0}}=A\kappa_{0}^{\gamma}$

and a characteristic pressure gradient

${g_{0}}=\frac{\mu{u_{0}}}{\kappa_{0}}$

can be defined and be used for non-dimensionalization. ${u_{0}}$ and ${g_{0}}$
represent the velocity and the pressure gradient at which the system would
change its behavior if the field were homogeneous (_i.e_ $\sigma_{f}=0$).

###  2.4 Flow examples

####  2.4.1 Homogenous permeability - heterogeneous critical velocities

In order to describe qualitatively the problem, a simplified version is
presented, where the permeability is homogeneous but not the critical velocity
${u_{c}}$. The latter varies spatially according to a log-normal distribution.

In figure 2 (top) shows the flow field for different mean pressure gradients
$\langle\nabla P\rangle$. For a small pressure difference, the flow field is
homogeneous as for a Newtonian Darcy flow in a homogeneous field. As the
pressure difference increases, flow field remains homogeneous until the first
location reaches its local threshold (Fig. 2.a). Around this first point,
since the viscosity has changed, the flow field is then perturbed in a way
similar to a dipole flow. In this figure, the second regime is shear-thinning
$n<1$ so the viscosity has locally decreased. The flow has thus increased
downstream and upstream, while it has decreased on the lateral sides.
Consequently, the change of regime at one point should favour or defavour the
change of behaviour in the vicinity. One can therefore expect some
correlations of the regions that have changed their flow regime. For shear-
thinning (resp. thickening) fluids the change of regime should thus be
correlated in the streamwise (resp. perpendicular) direction. Figure 2
(bottom) represents the correspondings regions which have change their regime.
It can be noted that this correlation effect seems in fact quite marginal.
This might be explained by the fact that the distribution of ${u_{c}}$ is
quite large $\sigma=1$. The dipole velocity perturbation might thus be unable
to compensate large fluctuation of the threshold value. As the pressure
difference increases, more and more regions change their behavior until the
entire domain is in the non-Newtonian regime.

Figure 3 presents the evolution of different interesting quantities. Figure
3.a diplays the mean velocity as a function of the mean pressure gradient. As
expected, for low pressure, the average flow rate varies linearly with the
pressure as expected. At a certain point, the curve begins to deviate slowly
from the linear trend to reach a power-law trend. This change of behavior can
be magnified by plotting the amount ${\langle u\rangle}/\langle\nabla
P\rangle$. We observe three flow regimes, the first regime is one where the
flow is linear. At a high flow rate, the flow follows a power law regime of
the exponent ${\langle u\rangle}\propto\langle\nabla P\rangle^{1/n}$. Between
these two regimes, there is a transitional regime which is the main subject of
this article. A convenient quantity to characterize this transition is the
percentage of regions being in the non-linear regime
$\mathcal{O}(\langle\nabla P\rangle)$ ranging from $0$ to $1$ (Fig. 3.c).
Another important quantity, shown in figure 3.d, is the relative standard
deviation of the velocity field, $\frac{\sigma_{u}}{{\langle u\rangle}}$,
which characterizes the flow heterogeneities and is of great importance in the
problem of transport of species. Here, we observe that the standard deviation
is zero at low pressure because the permeability is homogeneous. After the
transition, it reaches a constant non-zero value. The reason behind this
observation is that, even if the permeability is homogeneous, each location
changes its behavior at a different mean pressure gradient. It thus introduces
heterogeneities in the flow field.

Figure 2:  For a shear-thinning truncated rheology $n=1/2$ in a medium with a
homogeneous permeability field but with a heterogeneous threshold, $\sigma=1$.
Top: Velocity field colormap for different imposed mean pressure gradients
increasing from left to right. The average flow direction is from top to
bottom. Middle row: in black are represented the regions above the threshold.
Bottom: Distribution function of the velocity (blue) and threshold velocity
${u_{c}}$ (orange).

Figure 3:  For the medium of Fig. 3 as function of the imposed mean gradient
of pressure: a: Mean velocity $\langle u\rangle$. b: Mean velocity rescaled
with the Darcy’s law. c: Percentage of regions which have change their
behavior (_i.e_ $u>{u_{c}}$). d: relative standart deviation of the velocity
field $\frac{\sigma_{u}}{{\langle u\rangle}}$. Vertical dashed lines represent
the transitions of the flow regimes.

Figure 4:  For a shear-thinning truncated rheology $n=1/2$ in a medium with a
heterogeneous permeability field, with $\sigma_{f}=1$, but with a constant
velocity threshold $u_{c}$ (_i.e_ $\gamma=0$). Top: Velocity field colormap
for different imposed mean pressure gradients increasing from left to right.
The average flow direction is from top to bottom. Middle row: in black are
represented the regions above the threshold. Bottom: Evolution of the velocity
(blue) and threshold velocity (orange) distribution function.

####  2.4.2 Heterogeneous permeability - homogeneous critical velocities

Another interesting example is the opposite case, where the permeability
$\kappa(\vec{r})$ is heterogeneous but not the critical velocity ${u_{c}}$,
corresponding to $\gamma=0$. Figure 4 represents the flow field distribution
for different applied pressure difference. At low pressure drop, the velocity
field corresponds to the heterogeneous Darcy’s equation and is thus
heterogenous as well. As the mean flow rate increase, the first region to
change its viscosity will be the one with the highest velocity in the
Newtonian case. The velocity distribution is then drastically change above
${u_{c}}$ but the topology of the regions $u>{u_{c}}$ reflects the properties
of the velocity field of the heterogeneous Newtonian Darcy’s equation. In
particular, the correlations are in the streamwise direction as it will be
discussed later.

##  3 Analytical aspects

In this setion, we discuss some analytical aspects of the problem. It is
important to note that many of the results presented here are very general,
and can be applied to other non-linear rheologies than the truncated one.

###  3.1 Energy balance - Virtual work

The energy balance and the virtual work approaches lead to interesting
properties. Virtual work consists in evaluating the work of the real forces on
a virtual flow field, which is not necessarily the real one. Among all the
possible velocity field, the solution of the problem is the minimum of a
certain functional to be defined. This approach will be applied at two scales:
the microscopic scale, where the fluid obeys Stokes non-Newtonian flow, and
the macroscopic scale where the fluid obeys the non-Newtonian Darcy’s law.

At the microscopic level, this approach was for example introduced by Duvaut &
Lions (1976) for yield stress fluids. Here, the description is summarized as
it has already been discussed for instance for Hershel Bulkley in porous media
in the appendix of Bauer et al. (2019). The method is however extended to any
rheology. At the microscopic level, two limites, at high and low flow, can be
formally predicted. This can partialy justify the mean field approach in order
to generalize Darcy’s law.

At the macroscopic level, to our knowledge, it has never been applied, so the
method will presented in more details. Here again, it allows to predict
asymptotic behavior of the flow.

###  3.2 Microscopic description

A porous parallelepiped medium containing solid and porous regions is
considered. The entry section is $S$, the length $L$ and a characteristic size
(average pore, grain diameter, etc.) is named $\lambda$. A pressure $P_{in}$
and $P_{out}$ is applied to the inlet and outlet respectively. The total flow
rate through $S$ is denoted $Q$. No flow on the lateral side of the porous
media is assumed. In the porous regions, the rheology of the fluid is defined
by the function $\tau(\dot{\gamma})$, a monotonically increasing function.
Following Chevalier et al. (2014), a standard energy balance gives :

$Q\Delta P=\int\tau(\dot{\gamma})\dot{\gamma}\;dr^{3}.$ (24)

This equation indicates that the relationship between the mean flow rate and
the mean applied pressure can be deduced from the shear rate distribution
$\dot{\gamma}(\vec{r})$. Its determination requires to solve the Stokes
equation. Since pressure can be considered as an energy per unit volume,
$Q\Delta P$ represents the difference between the incoming and outgoing flux
of this energy. It should therefore be equal to the total viscous dissipation
within the medium given by the righthand term of eq. (24).

If the solution $\dot{\gamma}$ is determined by Stokes equation, it can also
be obtained from a minimum principle: among all admissible velocity fields
$\vec{v}$ (_i.e._ satisfying the free divergence and the boundary condition),
the solution $\vec{u}$ of the Stokes equation minimizes the functional:

$\vec{u}=\arg\min_{\vec{v}}\left\\{\int\phi(\dot{\gamma}[\vec{u}])\;dr^{3}-\Delta
P\;Q[\vec{v}]\right\\},$ (25)

where
$\phi(\dot{\gamma})=\int_{0}^{\dot{\gamma}}\tau(\dot{\gamma}^{\prime})d\gamma^{\prime}$.
$Q[\vec{v}]=\int_{inlet}\vec{v}.\vec{dS}$ the total flow rate of the $\vec{v}$
field. Physically, this relation can be considered as an extension of the
minimum dissipation principle for the Newtonian Stokes equation (Keller et
al., 1967; Happel & Brenner, 1983). It is important to note that this minimum
principle is applicable for an imposed pressure difference $\Delta P$. A
similar result can also be obtained for an imposed flow rate $Q$ (Bauer et
al., 2019).

####  3.2.1 Power-law fluids

Considering a simple power-law fluid, $\tau(\dot{\gamma})=K\dot{\gamma}^{n}$,
it follows that $\phi(\dot{\gamma})=K\frac{1}{{n}+1}\dot{\gamma}^{{n}+1}$.
Using eq. (25), the velocity field is the solution of:

$\vec{u}=\arg\min_{\vec{v}}\left\\{\int
K\frac{1}{{n}+1}\dot{\gamma}^{{n}+1}[\vec{v}]\;dr^{3}-\Delta
PQ[\vec{v}]\right\\}.$

It is then important to note that the result of the minimum should not change
by multiplying the function by any real number. Multiplying by
${\epsilon}^{n+1}$ leads to:

$\vec{u}=\arg\min_{\vec{v}}\left\\{\int
K\frac{1}{{n}+1}\dot{\gamma}^{{n}+1}[{\epsilon}\vec{v}]\;dr^{3}-({\epsilon}^{n}\Delta
P)Q[{\epsilon}\vec{v}]\right\\},$

due to the linearity of the operator $\dot{\gamma}[\vec{v}]$ and $Q[\vec{v}]$.
And thus:

$\vec{u}=\frac{1}{{\epsilon}}\arg\min_{\vec{v}}\left\\{\int
K\frac{1}{{n}+1}\dot{\gamma}^{{n}+1}[\vec{v}]\;dr^{3}-({\epsilon}^{n}\Delta
P)Q[\vec{v}]\right\\},$

by changing the variable. This indicates that the velocity field, solution of
the problem for the imposed pressure ${\epsilon}^{n}\Delta P$ is the same
solution at $\Delta P$, but multipied by ${\epsilon}$. It follows the scaling
relationships:

$\displaystyle\Delta P$ $\displaystyle\rightarrow$
$\displaystyle{\epsilon}^{n}\Delta P$ $\displaystyle\vec{u}$
$\displaystyle\rightarrow$ $\displaystyle{\epsilon}\vec{u}$ $\displaystyle Q$
$\displaystyle\rightarrow$ $\displaystyle{\epsilon}Q.$

It must be noted that this scaling relationship agrees with the first order
expansion obtain by homogenization (Shah & Yortsos, 1995; Auriault et al.,
2002). For any power-law rheology of index ${n}$, the field
${}_{{n}}{\vec{\omega}}=\frac{\vec{u}}{Q[\vec{u}]}$ is thus a constant
independant of the applied pessure. And the non-dimensional shear rate field,
defined as ${{}_{{n}}\delta}=\frac{\lambda S\dot{\gamma}}{Q}$ is also
invariant with the applied pressure.

From the energy balance equation, eq. (24), it follows:

$\displaystyle\Delta P$ $\displaystyle=$
$\displaystyle\frac{KQ^{n}}{\lambda^{{n}+1}S^{{n}+1}}\int{{}_{{n}}\delta}^{{n}+1}dr^{3}$
(26) $\displaystyle=$
$\displaystyle\frac{KQ^{n}L}{\lambda^{{n}+1}S^{n}}\langle({{}_{{n}}\delta})^{{n}+1}\rangle.$
(27)

And so

$\lambda\frac{\Delta P}{L}=K\left(\frac{Q}{\lambda
S}\right)^{n}\langle({{}_{{n}}\delta})^{{n}+1}\rangle.$ (28)

This equation shows that the mean field approach is thus able to predict
correctly flow rate - pressure curve, providing the introduction of
$\tau_{eff}=\lambda\frac{\Delta P}{L}$ and
$\dot{\gamma}_{eff}=\frac{Q}{\lambda
S}\sqrt[{n}]{\langle{{}_{{n}}\delta}^{{n}+1}\rangle}=\frac{\langle
u\rangle}{\lambda}\sqrt[{n}]{\langle{{}_{{n}}\delta}^{{n}+1}\rangle}$. The
crucial point for this approach to be valid, is that the quantity
$\sqrt[{n}]{\langle{{}_{{n}}\delta}^{{n}+1}\rangle}$ is a constant field
independant of the applied pressure (or flow rate). This field depends on the
porous structure, so it requires to solve the Stokes equation once.

####  3.2.2 Rheology with change of behavior

For more complex rheologies, such as the truncated one, there is no guarantee
that $\sqrt[{n}]{\langle{}_{{n}}g^{{n}+1}\rangle}$ is a constant field
independant of the boundary conditions. However, it is interesting to consider
the case where the rheological function $\tau(\dot{\gamma})$ admits power-law
asymptotes:

$\left\\{\begin{array}[]{lll}\tau(\dot{\gamma})\sim
K_{0}\dot{\gamma}^{a}&\;\;{\rm if}&\dot{\gamma}\rightarrow 0\\\
\tau(\dot{\gamma})\sim K_{\infty}\dot{\gamma}^{b}&\;\;{\rm
if}&\dot{\gamma}\rightarrow\infty.\end{array}\right.$ (29)

Since $\dot{\gamma}\sim Q$, it is expected that at high (resp. low) flow rates
the highest (resp. lowest) exponent will dominate in $\phi(\dot{\gamma})$, so
that

$\left\\{\begin{array}[]{lll}\phi(\dot{\gamma})\xrightarrow{Q\to
0}\dot{\gamma}^{a+1}\\\
\phi(\dot{\gamma})\xrightarrow{Q\to\infty}\dot{\gamma}^{b+1}.\\\
\end{array}\right.$ (30)

Using these limits in the minimization equation eq. (25), indicates that the
solutions are equivalent to solving the flow for a power-law fluid with the
dominatant exponent. It follows that the fields $u/Q$ and the non-dimensional
shear rate $\frac{\lambda S\dot{\gamma}}{Q}$ become independent of the applied
pressure at sufficiently high and low applied pressure (or flow rate). The
knowledge of these two quantities allows then to derivate the effective
Darcy’s law in the two limits, using the energy balance eq. (24)

As an example, we assume a polynomial rheology:

$\tau(\dot{\gamma})=A\dot{\gamma}^{3}+B\dot{\gamma}^{2}+C\dot{\gamma}+D.$ (31)

It follows:

$\phi(\dot{\gamma})=\frac{1}{4}A\dot{\gamma}^{4}+\frac{1}{3}B\dot{\gamma}^{3}+\frac{1}{2}\dot{\gamma}^{2}+D\dot{\gamma}.$

At high (resp. low) flow rate, the non-dimensional shear rate will thus
converge to the one obtained with a power-law fluid of
${{}_{3}\delta}(\vec{r})$ (resp. ${{}_{1}\delta}(\vec{r})$). Eq. (24) reads

$Q\Delta
P=A\int\dot{\gamma}^{4}dr^{3}+B\int\dot{\gamma}^{3}dr^{3}+C\int\dot{\gamma}^{2}dr^{3}+D\int\dot{\gamma}dr^{3}.$

Using the non-dimensional shear rate field ${{}_{{n}}\delta}=\frac{\lambda
S\dot{\gamma}}{Q}$ , it follows the two asymptotic behaviors:

$\displaystyle\lambda\frac{\Delta P}{L}=A\left(\frac{{\langle
u\rangle}}{\lambda}\right)^{3}\langle{{}_{3}\delta}^{4}\rangle+B\left(\frac{{\langle
u\rangle}}{\lambda}\right)^{2}\langle{{}_{3}\delta}^{3}\rangle+C\left(\frac{{\langle
u\rangle}}{\lambda}\right)\langle{{}_{3}\delta}^{2}\rangle+D\langle{{}_{3}\delta}\rangle$
$\displaystyle\;\;\;\rm{for}\;\;\;$ $\displaystyle{\langle
u\rangle}\rightarrow\infty$ $\displaystyle\lambda\frac{\Delta
P}{L}=A\left(\frac{{\langle
u\rangle}}{\lambda}\right)^{3}\langle{{}_{1}\delta}^{4}\rangle+B\left(\frac{{\langle
u\rangle}}{\lambda}\right)^{2}\langle{{}_{1}\delta}^{3}\rangle+C\left(\frac{{\langle
u\rangle}}{\lambda}\right)\langle{{}_{1}\delta}^{2}\rangle+D\langle{{}_{1}\delta}\rangle$
$\displaystyle\;\;\;\rm{for}\;\;\;$ $\displaystyle{\langle
u\rangle}\rightarrow 0,$

with the notation $\langle.\rangle=\frac{1}{SL}\int.dr^{3}$ and $Q=S{\langle
u\rangle}$. Here, both limits have a polynomial form, similar to the
rheological function $\tau(\dot{\gamma})$. The coefficients depend on the
k-$th$ moments $\langle{{}_{{n}}\delta}^{k}\rangle$ which are constant but
require solving a power law rheology with index $1$ and $3$.

This equation leads to a form similar to the mean field approach:

$\beta\lambda\frac{\Delta P}{L}=A(\frac{\alpha{\langle
u\rangle}}{\lambda})^{3}+B(\frac{\alpha{\langle
u\rangle}}{\lambda})^{2}+C\frac{\alpha{\langle u\rangle}}{\lambda}+D.$ (32)

However, there are two notable differences.

First of all, matching the main order term requires to define
$\dot{\gamma}_{eff}=\frac{{\langle
u\rangle}}{\lambda}\langle{{}_{3}\delta}^{4}\rangle^{1/3}$ at high flow rate
and $\dot{\gamma}_{eff}=\frac{{\langle
u\rangle}}{\lambda}(\langle{{}_{1}\delta}^{4}\rangle)^{1/3}$, at low flow
rate. The definition of $\dot{\gamma}_{eff}$ should thus change in the two
limits. Secondly, the other coefficients then do not necessarily match. For
example, there is a mismatch in the second term if
$\langle{{}_{3}\delta}^{4}\rangle^{1/3}\neq\langle{{}_{3}\delta}^{3}\rangle^{1/2}$.
It is therefore generally not possible to determine $\alpha$ and $\beta$ so
that all coefficients match, unless the polynomial has only two terms (as for
Herschel-Bulkley).

This approach is applicable when the rheological function $\tau(\dot{\gamma})$
has a power-law limit $\tau(\dot{\gamma})\sim K_{0}\dot{\gamma}^{a}$ (resp.
$\sim K_{0}\dot{\gamma}^{b}$) in the high (resp. low) flow rate limit. It is
thus applicable with the truncated rheology ($a=1$ and $b=n$) but also for
yield stress fluids (for Herschel-Bulkley, $a=0$ and $b=n$). It has been used
by Chevalier et al. (2014) to propose a mesoscopic Darcy’s law for Herschel-
Bulkley. The authors also used the Newtonian velocity distribution even for
$n\neq 1$ to estimate the coefficients. The agreement with the experiments
seems then to suggest that
$\langle{{}_{n}\delta}^{n+1}\rangle^{1/n}\sim\langle{{}_{1}\delta}^{2}\rangle$,
at least in the porous media used (bead packing).

###  3.3 Macroscopic description

In this section, a general nonlinear Darcy equation given by eqs. (16) is
assumed. It will be demonstrated that the solution of the heterogenous non-
linear Darcy’s equation is also the minimum of a functional $\Phi[\vec{v}]$,
among all the admissible field $\vec{v}$. For this geometry, the pressure is
prescribed at the inlet edge $P_{in}$ and outlet edge $P_{out}$. No flux flow
through the lateral boundaries is assumed.

First, we define $\Omega$ the set of admissible velocity fields satisfying the
divergent free and the boundary conditions, Multiplying eq. (17) by any
$\vec{v}\in\Omega$, and integrating over the domain, yields to:

$\forall\vec{v}\in\Omega\;,\;\int\vec{\nabla}P.\vec{v}\;dr^{3}=-\int
g(|\\!|\vec{u}|\\!|)\frac{\vec{u}.\vec{v}}{|\\!|\vec{u}|\\!|}dr^{3}$ (33)

Using the free divergence of $\vec{v}$ and the divergent theorem, it follows

$\displaystyle\forall\vec{v}\in\Omega\;,\;-\int
g(|\\!|\vec{u}|\\!|)\frac{\vec{u}.\vec{v}}{|\\!|\vec{u}|\\!|}dr^{3}$
$\displaystyle=$ $\displaystyle(P_{out}-P_{in})\int_{outlet}\vec{v}.\vec{dS}$
(34) $\displaystyle=$ $\displaystyle-Q[\vec{v}]\Delta P,$

where $Q[\vec{v}]=\int_{outlet}\vec{v}.\vec{dS}=-\int_{inlet}\vec{v}.\vec{dS}$
is the total flow rate of the field $\vec{v}$ and $\Delta P=P_{in}-P_{out}$.

It can be demonstrated that the field $\vec{u}$, solution of the Darcy’s
equation eqs. (16), is also the minimum among, all admissible
$\vec{v}\in\Omega$, of $\Phi[\vec{v};\Delta P]$ defined as:

$\Phi[\vec{v};\Delta P]=\int G(|\\!|\vec{v}|\\!|)\;dr^{3}-\Delta PQ[\vec{v}].$
(35)

Where $G(|\\!|\vec{v}|\\!|)$ is the primitive of $g(|\\!|\vec{v}|\\!|)$:

$G(|\\!|\vec{v}|\\!|)=\int_{0}^{|\\!|\vec{v}|\\!|}g(y)dy.$

It is equivalent to prove that:

$\forall\vec{v}\in\Omega\;,\;\Phi[\vec{v}+\vec{u};\Delta
P]-\Phi[\vec{u};\Delta P]\geq 0.$ (36)

Indeed, eq. (34) leads to:

$\forall\vec{v}\in\Omega\;,\;\Phi[\vec{v}+\vec{u};\Delta
P]-\Phi[\vec{u};\Delta
P]=\int\left\\{G(|\\!|\vec{v}+\vec{u}|\\!|)-G(|\\!|\vec{u}|\\!|)-\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}\vec{u}.\vec{v}\right\\}\;dr^{3}.$
(37)

Since $g$ is an increasing function, $G$ is convex, and thus:

$G(|\\!|\vec{v}+\vec{u}|\\!|)-G(|\\!|\vec{u}|\\!|)\geq
g(|\\!|\vec{u}|\\!|)(|\\!|\vec{u}+\vec{v}|\\!|-|\\!|\vec{u}|\\!|).$ (38)

It follows the required property:

$\displaystyle\forall\vec{v}\in\Omega\;,\;\Phi[\vec{v}+\vec{u};\Delta
P]-\Phi[\vec{u};\Delta P]$ $\displaystyle\geq$
$\displaystyle\int\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}\left[|\\!|\vec{u}+\vec{v}|\\!||\\!|\vec{u}|\\!|-|\\!|\vec{u}|\\!|^{2}-\vec{u}.\vec{v}\right]\;dr^{3}$
(39) $\displaystyle\geq$ $\displaystyle 0,$

where Cauchy-Schwartz inequality has been used
$(\vec{u}+\vec{v}).\vec{u}\leq|\\!|\vec{u}+\vec{v}|\\!||\\!|u|\\!|$.

Several remarks can be made:
Remark 1 As for the microscopic case, it is instructive to put the solution
$\vec{u}$ in eq. (34). Leading to the energy balance:

$Q[\vec{u}]\Delta P=\int g(|\\!|\vec{u}|\\!|)|\\!|\vec{u}|\\!|dr^{3}>0.$ (40)

This expression can be seen as a balance of energy and it also proves the
intutive result that the flow rate and the difference of pressure
$P_{in}-P_{out}$ have always the same sign.

Remark 2: Several examples can be given for different rheologies:

  * •

Newtonian (Darcy):

$\nabla P=-\frac{{{\mu}}}{\kappa(\vec{r})}\vec{u}.$ (41)

In this case,
$g(|\\!|\vec{u}|\\!|)=\frac{{{\mu}}}{\kappa(\vec{r})}|\\!|\vec{u}|\\!|$ yields
to:

$\Phi[\vec{v};\Delta
P]=\int\frac{{{\mu}}}{2\kappa(\vec{r})}|\\!|\vec{v}|\\!|^{2}-\Delta
PQ[\vec{v}].$

This expression is interesting because it shows the classical result that, in
order to minimize dissipation, it is more favorable to have a higher velocity
where permeability is high. However, it is important to note that the
allowable field $\vec{v}$ must satisfy the conservation of mass. This
constraint can cause low permeability regions to have a high velocity (and
vice versa). Conservation of mass tends to correlate the velocity field in the
streamwise direction.

  * •

Herschel-Bulkley: At Darcy’s scale, the velocity field can be described by
(see Kostenko & Talon (2019)):

$-\nabla
P=c(\vec{r})|\\!|\vec{u}|\\!|^{n-1}\vec{u}+{g_{c}}(\vec{r})\frac{\vec{u}}{|\\!|\vec{u}|\\!|}$
(42)

where ${g_{c}}(\vec{r})$ is the local critical pressure gradient below which
there is no flow, $c(\vec{r})$ is a prefactor that depend on the consistency
and the local geometry, and $n$ the flow index.

It then follows
$g(|\\!|\vec{u}|\\!|)=c(\vec{r})|\\!|\vec{u}|\\!|^{n}+{g_{c}}(\vec{r}).$ Thus:

$\Phi[\vec{v};\Delta
P]=\int\left[\frac{c(\vec{r})}{n+1}|\\!|\vec{v}|\\!|^{n+1}+{g_{c}}(\vec{r})|\\!|\vec{v}|\\!|\right]dr^{3}-\Delta
PQ[\vec{v}].$

  * •

Truncated rheology:

$\left\\{\begin{array}[]{lll}\vec{\nabla}P=-\frac{{{\mu}}}{\kappa}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|<{u_{c}}(\vec{r})\\\
\vec{\nabla}P=-\frac{{{\mu}}}{\kappa}\left[\frac{|\\!|\vec{u}|\\!|}{{u_{c}}}\right]^{n-1}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|>{u_{c}}(\vec{r})\\\ \end{array}\right.,$ (43)

leads to:

$\left\\{\begin{array}[]{lll}\Phi[\vec{v};\Delta
P]=\int\left[\frac{{{\mu}}}{2\kappa(\vec{r})}|\\!|\vec{v}|\\!|^{2}\right]dr^{3}-\Delta
PQ[\vec{v}]&\;\;\;{\rm if}&|\\!|\vec{v}|\\!|<{u_{c}}(\vec{r})\\\
\Phi[\vec{v};\Delta
P]=\int\left[\frac{{{\mu}}}{(n+1)\kappa|\\!|\vec{}{u_{c}}|\\!|^{n-1}}\left[|\\!|\vec{v}|\\!|^{n+1}-{u_{c}}^{n+1}\right]+\frac{{{\mu}}}{2\kappa(\vec{r})}{u_{c}}^{2}\right]-\Delta
PQ[\vec{v}]&\;\;\;{\rm if}&|\\!|\vec{v}|\\!|>{u_{c}}(\vec{r})\\\
\end{array}\right.,$ (44)

  * •

Power-law rheology: Darcy’s law now writes:

$-\vec{\nabla}P=c(\vec{r})|\\!|\vec{u}|\\!|^{{n}-1}\vec{u}.$ (45)

leading to:

$\vec{u}=\arg\min_{\vec{v}}\left\\{\int\frac{c(\vec{r})}{{n}+1}|\\!|\vec{v}|\\!|^{{n}+1}dr^{3}-\Delta
PQ[\vec{v}]\right\\}.$ (46)

Similarly to the microscopic scale, the invariance of the velocity field can
be proved. For any ${\epsilon}$, multiplying by ${\epsilon}^{{n}+1}$ the
previous equation gives:

$\displaystyle\vec{u}(\Delta P)$ $\displaystyle=$
$\displaystyle\arg\min_{\vec{v}}\left\\{\int\frac{c(\vec{r})}{{n}+1}|\\!|{\epsilon}\vec{v}|\\!|^{{n}+1}dr^{3}-{\epsilon}^{n}\Delta
PQ[{\epsilon}\vec{v}]\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{{\epsilon}}\arg\min_{\vec{v}}\left\\{\int\frac{c(\vec{r})}{{n}+1}|\\!|\vec{v}|\\!|^{{n}+1}dr^{3}-{\epsilon}^{n}\Delta
PQ[\vec{v}]\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{{\epsilon}}\vec{u}({\epsilon}^{n}\Delta P).$

It follows the scaling relationships:

$\displaystyle\Delta P$ $\displaystyle\rightarrow$
$\displaystyle{\epsilon}^{n}\Delta P$ $\displaystyle\vec{u}$
$\displaystyle\rightarrow$ $\displaystyle{\epsilon}\vec{u}$ $\displaystyle Q$
$\displaystyle\rightarrow$ $\displaystyle{\epsilon}Q.$

Similarly, the field ${{}_{{n}}\vec{\omega}}=\frac{\vec{u}}{\langle u\rangle}$
is thus a constant field, independant of the applied pressure or the mean
velocity.

As for the microscopic case, it follows that for any rheology with power-law
asymptotes, in the limit of a low or high flow rate, the
$\frac{\vec{u}(\vec{r})}{{\langle u\rangle}}$ distribution should tend towards
a field independent of the mean pressure or flow rate gradient. These constant
fields are determined by solving a power-law fluid with the corresponding
dominant exponent. This is true for Herschel-Bulkley or truncated rheology but
applies to any polynomial function (or series of power-laws).

One example could be the Forchheimer’s law, which corresponds to the
generalization of Darcy’s law including the influence of inertia. It usualy
proposed a relation of the form (Forchheimer, 1901):

$\nabla P=Au^{3}+Bu^{2}+Cu,$ (47)

and thus
$g(|\\!|\vec{u}|\\!|)=A|\\!|\vec{u}|\\!|^{3}+B|\\!|\vec{u}|\\!|^{2}+C|\\!|\vec{u}|\\!|$.
In the high flow rate regime, the velocity field
$\frac{\vec{u}(\vec{r})}{{\langle u\rangle}}$ should tend to the solution of a
shear-thickening fluid eq. (45) with the exponent ${n}=3$.

Using the energy balance equation eq. (40), it follows the two assymptotic
regimes:

$\displaystyle\frac{\Delta
P}{L}=A\langle{{}_{1}\vec{\omega}}^{4}\rangle\langle
u\rangle^{3}+BA\langle{{}_{1}\vec{\omega}}^{3}\rangle\langle
u\rangle^{2}+C\langle{{}_{1}\vec{\omega}}^{2}\rangle\langle
u\rangle^{1}\;\;\text{for}\;\;$ $\displaystyle\langle u\rangle\rightarrow 0$
$\displaystyle\frac{\Delta
P}{L}=A\langle{{}_{3}\vec{\omega}}^{4}\rangle\langle
u\rangle^{3}+BA\langle{{}_{3}\vec{\omega}}^{3}\rangle\langle
u\rangle^{2}+C\langle{{}_{3}\vec{\omega}}^{2}\rangle\langle
u\rangle^{1}\;\;\text{for}\;\;$ $\displaystyle\langle
u\rangle\rightarrow\infty,$

where the $j$th-moment $\langle{{}_{n}\vec{\omega}}^{j}\rangle$ requires to
determine the flow of a power-law fluid in the medium.

In conclusion, any nonlinear Darcy’s law with power law asymptotes converges
to the solution of a power law fluid in heterogeneous media at low or high
flow rate. The next section is devoted to analyzing the statistical properties
of a power law fluid in heterogeneous fracture media or porous media. In
particular, analytical solutions for the mean and standard deviation of the
velocity field will be derived in the limit of small amplitudes of
heterogeneities $\sigma_{f}\ll 1$.

###  3.4 Perturbation expansion for a power-law rheology

The objective of this section is to derive the mean properties of the flow
field. The approach is to extend the work of Gelhar & Axness (1983) to power
law fluids, which consists in expanding the solution around the mean value and
assuming sufficiently small spatial fluctuations.

Assuming that the flow field is solution of a power-law rheology in
hetergenous porous media:

$-\vec{\nabla}P=c(\vec{r})|\\!|\vec{u}|\\!|^{n-1}\vec{u},$ (48)

with the free divergence:

$\nabla.\vec{u}=0.$ (49)

In order to compare with the truncated rheology, $c(\vec{r})$ should be equal
to $c(\vec{r})=\frac{{{\mu}}}{\kappa(\vec{r}){u_{c}}^{n-1}}$ with $\kappa$ and
${u_{c}}$ distributed according to log-normal as described previously.
However, in order to remain as general as possible, the generic function
$c(\vec{r})$ will be kept as long as possible.

The principle is to decompose each field $\vec{u}$, $P$, and $c(\vec{r})$ into
a mean part and the spatially fluctuating part

$\displaystyle\vec{u}(\vec{r})$ $\displaystyle=$ $\displaystyle
U\vec{e}_{x}+\delta\vec{u}(\vec{r}),$ (50)
$\displaystyle\vec{\nabla}P(\vec{r})$ $\displaystyle=$ $\displaystyle
G\vec{e}_{x}+\vec{\nabla}\delta p(\vec{r}),$ (51) $\displaystyle c(\vec{r})$
$\displaystyle=$ $\displaystyle C_{0}+\delta c(\vec{r}),$ (52)

with $G=\langle\nabla P\rangle$ and $U=\langle u_{x}\rangle$. The mean flow is
assumed to be along the $\vec{e}_{x}$ axis.

Because $c(\vec{r})$ is distributed according to a log-normal, it is more
convenient, and without loss of generality, to introduce:

$c(\vec{r})=\exp{(g_{0}+\delta g)}=C_{0}\exp(\delta g).$

Expanding the different terms in eq. (48) to the second order leads to,

$\displaystyle c$ $\displaystyle=$ $\displaystyle C_{0}(1+\delta
g+\frac{1}{2}\delta g^{2})$ (53) $\displaystyle|\\!|\vec{u}|\\!|^{n-1}$
$\displaystyle=$ $\displaystyle|\\!|U\vec{e}_{x}+\delta\vec{u}|\\!|^{n-1}$
(54) $\displaystyle=$ $\displaystyle U^{n-1}+(n-1)U^{n-2}\delta
u_{x}+\frac{n-1}{2}U^{n-3}(\delta\vec{u}.\delta\vec{u})+\frac{(n-1)(n-3)}{2}U^{n-3}\delta
u_{x}^{2}.$ (55)

Taking the spatial average $\langle.\rangle$ of eq. (48) all the first order
terms vanish by definition, yielding to:

$\displaystyle-G\vec{e}_{x}$ $\displaystyle=$ $\displaystyle
C_{0}U^{n}\vec{e}_{x}+\frac{1}{2}C_{0}\langle\delta g^{2}\rangle
U^{n}\vec{e}_{x}+C_{0}\frac{n-1}{2}U^{n-2}\langle\delta\vec{u}.\delta\vec{u}\rangle\vec{e}_{x}+\frac{(n-1)(n-3)}{2}U^{n-2}\langle\delta
u^{2}\rangle$ $\displaystyle+C_{0}U^{n-1}\langle\delta
g\;\delta\vec{u}\rangle+C_{0}(n-1)U^{n-1}\langle\delta g\;\delta
u_{x}\rangle\;\vec{e}_{x}+C_{0}(n-1)U^{n-2}\langle\delta
u_{x}\;\delta\vec{u}\rangle.$

Along the $e_{x}$ axis, it follows:

$-\frac{G}{C_{0}U^{n}}=1+\frac{1}{2}\langle\delta
g^{2}\rangle+n\frac{1}{U}\langle\delta g\;\delta
u_{x}\rangle+\frac{n(n-1)}{2}\frac{1}{U^{2}}\langle\delta u_{x}\delta
u_{x}\rangle+\frac{n-1}{2}\frac{1}{U^{2}}\langle\delta u_{y}^{2}\rangle.$ (56)

This expression relates the mean gradient to the mean velocity of the flow if
we know the different cross-correlation terms of the spatially fluctuating
fields, which are determined next.

The first order of eq. (48) gives:

$-\vec{\nabla}\delta p=C_{0}U^{n}\delta g\;\vec{e}_{x}+C_{0}(n-1)U^{n-1}\delta
u_{x}\vec{e}_{x}+C_{0}U^{n-1}\delta\vec{u}.$ (57)

Thus,

$-\frac{1}{C_{0}U^{n-1}}\vec{\nabla}{\delta p}=U\delta
g\vec{e}_{x}+(n-1)\delta u_{x}\vec{e}_{x}+\delta\vec{u}.$ (58)

Taking the curl leads to:

$0=\partial_{x}(\delta u_{y})-\partial_{y}(U\delta g+n\delta u_{x}).$ (59)

It is now more convenient to write this equation in the Fourier space.
Defining $\hat{\vec{u}}(k_{x},k_{y})$ and $\hat{g}(k_{x},k_{y})$, respectively
the Fourier transform of $\delta\vec{u}$ and $\delta g$, gives:

$ik_{x}\hat{u}_{y}-ik_{y}(U\hat{g}+n\hat{u}_{x})=0.$ (60)

Using the free divergence in Fourier space:

$k_{x}\hat{u}_{x}+k_{y}\hat{u}_{y}=0,$ (61)

it follows the relationship between the fluctuation of $\delta\vec{u}$ and
$\delta g$, in Fourier space:

$\displaystyle\hat{u}_{x}$ $\displaystyle=$
$\displaystyle-\frac{k_{y}^{2}}{nk_{y}^{2}+k_{x}^{2}}\hat{g}$ (62)
$\displaystyle\hat{u}_{y}$ $\displaystyle=$
$\displaystyle\frac{k_{y}k_{x}}{nk_{y}^{2}+k_{x}^{2}}\hat{g}.$

From these expressions, it is possible to determine the different cross
correlation terms using Parseval’s formula:

$\iint F\;H\;dxdy=\iint\hat{F}\;\hat{H}^{*}\;dk_{x}dk_{y},$

for any field $F(\vec{r})$ and $H(\vec{r})$.

Thus,

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{U^{2}}\langle\delta u_{x}\delta
u_{x}\rangle=\iint\hat{u}_{x}\hat{u}_{x}^{*}\;dk_{x}dk_{y}=\iint\frac{k_{y}^{4}}{(nk_{x}^{2}+k_{y}^{2})^{2}}\;\hat{g}\hat{g}^{*}\;dk_{x}dk_{y}$
(63) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{U}\langle\delta u_{x}\delta
g\rangle=\iint\hat{u}_{x}\hat{g}^{*}\;dk_{x}dk_{y}=-\iint\frac{k_{y}^{2}}{nk_{x}^{2}+k_{y}^{2}}\;\hat{g}\hat{g}^{*}\;dk_{x}dk_{y}$
(64) $\displaystyle I_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{U^{2}}\langle\delta u_{y}\delta
u_{y}\rangle=\iint\hat{u}_{y}\hat{u}_{y}^{*}\;dk_{x}dk_{y}=\iint\frac{k_{y}^{2}\;k_{x}^{2}}{(nk_{x}^{2}+k_{y}^{2})^{2}}\;\hat{g}\hat{g}^{*}\;dk_{x}dk_{y}.$
(65)

These equations are very general and should apply to any distribution and
correlation of permeability field, providing that the amplitude of the
heterogeneities is small enough.

Using now the particular distribution field $g(\vec{r})$ with
$c(\vec{r})=\frac{{{\mu}}}{\kappa(\vec{r}){u_{c}}^{n-1}}$ and
${u_{c}}=A\kappa^{\gamma}$ leads to:

$g_{0}+\delta g=\log{{\mu}}+(1-n)\log A-(1+\gamma(n-1))(f_{0}+f).$ (66)

$\delta g$ is thus also a log-normal distribution with
$\sigma_{g}=\sqrt{\langle\delta g^{2}\rangle}=(1+\gamma(n-1))\sigma_{f}$ and

$C_{0}=\frac{{{\mu}}}{\kappa_{0}u_{0}^{n-1}}={{\mu}}\kappa_{0}^{-1-\gamma(n-1)}A^{1-n}.$
(67)

It follows that $g$ is also correlated with a Gaussian function with:

$\hat{g}\hat{g}^{*}=\beta^{2}e^{-2\frac{k_{x}^{2}+k_{y}^{2}}{k_{0}^{2}}},$
(68)

where $\beta$ is a normalisation prefactor determined by
$\sigma_{g}=\sqrt{\langle\delta g^{2}\rangle}$.

The different correlation functions can be derived after some manipluations,

$\displaystyle I_{1}(n)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{g}^{2}}{2\pi}\int_{0}^{2\pi}\frac{\sin^{4}\theta}{(n\cos^{2}\theta+\sin^{2}\theta)^{2}}\;d\theta=\frac{1}{2}\sigma_{g}^{2}\frac{1+2\sqrt{n}}{(1+\sqrt{n})^{2}n^{3/2}}$
(69) $\displaystyle I_{2}(n)$ $\displaystyle=$
$\displaystyle-\frac{\sigma_{g}^{2}}{2\pi}\int_{0}^{2\pi}\frac{\sin^{2}\theta}{n\cos^{2}\theta+\sin^{2}\theta}\;d\theta=-\sigma_{g}^{2}\frac{1}{\sqrt{n}+n}$
(70) $\displaystyle I_{3}(n)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{g}^{2}}{2\pi}\int_{0}^{2\pi}\frac{\sin^{2}\theta\cos^{2}\theta}{(n\cos^{2}\theta+\sin^{2}\theta)^{2}}\;d\theta=\sigma^{2}_{g}\;\Theta(n).$
(71)

Unfortunately, no analytical form of the last integral could be found. It must
therefore be computed numerically. Note that this integral has the property
$n^{2}I_{3}(n)=I_{3}(1/n)$. Table 1 gives the value of $\Theta(n)$ for the
indexes $n$ used in this manuscript. The most important correlation term is
the standard deviation of the velocity field, Eq. (69), which will be used
later:

$\frac{\sigma^{2}_{u_{x}}}{U^{2}}=\sigma_{f}^{2}\left[1+\gamma(n-1)\right]^{2}\frac{1+2\sqrt{n}}{2(1+\sqrt{n})^{2}n^{3/2}}.$
(72)

n | 1/3 | 1/2 | 1 | 2 | 3
---|---|---|---|---|---
$\Theta(n)$ | 0.348076 | 0.242641 | 0.125 | 0.0606602 | 0.0386751
Table 1:  Numerical values of the function $\Theta(n)$ defined in eq. 71

These expressions can be used in eq. (56) to determine the mean pressure
gradient $G$ by imposing the mean flow rate $U$:

$-\frac{G}{C_{0}U^{n}}=1+\sigma_{f}^{2}\left[1+\gamma(n-1)\right]^{2}\left[\frac{1}{2}-\frac{n}{n+\sqrt{n}}+\frac{n(n-1)}{4}\frac{1+2\sqrt{n}}{(1+\sqrt{n})^{2}n^{3/2}}+\frac{n-1}{2}\Theta(n)\right].$
(73)

A similar procedure could be made to express the mean flow rate $U$ as
function of the imposed gradient of pressure $G$. The full calculation is left
to the reader. The basic idea is to write the constitutive equation in the
form111Another main difference is that in the first order expansion the
divergent has to be taken instead of the curl, to eliminate the velocity and
then relating the fluctuation in pressure as function of the fluctuation in
$C$. $\vec{u}=-D|\\!|\vec{\nabla}P|\\!|^{\alpha-1}\vec{\nabla}P$, with
$\alpha=1/n$ and $D=C^{-\alpha}$. This results to :

$-\frac{U}{D_{0}G^{\alpha}}=1+\sigma_{f}^{2}\left[\alpha+\gamma(1-\alpha)\right]^{2}\left[\frac{1}{2}-\frac{\alpha}{\alpha+\sqrt{\alpha}}+\frac{\alpha(\alpha-1)}{4}\frac{1+2\sqrt{\alpha}}{(1+\sqrt{\alpha})^{2}\alpha^{3/2}}+\frac{\alpha-1}{2}\Theta(\alpha)\right],$
(74)

with
$D_{0}=C_{0}^{-\alpha}={{\mu}}^{-\alpha}\kappa_{0}^{\alpha+\gamma(1-\alpha)}A^{1-\alpha}$.

Surprisingly, this expression is very similar to eq. (73). This similarity
orginates from a symmetry property of 2D flow fields where the role of
pressure and velocity can be switched. More detail will be provided in the
next section. Eqs. (72) and (74) indicate that the mean flow rate is
independant of the disorder amplitude for Newtonian fluid ($n=\alpha=1$) as
demonstrated by Gelhar & Axness (1983). For power fluids, the dependence is
significantly more complicated.

Fig. 5 (left) represents the evolution of the average flow rate as a function
of $n$ for different $\gamma$ as predicted by eq. (74) and using
$D=G=\sigma_{f}=1$. For any $\gamma$, the average flow rate is decreasing with
$n$. For shear-thinning fluids, $n<1$, the curve is greater than one, meaning
that the average flow increases with the heterogeneity $\sigma_{f}$. In
addition, this effect is more pronounced for lower values of $\gamma$. For
shear thickening fluids, $n>1$, the curve is below one, so the average flow
rate decreases with $\sigma_{f}$. This effect is magnified by the higher value
of $\gamma$.

Fig. 5 shows the width of the velocity distribution. For most of the $\gamma$
values, the heterogeneity of the velocity field decreases with $n$. Shear-
thinning fluids have thus a much more heterogeneous velocity field than the
shear thickening one. And this effect is more important for the lowest value
of $\gamma$. Surprisingly, this trend seems to inverse above a certain value
$\gamma\gtrsim 0.9$. A change in monotony is then observed, with shear
thickening fluids being more heterogeneous than shear thinning ones.

Figure 5:  Left : Mean flow as function of the rheological index $n$ for
different $\gamma$ as predicted by eq. (74), using $D=G=\sigma_{f}=1$. Right:
Normalized standart deviation of the velocity field as function of $n$ given
by eq. (72) with $\sigma_{f}=1$.

###  3.5 Rotational symmetry

This argument orginates from Matheron (1967) in a two-dimensional flow field
for Newtonian Darcy’s law. It can however be generalized to non-Newtonian
fluids. This symmetry is also applicable to the 2D pore network model
(Straley, 1977; Talon & Hansen, 2020). As previously, a generic non-linear
Darcy equation is assumed:

$\vec{\nabla}.\vec{u}=0,$ (75)

and

$\vec{u}=-\frac{f(|\\!|\vec{\nabla}P|\\!|)}{|\\!|\vec{\nabla}P|\\!|}\vec{\nabla}P$
(76)

or

$\vec{\nabla}P=-\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}\vec{u}.$ (77)

The idea is to rotate the two fields $\vec{u}$ and $\vec{\nabla}P$ by
$90^{\circ}$. In a coordinate system $(x,y,z)$, where the flow takes place in
the plane $(x,y)$, this rotation is performed by making the cross product with
the vector $\vec{e}_{z}$:

$\vec{e}_{z}\times\vec{\nabla}P=-\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}(\vec{e}_{z}\times\vec{u}).$
(78)

Defining the rotated fields $\vec{q}=\vec{e}_{z}\times\vec{\nabla}P$ and
$\vec{Z}=\vec{e}_{z}\times\vec{u}$, it can be shown that
$\vec{\nabla}.\vec{q}=0$ and $\vec{\nabla}\times\vec{Z}=\vec{0}$. This means
that $\vec{q}$ is a flux vector and $\vec{Z}$ derives from a potential field
$\vec{Z}=\vec{\nabla}\Psi$. The two new fields satisfy:

$\vec{q}=-g(|\\!|\vec{u}|\\!|)\frac{\vec{\nabla}\Psi}{|\\!|\vec{\nabla}\Psi|\\!|}.$
(79)

Since $g(|\\!|\vec{u}|\\!|)=|\\!|\vec{\nabla}P|\\!|$, it follows
$|\\!|\vec{\nabla}\Psi|\\!|=g^{-1}(|\\!|\vec{q}|\\!|)$. Yielding to:

$\vec{\nabla}\Psi=-g^{-1}(|\\!|\vec{q}|\\!|)\frac{\vec{q}}{|\\!|\vec{q}|\\!|}.$
(80)

As a result, the rotated fields $\vec{q}$ and $\vec{\nabla}\Psi$ satisfy a
non-Newtonian Darcy’s equation but with a rheology inverse to the original
one. In particular, solving a shear-thinning fluid in one direction is then
equivalent to solving a shear thickening in the other direction.

Considering the truncated model:

$\left\\{\begin{array}[]{lll}\vec{\nabla}P=-\frac{{{\mu}}}{\kappa}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|<{u_{c}}\\\
\vec{\nabla}P=-\frac{{{\mu}}}{\kappa}\left[\frac{|\\!|\vec{u}|\\!|}{{u_{c}}}\right]^{n-1}\vec{u}&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|>{u_{c}}\\\ \end{array}\right.,$ (81)

is thus equivalent to:

$\left\\{\begin{array}[]{lll}\vec{\nabla}\Psi=-\frac{\kappa}{\mu}\vec{q}&\;\;\;{\rm
if}&|\\!|\vec{q}|\\!|<q_{c}\\\
\vec{\nabla}\Psi=-\frac{\kappa}{\mu}\left[\frac{|\\!|\vec{q}|\\!|}{q_{c}}\right]^{1/n-1}\vec{q}&\;\;\;{\rm
if}&|\\!|\vec{q}|\\!|>q_{c}\\\ \end{array}\right.,$ (82)

with $q_{c}=\frac{\mu\;{u_{c}}}{\kappa}$. Using the relation
${u_{c}}=A\kappa^{\gamma}$, the parameters have thus changed to:

$\displaystyle\tilde{n}$ $\displaystyle\rightarrow$ $\displaystyle 1/n$
$\displaystyle\tilde{\mu}$ $\displaystyle\rightarrow$ $\displaystyle 1/\mu$
$\displaystyle\tilde{\kappa}$ $\displaystyle\rightarrow$ $\displaystyle
1/\kappa$ $\displaystyle\tilde{\gamma}$ $\displaystyle\rightarrow$
$\displaystyle-\gamma+1$ $\displaystyle\tilde{A}$ $\displaystyle\rightarrow$
$\displaystyle\mu A.$

It is remarkable that for the most natural value $\gamma=1/2$ in porous media,
this coefficient is invariant with this transformation. It should also be
noted that the inverse of a lognormal distribution remains lognormal with
$\tilde{f}_{0}\rightarrow-f_{0}$ and the same $\sigma_{f}$. It follow that the
study can be limited to shear thinning fluids ($n<1$) without loss of
generality. Figure 6 illustrates this symmetry by comparing the average flow
rate as a function of the average pressure gradient for the two exponents
$n=1/2$ and $n=2$. A very good correspondance between the two curves can be
observed.

Figure 6:  Blue curve: mean flow rate $\langle u\rangle$ as function of the
mean gradient of pressure $\langle\nabla P\rangle$ for a truncated rheology
with $n=1/2$, $\gamma=1/2$ and $\sigma_{f}=1$. Red circles: mean gradient of
pressure as function of the mean velocity with $n=2$, $\gamma=1/2$ and
$\sigma_{f}=1$.

##  4 Numerical results

###  4.1 Flow rate

The following section is dedicated to the numerical results obtain for the
truncated theology. We will first analyse the two asymptotic limits and the
study the intermediate flow regime.

####  4.1.1 Influence of the disorder

Fig. 7.a diplays the mean flow-rate as a function of the mean pressure
gradient $\vec{\nabla}P$ for a sheared-thinning fluid, $n=1/2$, and for the
different amplitude of heterogeneities $\sigma_{f}$. As expected, and
similarly to the homogeneous permeability medium, a transient behavior is
observed around $u\sim{u_{0}}$ and $\nabla P\sim{g_{0}}$. At low pressure
drop, the flow rate follows a linear relation ${\langle
u\rangle}\propto\vec{\nabla}P$, while at high pressure difference, it follows
a power law ${\langle u\rangle}\propto\vec{\nabla}P^{1/n}$. The influence of
disorder is manifold. First, it smoothes out the transition and increases its
range. For $\sigma_{f}=0$, the transition is sharp, while for high
$\sigma_{f}$, the transition is smoother, extending over half a decade. This
effect is more pronounced after normalizing the average flow-rate by the
Newtonian Darcy’s law $\frac{\mu{\langle u\rangle}}{\kappa_{0}\nabla P}$ (Fig.
7.b). The second effect of heterogeneities is to increase the flow-rate in
both asymptotic regimes. While it seems very marginal in the low pressure
regime, the effect is quite significant in the high pressure regime, where the
flow rate is almost doubled between $\sigma_{f}=0$ and $\sigma_{f}=2.5$.

Figure 7:  For a truncated shear thinning fluids with $n=1/2$ and
$\gamma=1/2$: non-dimensional mean flow rate (a), mean flow rate normalized by
the Newtonian Darcy’s regime (b) and relative standart deviation of the
velocity field (c) as function of the non-dimensional mean gradient of
pressure for different amplitude of heterogeneities $\sigma_{f}$. (d) Relative
standart deviation of the velocity field for a shear-thickening fluid ($n=2$).

Figures 7.c and 7.c represent the relative standard deviation of the flow
field $\frac{\sigma_{ux}}{{\langle u\rangle}}$ as a function of the applied
pressure difference and for different $\sigma_{f}$. It displays a change of
plateau when changing the flow regime. The presence of the two plateaus can be
easily understood from the results of the section 3.1, where it has been shown
that the $\frac{\vec{u}}{{\langle u\rangle}}$ field tends towards a constant
in low and high flow regimes. Here again, the extension of the transition
between these two asymptotes depends on the heterogeneities of the porous
medium, sharp at low $\sigma_{f}$ while smoother at high $\sigma_{f}$. For
shear thinning fluids, flow heterogeneity is increased in the high flow rate
regime, while it is decreased for shear thickening fluids. In addition, flow
heterogeneity increases with the disorder of the medium $\sigma_{f}$. For
shear thining fluids, this increase is more important in the power-law regime
than in the linear one. Compared to the Newtonian regime the velocity field
can either be more heterogenous for shear-thinning or less heterogenous for
shear-thickening.

As described above, the flow field in both asymptotic regimes converges to the
solution of a simple power law fluid whose characteristics were predicted at
low $\sigma_{f}$ in the section 3.4. Fig. 8 compares the numerical mean and
standard deviation of the flow field at sufficiently low (resp. high) applied
pressure gradient $\nabla P/{g_{0}}=1e-2$ (resp. $\nabla P/{g_{0}}=100$)
against the prediction of eqs. (72) and (74). The two figures show a very good
agreement between the analytical prediction and the simulations, even for an
amplitude of heterogeneity as high as $\sigma_{f}\sim 2$, which is quite
significant for a log-normal distribution. An expansion of order 2 seems
therefore sufficient to predict both the mean flow and the standard deviation
of the flow field as the variation seems to remain parabolic over a wide range
of $\sigma_{f}$.

Figure 8:  For a shear thinning fluid $n=1/2$ and $\gamma=1/2$. Left: mean
velocity versus the amplitude of the heterogeneities $\sigma_{f}$ at low
applied pressure drop $\nabla P/{g_{0}}=1e-2$ (blue) and high pessure drop
$\nabla P/{g_{0}}=100$ (red). Right: relative standard deviation as function
of $\sigma_{f}$. Circles represent the simulations and the dashed line is the
perturbation expansion prediction eqs. (72) and (74).

####  4.1.2 Range of the transitional flow regime

As described above, this range of the transitional flow regime depends
significantly on $\sigma_{f}$. To understand this transition, it is useful to
examine the evolution of the velocity distribution with the mean pressure
gradient shown in Fig. 9.a. For a small pressure gradient, as long as the
distribution is away from the mean critical velocity ${u_{0}}$, its shape
remains constant while being shifted as $\nabla P$ increases.

As the distribution approaches the mean critical velocity ${u_{0}}$, it begins
to deform. This reflects the fact that some regions begin to reach their
critical velocity and thus changing their viscosity. At sufficiently high
pressure, once this distribution is far enough from ${u_{0}}$, the
distribution adopts another constant shape which is only shifted when
increasing $\nabla P$. In Fig.9.b the two asymptotic forms
${{}_{1}\vec{\omega}}$ and ${{}_{1/2}\vec{\omega}}(\vec{r})$ are plotted.
Notably, the shear thinning fluid $n=1/2$ has a wider velocity distribution
but is also slightly more asymmetrical towards the low velocity.

The invariance of the ${{}_{n}\vec{\omega}}$ distribution can be exploited to
estimate the range of the transition regime. Starting at a very low flow rate,
all regions are in Newtonian regime and $|\\!|\vec{u}|\\!|={\langle
u\rangle}|\\!|{{}_{1}\vec{\omega}}|\\!|$. Increasing the average flow rate,
the first location to change its rheology occurs where
$|\\!|\vec{u}|\\!|={u_{c}}$. Thus:

$\max{\frac{{\langle
u\rangle}\;|\\!|{{}_{1}\vec{\omega}}(\vec{r})|\\!|}{{u_{c}}(\vec{r})}}=1,$
(83)

which means that the change of regime will start at the average velocity
${\langle u\rangle}=u_{1}$ satisfying:

$\frac{1}{u_{1}}=\max{\frac{|\\!|{{}_{1}\vec{\omega}}(\vec{r})|\\!|}{{u_{c}}(\vec{r})}}.$
(84)

Similarly, starting at high pressure drop, all the regions are in the power-
law regime and $\vec{u}(\vec{r})/{\langle
u\rangle}={{}_{n}\vec{\omega}}(\vec{r})$ is constant. While decreasing
${\langle u\rangle}$, the first point to change its behavior occurs at:

$\frac{1}{u_{2}}=\min{\frac{|\\!|{{}_{n}\vec{\omega}}(\vec{r})|\\!|}{{u_{c}}(\vec{r})}}.$
(85)

Assuming that at order zero $\nabla
P_{1}\simeq\frac{{{\mu}}}{\kappa_{0}}u_{1}$ and $\nabla
P_{2}\simeq{{\mu}}{\kappa_{0}}^{-1-\gamma(n-1)}A^{1-n}u_{2}^{n}$, it leads to
an estimation for the bounds of the transient regime:

Lower bound:

$\nabla
P_{1}=\frac{{{\mu}}}{\kappa_{0}\max{\frac{{{}_{1}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}}}$
(86)

Upper bound:

$\nabla
P_{2}=\frac{{{\mu}}{\kappa_{0}}^{-1-\gamma(n-1)}A^{1-n}}{\left(\min{\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}}\right)^{n}}.$
(87)

The range of the transition hinges thus on the amplitude of the
heterogeneities of the velocity field. If ${{}_{1}\vec{\omega}}$ is broad, it
is expected that the transition will start at lower pressure gradient because
the maximum is larger. And if ${{}_{n}\vec{\omega}}$ is broad, the transition
should end at higher pressure gradient. It is important to emphasize, however,
that the relevant quantity is in fact the width of
$\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}$ which does not
necessarily follow the width of ${{}_{n}\vec{\omega}}$, as will be discussed
below.

It can also be noted that, since the extension of the transitional regime
depends on the extrema of ${{}_{n}\vec{\omega}}/{u_{c}}$ and
${{}_{1}\vec{\omega}}/{u_{c}}$, this range is expected to depend on the width
of these two distributions but also on the size of the system. Indeed, for a
field of size $L_{x}L_{y}$ with a correlation length $\lambda$, it can
approximately be considered as $N=L_{x}L_{y}/\lambda^{2}$ independent random
numbers. On average, theory of extreme statistics suggests that:

$\langle\max_{\vec{r}}\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}\rangle=\mathbb{P}^{-1}_{n}(1-1/N)$
(88)

and

$\langle\min_{\vec{r}}\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}\rangle=\mathbb{P}^{-1}_{n}(1/N),$
(89)

where $\mathbb{P}_{n}$ is the cumulative distribution function of
${{}_{n}\vec{\omega}}/{u_{c}}$. The range of the transitional regime $[\nabla
P_{1};\nabla P_{2}]$ is thus expected to increase with the size of the system,
which seems quite intuitive. The larger the size of the system, the higher the
maximum will be and the lower the minimum. The exact dependence of these
extrema with the system size depends on the specific form of the distribution,
particularly on the tails.

Fig. 10 represents different phase diagrams of the flow regime as a function
of the mean pressure gradient and varying $\sigma_{f}$, $\gamma$ or $n$. A
cross represents system which are in the intermediate regime. In addition, the
bounds predicted by eqs. (86) and (87) have been plotted and they show good
agreement.

Fig. 10.a represents the evolution of the transient regime as a function of
the amplitude of the disorder $\sigma_{f}$. As expected, the amplitude of the
transition increases with the disorder. The upper bound $\nabla P_{2}$
increases however faster than the lower one $\nabla P_{1}$, this reflects the
fact that the width of the distribution of ${{}_{1/2}\vec{\omega}}$ is larger
than that of ${{}_{1}\vec{\omega}}$.

Fig. 10.b shows the dependence of the extension with the rheological index
$n$. Albeit the width of the ${{}_{n}\vec{\omega}}$ distribution is decreasing
with $n$ (see figure 11), the range increases with $n$. This inversion is in
fact due to the presence of the exponent $n$ in eq. (87) which inverses the
expected trend.

The influence of $\gamma$ (Fig. 10.c) on the transition range is more
surprising as it exhibits a non-monotonic behavior. For low $\gamma$ values,
the range is decreasing with $\gamma$ but it becomes increasing beyond a
certain $\gamma$ value. Here also, the transition doesn’t follow necessarily
the amplitude of the velocity distribution which is decreasing with $\gamma$
(see Fig. 11 or the perturbation expansion results). This shows that the
standard deviation of
$\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}=\frac{{{}_{n}\vec{\omega}}(\vec{r})}{A\kappa^{\gamma}(\vec{r})}$
does not necessarily follow the one of ${{}_{n}\vec{\omega}}(\vec{r})$. This
is due to the fact that the velocity and permeability field are not entirely
correlated: higher velocity regions do not necessarily correspond to high
permeability (${u_{c}}$) ones.

This can be understood by considering the correlations between the two fields.
While the permeability field is anisotropic, the velocity field is much more
correlated in streamwise direction due to the conservation of mass. A region
of low permeability for instance can thus have a high velocity if it is
located in front or behind a region of high permeability. The correlation of
the velocity field is thus a crucial characteristic to understand the
evolution of the transition regime. This point will be discussed in more
details in the next section.

Figure 9:  For a shear-thinning fluid, $n=1/2$ and $\gamma=1/2$. Left:
Evolution of the probability distribution function (PDF) of the velocity as
function of the applied mean pressure gradient. The vertical dashed line
represents the mean threshold velocity ${u_{0}}$. Right: normalized velocity
distribution ${{}_{1}\vec{\omega}}$ and ${{}_{n}\vec{\omega}}$, obtained
respectively at $\nabla P/{g_{0}}=1e-2$ and $\nabla P/{g_{0}}=1e2$.

Figure 10:  Flow regime diagram of the system as function of the amplitude of
the hetergogenities $\sigma_{f}$ and the mean pressure gradient. Crosses
represent systems in the intermediate flow regime defined as ${\cal
O}\in[10^{-5},1-10^{-5}]$ for different parameters: left
$(n,\gamma)=(1/2,0.5)$, center: $(n,\sigma)=(1/2,1)$, and right
$(\sigma,\gamma)=(1,0.5)$. The red and blue dashed line represent respectively
the bounds predicted by eqs. (86) and (87).

Figure 11:  Standart deviation of the field
$\frac{{{}_{n}\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}$ (orange) and
${{}_{n}\vec{\omega}}(\vec{r})$ (blue) obtain from the simulations as function
of $\gamma$ (left) and $n$ (right). Left: $(\sigma;n)=(1,0.5)$; right:
$(\sigma;\gamma)=(1,0.5)$;

###  4.2 Statistical propertie of the flow field in the transitional regime

In this section, the properties of the flow field within the transitional
regime are analysed, which diplays some interesting statistical properties.
One way to apprehend it is to notice a similarity with the percolation
problem. Indeed, as shown in figure 12, when ${\langle u\rangle}$ is
increased, more and more regions meet the criteria
$|\\!|\vec{u}|\\!|>{u_{c}}$. This also allows to define cluster of connected
regions. If the field $|\\!|\vec{u}|\\!|/{\langle
u\rangle}={\vec{\omega}}(\vec{r})$ were constant, the transition would occur
at

$\frac{{\vec{\omega}}(\vec{r})}{{u_{c}}(\vec{r})}=\frac{1}{{\langle
u\rangle}},$

which would correspond to a percolation problem. It is thus expected to
observe similar behaviors such as fractal properties of these clusters and the
presence of criticality. It is however important to remember that the problem
is not strictly equivalent to percolation because the $\frac{\vec{u}}{{\langle
u\rangle}}$ field is not constant in the transient regime. In particular, the
change in viscosity introduces correlations in the velocity field and thus
modifies the shape of the clusters as described previously.

Figure 12:  Cluster of regions (in black) above their threshold
$|\\!|\vec{u}|\\!|>{u_{c}}$ for a shear-thinning fluid: $n=0.5$ and
$\gamma=0.5$, $\sigma=1$ and for different applied pressure drop.

The shape of the clusters are characterized by their total size $S$ and the
two dimensions of the bounding rectangle containing it: $L$ along the flow
direction and $W$ transversely to it.

####  4.2.1 Size distribution

Figure 13:  Left: Probability distribution function of the cluster size for
different applied pressure drops using the parameters
$(n,\gamma,\sigma_{e})=(1/2,1/2,1)$. Dashed black line represents the fitted
power-law exponent. Right : Size of the largest cluster $S_{0}$ versus the
pressure gradient. The dotted vertical line represents $\nabla P_{c}$. Inset:
$S_{0}$ as a function of $\nabla P-\nabla P_{c}$ on a logarithmic scale, the
line represents the fitted power law. Bottom: Size probability distribution
function for different pressure gradients normalized according to eqs. (90)
and (91).

Fig 13.a displays the size distribution $P(S)$ for different applied pressure
drops and for the parameters $n=1/2$ and $\sigma=1$. For any applied pressure,
the distribution follows a decaying power-law over a large range of sizes.
There is, however, a large-scale cut-off, $S_{0}$, which varies with the
pressure gradient. Reporting the variation of $S_{0}$ as a function of the
pressure (Fig. 13.b) shows that $S_{0}$ diverges at a certain value of the
mean pressure gradient $\nabla P=\nabla P_{c}$ according to a power law :

$S_{0}\propto|\nabla P-\nabla P_{c}|^{-\nu_{S}}.$ (90)

Combining the two observations, it follows the scaling law for the size
distribution:

$p(S)\propto S^{-\tau_{S}}f(\frac{S}{S_{0}}),$ (91)

which is confirmed by the good collapse of the distributions rescaled
according to eqs. (90) and (91) plotted in Fig. 13.c.

This scaling law is thus similar to the one found in other problems with a
critical transition such as percolation (Stauffer & Aharony, 1991) or
avalanches (e.g Amaral et al. (1995); Santucci et al. (2011); Chevalier et al.
(2017)). The sizes are distributed on many scales up to a size limit. And this
size limit diverges as the control parameter $\nabla P$, approaches the
critical value $\nabla P=\nabla P_{c}$. The system is then scale invariant
(_e.g._ fractal). These scaling laws are characterized by the two exponents
$\tau_{S}$ and $\nu_{S}$. It is important to note that a similar scaling law
can be observed for the length $L$ of the clusters:

$\displaystyle p(S)\propto S^{-\tau_{L}}f(\frac{L}{L_{0}})$
$\displaystyle\rm{with}$ $\displaystyle L_{0}\propto|\nabla P-\nabla
P_{c}|^{-\nu_{L}},$ (92)

allowing to define the exponents $\tau_{L}$ and $\nu_{L}$.

The measured exponents $\tau_{S}$, $\tau_{L}$, $\nu_{S}$ and $\nu_{L}$ for the
different sets of parameters are reported in tables 2-4. The most remarkable
result is the fact that $\tau_{S}$ seems to be indeed independent of the
parameters.

It is important to recall that in many critical systems, some exponents are
independant of the details of the disorder distribution. This behavior is
often refered as “universal”. This is the case for instance for the exponents
$\tau$ and $\nu$ in percolation or avalanches of an elastic line in a random
media (Barabasi & Stanley, 1995). The results here seem then to suggest an
universal behavior for the exponent $\tau_{S}$. However, the value
$\tau_{S}=1.65\pm 0.05$ is quite different from the standard percolation
problem ($\tau_{perc}=2.05$), which indicates that it would be of a different
universality class.

The trend is less clear for the exponent $\nu_{S}$, which seems to vary with
the rheological index $n$ and also the heterogeneities $\sigma_{f}$. However,
it is worth noting that the determination of the exponent $\nu_{S}$ is
generally more prone to errors as it requires the determination of $\nabla
P_{c}$, which is also subject to uncertainties. The error can be estimated to
be about $10\%$. It is then difficult to conclude on the universality of this
exponent.

####  4.2.2 Cluster’s shape

The shape of the clusters can also be characterized by the aspect ratio $L/W$.
Fig. 14.a represents the width of each cluster according to their length $L$
for two different parameter sets. A power-law relationship is observed:

$\langle W\rangle_{L}\propto L^{\zeta},$ (93)

which is a characteristic of the self-affine structure generally found in
anisotropic critical systems (e.g. avalanches, directed percolation, front
propagations, etc.). It characterizes the fact that, if many different cluster
sizes are present, the aspect ratio is not the same at each scale. If
$\zeta<1$, bigger clusters are more elongated than the smaller ones as
represented in Fig. 14.a (with $\gamma=0.25$). For $\zeta>1$, as for
$\gamma=0.75$ in this figure, bigger clusters are elongated in the direction
transverse to the flow.

Because clusters are not compact, viz. they may contain holes, another
interesting quantity to analyze is the area (mass) as a function of the
enclosing box size $WL$, as shown in figure 14.b. Here, again the relationship
observed is a power-law:

$S\propto(LW)^{\beta},$

which is a characteristic of fractal structure.

The measured exponents for different parameter sets, $\sigma$, $n$ and
$\gamma$, are displayed in table 2-4. The size exponent $\beta$ appears to be
almost constant, within the error bar, for any parameter value: $\beta\simeq
0.77\pm 0.05$. This exponent seems to be universal. More surprising is the
evolution of the aspect ratio exponent $\zeta$ which takes only two values:
either $\zeta\simeq 0.85\pm 0.05$ or $\zeta\simeq 1.15\pm 0.05$ depending on
the parameters. The two cases display in Figure 15 are in fact the only two
observable exponents. In this figure, the two cases appear very similar if
they were rotated by $90^{\circ}$. This is reminiscent to the rotational
symmetry described above. Indeed, concerning the influence of the rheological
index $n$, this change of correlation direction can be apprehended from the
rotational symmetry argument described in section 3.5.

Indeed, inversing the rheological index is equivalent to a rotation by
$90^{\circ}$. One can thus expect that, by inversing $n\rightarrow 1/n$, the
shape of the clusters are remaining identical providing that we inverse the
role of $W$ and $L$. This argument leads to the relation:

$\zeta(n)=\frac{1}{\zeta(1/n)},$

which seems to be satisfied by the two observed exponents, so the symmetry is
respected. The value $\zeta(n=1)\simeq 1$ appears thus to be a marginal value.
It remains quite remarkable that this exponent is constant for any shear-
thinning (or shear-thickening) $n$ exponent.

More unexpected is the similar change when $\gamma$ is varied (table 4).
$\zeta$ is constant up to a certain value $\gamma\sim 0.6$, where it switches
to its inverse value. Some conjecture can be made to interpret this switching.

First of all, this change of correlation confirms the fact that the spatial
correlation of the velocity field, correlated in the direction of the current,
doesn’t necessarily follow the correlation of the
$\vec{u}(\vec{r})/{u_{c}}(\vec{r})$ field. The value of $\gamma$ at which the
change occurs, seems to correspond to the change in trend observed in Fig. 11
for the evolution of the standard deviation $\sigma_{\omega/{u_{c}}}$.

To understand this change, it is usefull to consider the two extreme cases
$\gamma=0$ and $\gamma=1$. It should be noted that the criterion
$|\\!|\vec{u}|\\!|>{u_{c}}(\vec{r})$ also corresponds to a criterion for local
pressure gradient $|\\!|\vec{\nabla}p(\vec{r})|\\!|>{g_{c}}(\vec{r})$, with
${u_{c}}=A\kappa^{\gamma}$ and ${g_{c}}=\mu A\kappa^{\gamma-1}$. If
$\gamma=0$, ${u_{c}}$ is a constant: it is then expected that the cluster
defined by $|\\!|\vec{u}|\\!|>uc$, follows the correlation of the velocity
field. The cluster should then be more elongated in the streamwise direction.

But for $\gamma=1$, the situation is different because ${u_{c}}(\vec{r})$ is
now a random variable but ${g_{c}}$ is a constant. The clusters, also defined
by $|\\!|\vec{\nabla}p(\vec{r})|\\!|>{g_{c}}$, should have a correlation
direction similar to that of $\vec{\nabla}p$. If the velocity field is
naturally correlated in the flow direction, the pressure gradient tends to be
correlated in the transverse direction222This can be seen from the rotational
symmetry or using the expansion method for power-law fluids.. It is then
expected that the clusters, $|\\!|\vec{\nabla}p(\vec{r})|\\!|>{g_{c}}$, have a
correlation direction similar to that of $\vec{\nabla}p$, transverse to the
flow. For intermediate $\gamma$, there is thus a balance between these two
effects. The argument remains qualitative because, as discussed earlier,
shear-thinning fluids should promote slighlty the correlation in the
streamwise direction.

If the cluster elongation changes with $\gamma$, it is still quite surprising
that the correlation exponent $\zeta$ switches from one value to another. This
seems to suggest that this exponent is universal but a peferential direction
of correlation is determined by the value of $\gamma$.

$\sigma_{f}$ | $0.5$ | $1.0$ | $2.0$
---|---|---|---
$\tau_{S}$ | 1.65 | 1.65 | 1.67
$\nu_{S}$ | 3.0 | 3.25 | 3.2
$\nu_{L}$ | 1.8 | 1.8 | 1.75
$\zeta$ | 0.85 | 0.85 | 0.85
$\beta$ | 0.77 | 0.77 | 0.77
Table 2:  Measured exponents for different amplitude of heterogeneity $\sigma_{f}$ with $(n,\gamma)=(0.5,0.5)$. $n$ | $1/3$ | $0.5$ | $2/3$ | $1$ | $2$ | $3$
---|---|---|---|---|---|---
$\tau_{S}$ | 1.65 | 1.65 | 1.65 | 1.65 | 1.65 | 1.65
$\nu_{S}$ | 3.3 | 3.25 | 3.5 | 3.6 | 3.7 | 3.9
$\nu_{L}$ | 1.95 | 1.8 | 1.8 | 1.8 | 1.9 | 1.9
$\zeta$ | 0.83 | 0.85 | 0.85 | 0.95 | 1.15 | 1.17
$\beta$ | 0.77 | 0.77 | 0.8 | 0.8 | 0.78 | 0.78
Table 3:  Measured exponents for different amplitude of rheological index $n$ with $(\sigma_{f},\gamma)=(1,0.5)$.  $\gamma$ | $0$ | $1/4$ | $1/2$ | $2/3$ | $3/4$ | $1$
---|---|---|---|---|---|---
$\tau_{S}$ | 1.65 | 1.65 | 1.65 | 1.65 | 1.65 | 1.65
$\nu_{S}$ | 3.35 | 3.3 | 3.25 | 3.15 | 3.3 | 3.4
$\nu_{L}$ | 2.1 | 2.1 | 1.8 | 1.9 | 1.9 | 1.9
$\zeta$ | 0.85 | 0.85 | 0.85 | 1.15 | 1.15 | 1.13
$\beta$ | 0.77 | 0.77 | 0.77 | 0.77 | 0.77 | 0.7
Table 4:  Measured exponents for different $\gamma$ with
$(n,\sigma_{f})=(0.5,1)$.

Figure 14:  Left: Width ($W$) of clusters according to their length ($L$).
Right: surface of the cluster $S$ according to the bounding size $WL$. Blue:
$(n,\gamma,\sigma_{f})=(0.5,0.75,1)$ and green:
$(n,\gamma,\sigma_{f})=(0.5,0.25,1)$. The red lines correspond to the power-
law fit.

Figure 15:  Regions above the threshold $u>{u_{c}}$ (in black) for two
different sets of parameters and close to the critical point. The mean flow
direction is from top to bottom. Left: $(n,\gamma,\sigma_{f})=(0.5,0.25,1)$
and right: $(n,\gamma,\sigma_{f})=(0.5,0.75,1)$. The $\gamma$ parameter
drastically affects the orientation of the correlation from the direction of
the flow to perpendicular to it.

##  5 Conclusion

In this paper, we have investigated the flow in macroscopic heterogeneous
porous media with a non-linear rheology displaying a change of behavior such
as the truncated rheology. First, we have discussed the method to derive an
extension of the Darcy’s law for such rheology. The classical approach, which
consists in introducing an effective viscosity or effective stress as function
of an effective shear rate, has been discussed. In particular, it has been
shown that, if for power-law fluid the approach can be correct, it becomes
less accurate for more complicated rheological laws. Using the energy
balance/virtual work approach, it has been shown that the method predicts
qualitatively the asymptotic behaviors. However, one should expect a different
definition of the effective shear rate at low and high flow rate.

The influence of a macroscopic heterogeneous permeability and velocity
threshold fields has then been considered. In this case, an energy
balance/virtual work approach was able to predict some mathematical properties
of the flow field. It shows also the existence of two asymptotic behaviors, at
high and low flow rate. In both limits, the normalized velocity field becomes
constant and corresponds to the solution of a power-law fluid with the
dominant exponent.

The flow field properties of such power-law fluid in a heterogeneous porous
media was studied analytically using a perturbation approach. This allows us
to predict the mean and the standard deviation of the velocity field in both
asymptotic regimes. These predictions were found to be in good agreement with
the simulations. Another usefull theoretical feature is the rotational
symmetry of a 2D porous media, which shows that solving a shear thinning fluid
in one direction is similar to solving a shear thickening in the other
direction. This symmetry property is very interesting in this problem because
it allowed to understand the change of correlation observed in the
simulations.

Then the numerical simulations were presented. First, it shows that two
asymptotic behaviors can be well understood from the power-low flow field with
the dominant exponent. The range of the transitional regime was investigated.
Qualitatively, it is expected that the transition range is related to the
width of the velocity distribution of the power-law fluids in the same
permeability field. For wider velocity distribution should begin the
transition at a lower mean velocity and finish the transition at a higher
velocity. If this trend has been observed, the results have shown that it is
more complicated to make quantitative predictions. Indeed, one very important
aspect of the problem is the relationship between the velocity, the
permeability and the local threshold.

In this paper, a power-law dependency between the permeability and the
threshold have been assumed. The study revealed that the exponent $\gamma$
modifies drastically the correlation of the clusters. The reason behind this
observation is that there are different mechanisms at work. First, high
permeability regions have a higher velocity field but also a higher velocity
threshold. If the velocity threshold vary weakly with the permeability, the
correlation should follow the one of the velocity field, in the streamwise
direction. Conversly, if the velocity threshold varies greatly with the
permeability (_i.e._ $\gamma$ high), the transition might not corresponds to
the highest velocities. At the same time, a similar reasoning considering the
gradient of pressure field can be made: low permeability regions are more
likely to have a higher pressure gradient but also a lower pressure threshold.
So the correlations might be dominated by the pressure gradient field,
transversely to the flow direction. As a result, there is thus a competition
between this two effects which is balanced by the exponent $\gamma$. The
relationship between the permeability and the two thresholds seems then to be
very important. It was assumed a power-law but it is important to stress that
the dependence is probably more complicated than this.

A remark can be made. In the literature, a very common approach is to model a
porous media or a fracture by a bundle of parralel tubes or layers (Di
Federico, 1998; Chen et al., 2005; Nash & Rees, 2016; Rodríguez de Castro &
Radilla, 2017; Felisa et al., 2018). It is then important to note that, in
such models, the flow field is, by construction infinitely correlated in the
streamwise direction, while the pressure gradient is not correlated in the
transverse direction. These models are thus expected to be unable to capture
behaviors related to the change of correlation. They must consequently be
taken with caution while applied to 2D or 3D medium.

The statistical properties of the flow field in the transitional regime
appeared to be quite rich also. Indeed, in this regime, the regions above
their threshold define clusters which display properties found in critical
systems. The cluster sizes are distributed according to a power-law with a
cut-off. And the latter diverges at a particular critical mean gradient of
pressure. The shape of the cluster also displays self-affine fractality. An
important result is that some of the exponents (size distribution, shape
exponents) do not vary with the parameters and the amplitude of the
heterogeneity, which tends to suggest the presence of a universal behaviors.
One very interesting feature is that the shape exponent $\zeta$ is constant
but only switches to its inverse depending on the value of $\gamma$.

This statistical feature is thus very reminiscent to related problems such as
percolation Stauffer & Aharony (1991) and yield stress fluid in porous media
Kostenko & Talon (2019). The flow structure of yield stress fluid is, however,
quite different because below the threshold there is no flow. Regions above
the threshold are then necessarily channel paths connecting the inlet to the
outlet. The exponents are thus expected to be different. It was found
$\tau_{S}\simeq 1.15$ and $\zeta\simeq 0.75$. The size distribution exponent
is different while the roughness exponent is similar to the present case. More
closely would be the percolation problem (directed or not), but as discussed
in the paper, the problem is expected to be a little bit different because the
correlation of the velocity field evolve with the applied gradient of
pressure. The observed exponent are expected to be different: $\tau_{S}\simeq
2.1$ for percolation and $\tau_{S}\simeq 1.26$ for directed percolation. The
present case seems to be intermediate, the exponent are “universal” they fall
in a different universality class.

There are many interesting directions to pursue this work. One important
question is how this generalize to 3D permeability field. While the rotational
symmetry is no longer applicable in 3D, it remains the important fact that the
velocity field is correlated along the streamwise direction while the pressure
gradient correlates in the two transverse directions. A change of correlation
should thus be expected depending on $\gamma$. Critical behaviors should also
be probably observed but with different exponents. Another interesting study
would be to investigate the problem of transport of species in such system.
Indeed, the dispersion of a tracer hinges on the heterogeneities of the
velocity field and its correlation. It is thus expected to observe a change of
behavior due to the change of rheology. In addition, since the molecular
diffusion coefficient is usually linked to the viscosity of the fluid, the
coefficient is expected to be different when the rheology is below and above
the threshold. Finally, another direction of investigation would be the
similarities with others in porous media displaying similar behaviors. For
instance, two phase flow problem (Tallakstad et al., 2009; Yiotis et al.,
2013; Sinha et al., 2017; Yiotis et al., 2019), emulsion driven in porous (Le
Blay et al., 2020) or erosion of granular bed (Aussillous et al., 2016),
display similar critical behavior with the occurrence of preferential flow
paths dependening on the flow rate. These problems have in common the property
that at a certain critical velocity, the local flow condition is drastically
modified, because bubbles are mobilized or because the grains rearrange with
each other. It would thus be very interesting to investigate the similarities
and differences between these problems.

###  5.1 Acknowledgements. —

This work is supported by ”Investissement d’Avenir” LabEx PALM
(ANR-10-LABX-0039-PALM). I would like to thanks R. Kostenko, A. Hansen and A.
Rosso for fruitfull discussions.

##  Appendix A Numerical method

###  A.1 Augmented Lagragian method

This Augmented Lagragian method has been introduced to solve non-Newtonian
Stokes equation and has been used by many authors (see for instance Glowinski
& Le Tallec (1989); Roquet & Saramito (2003)). In this paper the method was
adapted to solve the non-linear Darcy’s equation:

$\displaystyle\vec{\nabla}P$ $\displaystyle=$
$\displaystyle-g(|\\!|\vec{u}|\\!|)\frac{\vec{u}}{|\\!|u|\\!|},$ (94)
$\displaystyle\vec{\nabla}.\vec{u}$ $\displaystyle=$ $\displaystyle 0.$

As for boundary conditions, pressure is imposed at the inlet and outlet,
$P_{in}$ and $P_{out}$ respectively. Periodic conditions are assumed at the
lateral sides. The mean flow is taken in the $x-$direction. As described in
the paper, the solution of such system of equation is equivalent to finding,
among all admissible velocity fields the functional:

$\Phi[\vec{u}]=\int
G(|\\!|\vec{u}|\\!|)dr^{2}-P_{in}\int_{in}\vec{u}_{x}{dy}+P_{out}\int_{out}u_{x}{dy},$
(95)

with

$G(|\\!|\vec{u}|\\!|)=\int_{0}^{|\\!|\vec{u}|\\!|}g(y)dy.$

The main idea of the augmented Lagragian method is to introduce a secondary
field $\vec{v}$ in order to discouple the non-linear problem
$G(|\\!|\vec{u}|\\!|)$ to the flow equation. The equality $\vec{u}=\vec{v}$ is
then garantied by the introduction of a Lagragian vector field $\vec{\mu}$. An
extra term $\frac{\alpha}{2}(\vec{u}-\vec{v})^{2}$ is also added to enhance
the convergence, where $\alpha$ is a small parameter. Another Lagragian field
$\lambda$ is introduced to impose the free divergence, the problem can thus be
recast into a saddle point determination:

$\min_{\vec{u},\vec{v}}\max_{\lambda,\vec{\mu}}\Psi[\vec{u},\vec{v},\vec{\mu},\lambda],$
(96)

with

$\Psi[\vec{u},\vec{v},\vec{\mu},\lambda]=\int\left[G(|\\!|\vec{u}|\\!|)+\vec{\mu}.(\vec{v}-\vec{u})+\frac{\alpha}{2}(\vec{u}-\vec{v})^{2}-\lambda\vec{\nabla}.\vec{v}\right]dr^{2}-P_{in}\int_{in}v_{x}{dy}+P_{out}\int_{out}v_{x}{dy}.$
(97)

If we now differentiate this functional, we obtain:

$\forall\delta\vec{v},\;\;\;\frac{\delta\Psi}{\delta\vec{v}}.\delta\vec{v}=\int\left[\vec{\mu}+\alpha\vec{v}-\alpha\vec{u}\right].\delta\vec{v}-\int\lambda\vec{\nabla}\delta\vec{v}\;dr^{2}-P_{in}\int_{in}\delta
v_{x}{dy}+P_{out}\int_{out}\delta v_{x}{dy},$ (98)
$\forall\delta\lambda,\;\;\;\frac{\delta\Psi}{\delta\lambda}\delta\lambda=-\int\delta\lambda\;\vec{\nabla}\vec{v}\;dr^{2}$
(99)
$\forall\delta\vec{u},\;\;\;\frac{\delta\Psi}{\delta\vec{u}}.\delta\vec{u}=\int\left[g(|\\!|\vec{u}|\\!|)\frac{\vec{u}}{|\\!|\vec{u}|\\!|}-\vec{\mu}+\alpha\vec{u}-\alpha\vec{v}\right].\delta\vec{u}\;dr^{2}$
(100)
$\forall\delta\vec{\mu},\;\;\;\frac{\delta\Psi}{\delta\vec{\mu}}.\delta\vec{\mu}=\int(\vec{v}-\vec{u}).\delta\vec{\mu}dr^{2}.$
(101)

This set of equations is quite cumbersome. However, the main advantage of this
approach lies in the fact that, for a given $\vec{u}$ and $\vec{\mu}$, finding
the saddle point for $\vec{u}$ and $\lambda$ eqs (98)-(99), is equivalent to
solving:

$\displaystyle\vec{v}=-\frac{1}{\alpha}(\vec{\nabla}\lambda+\vec{\mu})+\vec{u}$
(102) $\displaystyle\vec{\nabla}.\vec{v}=0,$ (103)

with the boundary conditions $\lambda=P_{in}$ and $\lambda=P_{out}$. This is a
standart linear Darcy’s equation with a source term, which can thus be solved
using standart methods. Here, a second order finite difference method was
used. It is important note that in these equation the permeability is
homogeneous and constant, which allows a very fast solution at each step.

For the given field $\vec{v}$, $\lambda$ and $\vec{\mu}$, the minimization of
$\Psi$ with respect to $\vec{u}$, is equivalent to solve:

$(\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}+\alpha)\vec{u}=\vec{\mu}+\alpha\vec{v},$
(104)

which represents an implicit problem. This can be solve numerically or
analytically. With our particular function $g(|\\!|u|\\!|)$, an analytical
solution can be found for some exponents $n=1/4,1/3,1/2,3/4\;\;\rm{and}\;\;1$.
Solutions are given in the next section.

After defining initial fields $\vec{u}_{0}$, $\vec{v}_{0}$, $\vec{\mu}_{0}$
and $\lambda_{0}$, the alghorithm is thus decomposed in the following step:

  * 1.

Solve $\vec{v}_{n+1}$ and $\lambda_{n+1}$ with the Darcy’s equation:

$\displaystyle\vec{v}_{n+1}=-\frac{1}{\alpha}(\vec{\nabla}\lambda_{n+1}+\vec{\mu}_{n})+\vec{u}_{n}$
(105) $\displaystyle\vec{\nabla}.\vec{v}_{n+1}=0,$ (106)

  * 2.

Determine $\vec{u}_{n+1}$, by solving

$(\frac{g(|\\!|\vec{u}_{n+1}|\\!|)}{|\\!|\vec{u}_{n+1}|\\!|}+\alpha)\vec{u}_{n+1}=\vec{\mu}_{n}+\alpha\vec{v}_{n+1},$
(107)

  * 3.

Advancing $\vec{\mu}$ toward the gradient eq. (101):

$\vec{\mu}_{n+1}=\vec{\mu}_{n}+(\vec{v}_{n+1}-\vec{u}_{n+1})d\mu,$ (108)

where $d\mu$ is a small parameter, taken equal to $\alpha$ for simplicity.

###  A.2 Solutions of eq. (104)

For the given fields $\vec{v}$ and $\vec{\mu}$, the alghorithm requires to
find $\vec{u}$ satisfying:

$(\frac{g(|\\!|\vec{u}|\\!|)}{|\\!|\vec{u}|\\!|}+\alpha)\vec{u}=\vec{\mu}+\alpha\vec{v}\equiv\vec{B},$
(109)

with

$\left\\{\begin{array}[]{lll}g(|\\!|\vec{u}|\\!|)=\frac{\mu}{\kappa}|\\!|\vec{u}|\\!|&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|<{u_{c}}\\\
g(|\\!|\vec{u}|\\!|)=\frac{{{\mu}}}{\kappa}\left[\frac{|\\!|\vec{u}|\\!|^{n}}{{u_{c}}^{n-1}}\right]&\;\;\;{\rm
if}&|\\!|\vec{u}|\\!|>{u_{c}}\\\ \end{array}\right.$ (110)

It must be noted that $\vec{u}$ and $\vec{B}$ are colinear and with the same
orientation because the left term in eq. (109) is positive. It is thus
sufficient to determine the norm of $|\\!|\vec{u}|\\!|=u$. The equations then
become:

$\left\\{\begin{array}[]{lll}(\frac{\mu}{\kappa}+\alpha)u=B&\;\;\;{\rm
if}&B<(\frac{\mu}{\kappa}+\alpha){u_{c}}\\\ Au^{n}+\alpha u=B&\;\;\;{\rm
if}&B>(\frac{\mu}{\kappa}+\alpha){u_{c}}\\\ \end{array}\right.,$ (111)

where $A=\frac{{{\mu}}}{\kappa{u_{c}}^{n-1}}$.

If the first equation is trivial, the second one has an analytical solution
only for specific value of $n$. We give here the ones used in this work.

  * •

$n=\frac{1}{3}$:

$u=-\frac{\sqrt[3]{\frac{2}{3}}A^{3}}{\Delta}+\frac{\Delta}{18^{1/3}\alpha^{3}}+\frac{B}{\alpha},$

with

$\Delta=\left(\sqrt{3}\sqrt{4A^{9}\alpha^{9}+27A^{6}\alpha^{10}B^{2}}-9A^{3}\alpha^{5}B\right)^{1/3}$

  * $n=\frac{1}{2}$:

$u=\left(\frac{-A+\sqrt{A^{2}+4\alpha B}}{2A}{}\right)^{2}$

  * •

$n=\frac{2}{3}$:

$\displaystyle\text{if}\;\;B<\frac{4A^{3}}{27\alpha^{2}},$ $\displaystyle
u=\frac{1}{\alpha^{3}}\left[-(A^{3}-3\alpha^{2}B)+2A^{3/2}\sqrt{A^{3}-6\alpha^{2}B}\cos{(\phi/3)}\right],\text{
with }$
$\displaystyle\phi=\arg{\left[-2A^{6}+18A^{3}\alpha^{2}B-27\alpha^{4}B^{2}+i3\sqrt{3}\alpha^{3}\sqrt{4A^{3}B^{3}-27\alpha^{2}B^{4}}\right]}$
$\displaystyle\text{if}\;\;B>\frac{4A^{3}}{27\alpha^{2}},$ $\displaystyle
u=\frac{1}{3\alpha^{3}}\left[-(A^{3}-3\alpha^{2}B)-2^{1/3}(-A^{6}+6A^{3}\alpha^{2}B)\frac{1}{\Delta}+\frac{\Delta}{2^{1/3}}\right],\text{
with }$
$\displaystyle\Delta=-2A^{6}+18A^{3}\alpha^{2}B-27\alpha^{4}B^{2}+3\sqrt{3}\alpha^{3}\sqrt{-4A^{3}B^{3}+27\alpha^{2}B^{4}}$

  * •

$n=2$:

$u=\frac{-\alpha+\sqrt{\alpha^{2}+4AB}}{2A}$

  * •

$n=3$:

$u=\frac{\Delta}{18^{1/3}A}-\frac{\sqrt[3]{\frac{2}{3}}\alpha}{\Delta},$

with

$\Delta=\left(9A^{2}B+\sqrt{3}\sqrt{27A^{4}B^{2}+4A^{3}\alpha^{3}}\right)^{1/3}.$

## References

  * Al-Fariss & Pinder (1987) Al-Fariss, T. & Pinder, K. L. 1987 Flow through porous media of a shear-thinning liquid with yield stress. Can. J. Chem. Eng. 65 (3), 391–405.
  * Amaral et al. (1995) Amaral, L. A. N., Barabasi, A. L., Buldyrev, S. V., Harrington, S. T., Havlin, S., Sadrlahijany, R. & Stanley, H. E. 1995 Avalanches and the directed percolation depinning model - experiments, simulations, and theory. Phys. Rev. E 51 (5), 4655–4673.
  * Andrade Jr et al. (1999) Andrade Jr, JS, Costa, UMS, Almeida, MP, Makse, HA & Stanley, HE 1999 Inertial effects on fluid flow through disordered porous media. Physical Review Letters 82 (26), 5249\.
  * Auradou et al. (2008) Auradou, H., Boschan, A., Chertcoff, R., Gabbanelli, S., Hulin, J.-P. & Ippolito, I. 2008 Enhancement of velocity contrasts by shear-thinning solutions flowing in a rough fracture. Journal of non-newtonian fluid mechanics 153 (1), 53–61.
  * Auriault et al. (2002) Auriault, J.-L., Royer, P. & Geindreau, C. 2002 Filtration law for power-law fluids in anisotropic porous media. International Journal of Engineering Science 40 (10), 1151 – 1163.
  * Aussillous et al. (2016) Aussillous, P., Zou, Z., Guazzelli, É., Yan, L. & Wyart, M. 2016 Scale-free channeling patterns near the onset of erosion of sheared granular beds. Proceedings of the National Academy of Sciences , arXiv: http://www.pnas.org/content/early/2016/10/04/1609023113.full.pdf.
  * Balhoff (2005) Balhoff, M. 2005 Modeling the flow of non-newtonian fluids in packed beds at the pore scale. PhD Thesis .
  * Balhoff & Thompson (2004) Balhoff, M. T. & Thompson, K. E. 2004 Modeling the steady flow of yield-stress fluids in packed beds. AIChE J. 50 (12), 3034–3048.
  * Barabasi & Stanley (1995) Barabasi, A-L & Stanley, H. E. 1995 Fractal concepts in surface growth. Cambridge university press.
  * Barnes et al. (1989) Barnes, H.A., Hutton, J.F. & Walters, K. 1989 An introduction to rheology, , vol. 3. Elsevier Science Limited.
  * Bauer et al. (2019) Bauer, D., Talon, L., Peysson, Y., Ly, H. B., Batôt, G., Chevalier, T. & Fleury, M. 2019 Experimental and numerical determination of darcy’s law for yield stress fluids in porous media. Phys. Rev. Fluids 4, 063301.
  * Bessonov et al. (2016) Bessonov, N., Sequeira, A., Simakov, S., Vassilevskii, Y. & Volpert, V. 2016 Methods of blood flow modelling. Mathematical modelling of natural phenomena 11 (1), 1–25.
  * Bird (1976) Bird, R B 1976 Useful non-newtonian models. Annual Review of Fluid Mechanics 8 (1), 13–34, arXiv: https://doi.org/10.1146/annurev.fl.08.010176.000305.
  * Bird et al. (1987) Bird, R B., Armstrong, R.C & Hassager, O. 1987 Dynamics of polymeric liquids. vol. 1: Fluid mechanics .
  * Boschan et al. (2008) Boschan, A., Ippolito, I., Chertcoff, R., Auradou, H., Talon, L. & Hulin, J. P. 2008 Geometrical and taylor dispersion in a fracture with random obstacles: An experimental study with fluids of different rheologies. Water Resour. Res. 44 (6), W06420.
  * Boyd et al. (2007) Boyd, J., Buick, J. M. & Green, S. 2007 Analysis of the casson and carreau-yasuda non-newtonian blood models in steady and oscillatory flows using the lattice boltzmann method. Physics of Fluids 19 (9), 093103, arXiv: https://doi.org/10.1063/1.2772250.
  * Chauveteau (1982) Chauveteau, G. 1982 Rodlike polymer solution flow through fine pores: Influence of pore size on rheological behavior. Journal of Rheology 26 (2), 111–142, arXiv: https://doi.org/10.1122/1.549660.
  * Chen et al. (2005) Chen, M., Rossen, W. & Yortsos, Y. C. 2005 The flow and displacement in porous media of fluids with yield stress. Chem. Eng. Sci. 60 (15), 4183 – 4202.
  * Chevalier et al. (2017) Chevalier, T., Dubey, A. K., Atis, S., Rosso, A., Salin, D. & Talon, L. 2017 Avalanches dynamics in reaction fronts in disordered flows. Phys. Rev. E 95, 042210\.
  * Chevalier et al. (2014) Chevalier, T., Rodts, S., Chateau, X., Chevalier, C. & Coussot, P. 2014 Breaking of non-newtonian character in flows through a porous medium. Phys. Rev. E 89, 023002.
  * Chevalier & Talon (2015) Chevalier, T. & Talon, L. 2015 Generalization of Darcy’s law for Bingham fluids in porous media: From flow-field statistics to the flow-rate regimes. Phys. Rev. E 91, 023011\.
  * Christopher & Middleman (1965) Christopher, R. H & Middleman, S. 1965 Power-law flow through a packed tube. Industrial & Engineering Chemistry Fundamentals 4 (4), 422–426.
  * Coussot (2005) Coussot, P. 2005 Rheometry of pastes, suspensions, and granular materials: applications in industry and environment. John Wiley and Sons.
  * Darcy (1856) Darcy, H. 1856 Les fontaines publiques de la ville de Dijon: exposition et application…. Victor Dalmont.
  * Di Federico (1998) Di Federico, V. 1998 Non-newtonian flow in a variable aperture fracture. Transport in Porous Media 30, 75–86.
  * Duvaut & Lions (1976) Duvaut, Georges & Lions, Jacques Louis 1976 Inequalities in mechanics and physics. Springer.
  * Eberhard et al. (2019) Eberhard, U., Seybold, H. J., Floriancic, M., Bertsch, P., Jiménez-Martínez, J., Andrade, J. S. & Holzner, M. 2019 Determination of the effective viscosity of non-newtonian fluids flowing through porous media. Frontiers in Physics 7, 71\.
  * Entov (1967) Entov, V.M. 1967 On some two-dimensional problems of the theory of filtration with a limiting gradient. Prikl. Mat. Mekh. 31, 820–833.
  * Felisa et al. (2018) Felisa, G., Lenci, A., Lauriola, I., Longo, S. & Federico, V. Di 2018 Flow of truncated power-law fluid in fracture channels of variable aperture. Advances in Water Resources 122, 317 – 327.
  * Forchheimer (1901) Forchheimer, P. 1901 Wasserbewegung durch boden. Z. Ver. Deutsch, Ing. 45, 1782–1788.
  * Gelhar & Axness (1983) Gelhar, L.W. & Axness, C.L. 1983 Three-Dimensional Stochastic Analysis of Macrodispersion in Aquifers. Water Resour. Res. 19, 161–180.
  * Glowinski & Le Tallec (1989) Glowinski, R. & Le Tallec, P. 1989 Augmented Lagrangian and operator-splitting methods in nonlinear mechanics, , vol. 9. SIAM.
  * Géraud et al. (2016) Géraud, B., Jones, S. A., Cantat, I., Dollet, B. & Méheust, Y. 2016 The flow of a foam in a two-dimensional porous medium. Water Resources Research 52 (2), 773–790, arXiv: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2015WR017936.
  * Happel & Brenner (1983) Happel, J. & Brenner, H. 1983 Low Reynolds number hydrodynamics. D. Reidel Publishing Co., Hingham, MA.
  * Hewitt et al. (2016) Hewitt, D. R., Daneshi, M., Balmforth, N. J. & Martinez, D. M. 2016 Obstructed and channelized viscoplastic flow in a hele-shaw cell. Journal of Fluid Mechanics 790, 173–204.
  * Hirasaki et al. (1974) Hirasaki, GJ, Pope, GA & others 1974 Analysis of factors influencing mobility and adsorption in the flow of polymer solution through porous media. Society of Petroleum Engineers Journal 14 (04), 337–346.
  * Keller et al. (1967) Keller, J. B., Rubenfeld, L. A & Molyneux, J. E. 1967 Extremum principles for slow viscous flows with applications to suspensions. Journal of Fluid Mechanics 30 (1), 97–125.
  * Kostenko & Talon (2019) Kostenko, R. & Talon, L. 2019 Numerical study of bingham flow in macroscopic two dimensional heterogeneous porous media. Physica A: Statistical Mechanics and its Applications 528, 121501\.
  * Lavrov (2013a) Lavrov, A. 2013a Numerical modeling of steady-state flow of a non-newtonian power-law fluid in a rough-walled fracture. Computers and Geotechnics 50, 101 – 109.
  * Lavrov (2013b) Lavrov, A. 2013b Redirection and channelization of power-law fluid flow in a rough-walled fracture. Chemical Engineering Science 99, 81 – 88.
  * Lavrov (2015) Lavrov, A. 2015 Flow of truncated power-law fluid between parallel walls for hydraulic fracturing applications. Journal of Non-Newtonian Fluid Mechanics 223, 141 – 146.
  * Le Blay et al. (2020) Le Blay, M., Adda-Bedia, M. & Bartolo, D. 2020 Emergence of scale-free smectic rivers and critical depinning in emulsions driven through disorder. Proceedings of the National Academy of Sciences 117 (25), 13914–13920, arXiv: https://www.pnas.org/content/117/25/13914.full.pdf.
  * Liu et al. (2019) Liu, C., De Luca, A., Rosso, A. & Talon, L. 2019 Darcy’s law for yield stress fluids. Phys. Rev. Lett. 122, 245502.
  * Lopez et al. (2003) Lopez, X., Valvatne, P. H. & Blunt, M. J. 2003 Predictive network modeling of single-phase non-newtonian flow in porous media. Journal of Colloid and Interface Science 264 (1), 256 – 265.
  * Matheron (1967) Matheron, G. 1967 Eléments pour une théorie des milieux poreux. Masson Paris.
  * McKinley et al. (1966) McKinley, R. M., Jahns, H. O., Harris, W. W. & Greenkorn, R. A. 1966 Non-newtonian flow in porous media. AIChE Journal 12 (1), 17–20, arXiv: https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.690120106.
  * Mei & Auriault (1991) Mei, C. C. & Auriault, J.-L. 1991 The effect of weak inertia on flow through a porous medium. Journal of Fluid Mechanics 222, 647–663.
  * Mourzenko et al. (1995) Mourzenko, V.V., Thovert, J.-F. & Adler, P.M. 1995 Permeability of a single fracture; validity of the reynolds equation. J. Phys. II France 5 (3), 465–482.
  * Nash & Rees (2016) Nash, Sarah & Rees, D Andrew S 2016 The effect of microstructure on models for the flow of a bingham fluid in porous media. Transp. Porous Media .
  * P. et al. (2006) P., Christian L., Tardy, P. M.J., Sorbie, K. S. & Crawshaw, J. C. 2006 Experimental and modeling study of newtonian and non-newtonian fluid flow in pore network micromodels. Journal of Colloid and Interface Science 295 (2), 542 – 550\.
  * Pearson & Tardy (2002) Pearson, JRA & Tardy, PMJ 2002 Models for flow of non-newtonian and complex fluids through porous media. J. Non-Newtonian Fluid Mech. 102 (2), 447–473.
  * Reynolds (1886) Reynolds, Osborne 1886 On the theory of lubrication and its application to mr. beauchamp tower’s experiments, including an experimental determination of the viscosity of olive oil. Phil. Trans. R. Soc. Lond. 177, 157–234.
  * Rodríguez de Castro & Radilla (2017) Rodríguez de Castro, A. & Radilla, G. 2017 Non-darcian flow of shear-thinning fluids through packed beads: Experiments and predictions using forchheimer’s law and ergun’s equation. Advances in Water Resources 100, 35 – 47.
  * Roquet & Saramito (2003) Roquet, N. & Saramito, P. 2003 An adaptive finite element method for Bingham fluid flows around a cylinder. Computer Methods in Applied Mechanics and Engineering 192 (192), 3317–3341.
  * Rossen (1990) Rossen, William R 1990 Theory of mobilization pressure gradient of flowing foams in porous media: I. incompressible foam. J. Colloid Interface Sci. 136 (1), 1 – 16.
  * Roustaei et al. (2016) Roustaei, A., Chevalier, T., Talon, L. & Frigaard, I. A. 2016 Non-Darcy effects in fracture flows of a yield stress fluid. Journal of Fluid Mechanics 805, 222–261.
  * Sadowski & Bird (1965) Sadowski, T. J. & Bird, R. B. 1965 Non‐newtonian flow through porous media. i. theoretical. Transactions of the Society of Rheology 9 (2), 243–250, arXiv: https://doi.org/10.1122/1.549000.
  * Santucci et al. (2011) Santucci, S., Planet, R., Maloyy, K. Jorgen & Ortin, J. 2011 Avalanches of imbibition fronts: Towards critical pinning. Europhys. Lett. 94 (4), 46005.
  * Shah & Yortsos (1996) Shah, C.B. & Yortsos 1996 The permeability of strongly disordered systems. Phys. Fluids 8, 280–282.
  * Shah & Yortsos (1995) Shah, C. B. & Yortsos, Y. C. 1995 Aspects of flow of power-law fluids in porous media. AIChE J. 41 (5), 1099–1112.
  * Shahsavari & McKinley (2016) Shahsavari, Setareh & McKinley, Gareth H. 2016 Mobility and pore-scale fluid dynamics of rate-dependent yield-stress fluids flowing through fibrous porous media. Journal of Non-Newtonian Fluid Mechanics 235, 76 – 82.
  * Sinha et al. (2017) Sinha, Santanu, Bender, Andrew T., Danczyk, Matthew, Keepseagle, Kayla, Prather, Cody A., Bray, Joshua M., Thrane, Linn W., Seymour, Joseph D., Codd, Sarah L. & Hansen, Alex 2017 Effective rheology of two-phase flow in three-dimensional porous media: Experiment and simulation. Transport in Porous Media 119 (1), 77–94.
  * Slattery (1967) Slattery, John C. 1967 Flow of viscoelastic fluids through porous media. AIChE Journal 13 (6), 1066–1071, arXiv: https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.690130606.
  * Sochi & Blunt (2008) Sochi, Taha & Blunt, Martin J. 2008 Pore-scale network modeling of ellis and herschel-bulkley fluids. J. Pet. Sci. Eng. 60 (2), 105 – 124.
  * Sorbie et al. (1989) Sorbie, K.S, Clifford, P.J & Jones, E.R.W 1989 The rheology of pseudoplastic fluids in porous media using network modeling. Journal of Colloid and Interface Science 130 (2), 508 – 534.
  * Sorbie (2013) Sorbie, Kenneth S 2013 Polymer-improved oil recovery. Springer Science & Business Media.
  * Stauffer & Aharony (1991) Stauffer, Dietrich & Aharony, Amnon 1991 Introduction to percolation theory. Taylor and Francis.
  * Straley (1977) Straley, Joseph P 1977 Critical exponents for the conductivity of random resistor lattices. Physical Review B 15 (12), 5733.
  * Tallakstad et al. (2009) Tallakstad, Ken Tore, Knudsen, Henning Arendt, Ramstad, Thomas, Løvoll, Grunde, Måløy, Knut Jørgen, Toussaint, Renaud & Flekkøy, Eirik Grude 2009 Steady-state two-phase flow in porous media: Statistics and transport properties. Phys. Rev. Lett. 102 (7), 074502–.
  * Talon et al. (2014) Talon, L., Auradou, H. & Hansen, A. 2014 Effective rheology of Bingham fluids in a rough channel. Front. Physics 2 (24), 24.
  * Talon & Bauer (2013) Talon, L. & Bauer, D. 2013 On the determination of a generalized Darcy equation for yield-stress fluid in porous media using a lattice-Boltzmann trt scheme. Eur. Phys. J. E 36 (12), 139\.
  * Talon & Hansen (2020) Talon, Laurent & Hansen, Alex 2020 Effective rheology of bi-viscous non-newtonian fluids in porous media. Front Phys 7, 225.
  * Tongwa et al. (2013) Tongwa, Paul, Nygaard, Runar, Blue, Aaron & Bai, Baojun 2013 Evaluation of potential fracture-sealing materials for remediating co2 leakage pathways during co2 sequestration. International Journal of Greenhouse Gas Control 18, 128–138.
  * Yiotis et al. (2019) Yiotis, A. G., Dollari, A., Kainourgiakis, M. E., Salin, D. & Talon, L. 2019 Nonlinear darcy flow dynamics during ganglia stranding and mobilization in heterogeneous porous domains. Phys. Rev. Fluids 4, 114302.
  * Yiotis et al. (2013) Yiotis, A. G., Talon, L. & Salin, D. 2013 Blob population dynamics during immiscible two-phase flows in reconstructed porous media. Phys. Rev. E 87, 033001.
  * Zami-Pierre et al. (2016) Zami-Pierre, F., de Loubens, R., Quintard, M. & Davit, Y. 2016 Transition in the flow of power-law fluids through isotropic porous media. Phys. Rev. Lett. 117, 074502.
  * Zimmerman et al. (1991) Zimmerman, R.W., Kumar, S. & Bodvarsson, G.S. 1991 Lubrication theory analysis of the permeability of rough-walled fractures. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28 (4), 325–331.