{'fraction_non_alphanumeric': 0.038243067334621905, 'fraction_numerical': 0.008954281882371504, 'mean_word_length': 4.3807704552690225, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the statistical properties of fluid flows with transitional power-law rheology in heterogeneous porous media

In this work, we study non-Newtonian fluid flow in heterogeneous porous media. We are interested in fluids presenting a specific change in rheology: Newtonian below a certain shear rate and power law above. Since porous media generally exhibit strong spatial heterogeneity at large geological scales, we study the interaction between such inhomogeneity and the nonlinear rheology of the fluid. The coupling between permeability heterogeneity and nonlinear rheology significantly affects the flow. We are particularly in the statistical properties of the velocity field (mean, variance, correlation, etc). Depending on the imposed mean pressure gradient, three macroscopic flow regimes are identified. For a low or high average pressure gradient, the average flow rate increases linearly or according to a power law, respectively. In the latter regime, we observe that the velocity field is more heterogeneous for shear-thinning fluids than for shear-thickening fluids. This is corresponding to a channeling effect of shear-thinning fluids. The intermediate regime corresponds to a progressive and inhomogeneous change of the local rheology. This transient regime is then characterized in terms of pressure gradient range. The flow field is also analyzed statistically. The spatial distribution of the regions above the rheology threshold shows interesting statistical properties. For instance, they exhibit multiscale characteristics (fractal), similar to other critical systems (percolation, avalanches, etc.). More surprisingly, even though some statistical properties are independent of the parameters, an interesting abrupt rotation of the correlations is found for a particular set of parameters. This is explained by using some symmetries of the problem.


Introduction
Many natural or industrial fluids exhibit non-Newtonian behaviors, [1][2][3], they are thus found in many applications related to porous or fractured media. A very important application is certainly enhanced oil recovery (EOR) (see [4]). The most standard way to recover oil is to inject water to force the oil to move (waterflooding). The efficiency of waterflooding is however related to the uniformity of the displacement front. The more uniform the front, the greater the amount of oil to be moved. In contrast, if the preferential paths are formed, they will bypass certain regions of the porous medium leaving the oil in place. Two main reasons can cause this bypass. The first one comes from the heterogeneity of the permeability field. The fluid indeed tends to flow into the highest permeability regions and avoid the lower permeability ones. The second reason comes from the viscosity contrast. If the displacing fluid has a lower viscosity than the displaced one, the front can destabilize and form fingerings as described by Saffman and Taylor [5] or Homsy [6]. To limit this bypass effect, a solution used by the oil industry is to add a polymer (e.g. Xanthan) to the displacing water (see [4]). Another advantage of improving the viscosity ratio lies in the fact that it also reduces the effects of heterogeneity in the permeability field. The main difficulty of this process is that most of the injected polymers have a non-Newtonian rheology. This is the case for example of Xanthan which is Newtonian at low shear rates but is shear thinning at higher shear rates. However, as has been shown in pore networks [7] or in fractures [8], nonlinear rheology can amplify or dampen heterogeneities. The objective of this paper is to further investigate the interplay between the macroscopic inhomogeneity of a permeability field and the flow of a nonlinear fluid. We focus particularly on the change in behavior of the viscosity.
Another interesting application of non-Newtonian fluids is the understanding of blood flow in capillary networks. Indeed, blood, which must be considered as a suspension, presents a non-Newtonian viscosity: shear thinning or yield stress. [8,9]. Non-Newtonian fluid is also present in the proppants used for hydraulic fracturing or the mud produced by drilling wells [10,11]. They are also commonly used for fracture sealing (e.g cements, polymers, etc.) [12].
A very recurrent problem when dealing with porous media is that of upscaling. If the equations of motion are generally well known at the pore scale (typically ∼ 10 −3 m), a particular interest is to understand the flow at much larger scales (∼ 1 − 10 3 m) This is usually done by deriving constitutive equations for average quantities at an intermediate scale and is illustrated by the famous Darcy's law for Newtonian fluids, which relates linearly the mean flow rate to the macroscopic gradient of pressure.
At the microscopic level, Newtonian fluids obey the Stokes equation (neglecting the inertia) : where v is the fluid velocity, p is the pressure and µ the viscosity. Averaging the velocity and pressure field over a large number of pores results in Darcy's law [13]: where u is the mean velocity, ∇P an averaged pressure gradient and κ the permeability of the porous medium which depends on its structure. At the geological scale, the type of rock may however spatially vary leading to a macroscopic heterogeneous permeability field. The large-scale flow obeys then to the heterogeneous Darcy's law: It is also worth recalling that a very similar equation is used for solving flow in rough fractures, usually referred as the Reynolds equation [14][15][16]. Indeed, in a fracture with varying opening and under the lubrication approximation (small thickness and small variation of the opening), the flow obeys: where the fracture is in the (x, y)-plane, b(x, y) represents the local opening, κ(x, y) = b 3 (x,y) 12 and q(x, y) = b(x,y) 0 u dz is the flow rate. 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 the permeability has the dimension of m 2 , whereas in a fracture, q is a volumetric flux per unit length (m 2 .s −1 ) and the "permeability" has the dimension of m 3 .
The influence of the heterogeneity of a permeability field (or fracture) has also been the subject of a considerable amount of work, starting with the work of Matheron [17], Gelhar and Axness [18] or Dagan [19]. One can also mention the review by Renard and de Marsily [20]. Since the fluid prefers to flow in high permeability regions and avoid low permeability ones, the heterogeneity strongly influences the velocity field. The central question is to understand how the permeability distribution affects the velocity field, in particular its average rate and the magnitude of its heterogeneity. The latter is specifically important for describing the transport of species (e.g., pollutants) in the subsurface [21].
All studies mentioned above apply to Newtonian fluids. Therefore, the following question naturally arises: how should these approaches be modified when considering non-Newtonian fluids? Although there exists a very large variety of non-Newtonian fluids [1,3,22], several similar approaches are commonly used in the case where it exists a simple relationship between shear rateγ = 1 2 tr [( ∇ v + ∇ v T ) 2 ] and viscosity 2 τ γ τ 0 γ 0 ∇P u g c u c Figure 1: Left: Relationship between shear rate and shear stress for a truncated rheology. For a shear rate lower than γ 0 (or a stress lower than τ 0 ), the fluid is Newtonian. Beyond this thresholds the rheology is a power law. Right: General Darcy's law, mean flow rate u as function of the average gradient of pressure ∇P in porous media for the truncated rheology using the "mean field" rheology approach. µ(γ), or equivalently betweenγ and shear stress τ (γ), with τ (γ) = µ(γ)γ. One approach (see for instance [23][24][25][26][27]) consists in determining an effective shear rateγ pm to derive an effective viscosity. Other ones are to establish an effective stress [30], or an average viscosity [31]. The common feature of these approaches is that they are based on the determination of mean effective quantities. They can be synthesized using scaling arguments. Indeed, by defining a typical length scaleλ (pore size, grain diameter, √ κ, etc.) and using the average flow rate u, a typical shear rateγ ef f ∝ u/λ can be defined. A typical shear stress τ ef f ∝λ∇P is also deduced from the pressure gradient. These quantities can then be used in the rheological functioṅ γ = f (τ ) to derive a generalization of Darcy's law in the form : 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 (see Fig. 2). Moreover, Shah and Yortsos [7] and Auriault [32] proposed a theoretical approach to homogenize the flow for power law fluids in a periodic porous media. In both cases, the derived Darcy's law is also a power-law which is in agreement with the effective quantities approach. The objective of this paper is to study rheologies that exhibit a change in behavior. Indeed, many non-Newtonian fluids, such as Xanthan, exhibit non-linear behavior only at high shear rate (or stress). At low shear rates, they still behave like a Newtonian fluid. To describe this rheology, a simple model is chosen which is the "truncated rheology", where the transition is sharp, at a given value of shear rate (or shear stress): where n is the flow index. Using the "mean field" approach, as studied numerically by [33,34] for example, Darcy's law can be written as follows: where κ is the permeability and u c a velocity threshold that depends onγ 0 and the porous structure. The prefactor in the second equation is determined by continuity. κ and u c are thus parameters depending on the topology of the porous medium (pore size distribution, porosity, etc.). At macroscopic geological scales, however, the structure of porous media is expected to vary, resulting in a heterogeneous permeability and u c ( r) field. The main objective of this paper is therefore to study the influence of this field inhomogeneity on the flow for a fluid presenting a transient rheology. Another application could be also to evaluate the use of non-Newtonian fluids to characterize the degree of heterogeneity of a field. Indeed, as each location changes its viscosity behavior at a different flow rate, recording the evolution of the average velocity could potentially give indications on the permeability heterogeneity. This article is structured as follows. Section 2 is devoted to the presentation of the problem. Section 3 contains the numerical results, where the different flow regimes are analyzed. In particular, the transient regime presents interesting statistical properties. Section 4 is dedicated to discussions and conclusions. Different appendices provide some mathematical properties of the nonlinear Darcy's law (perturbation expansion, symmetry) and also the numerical method used.

Problem description -Governing equations
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: In addition, the velocity field must satisfy mass conservation: The flow field is solved using a second order finite difference method combined with an augmented Lagrangian approach described in the appendix. In practice, the flow rate is determined by imposing a pressure difference ∆P between the inlet and the outlet (see Fig. 2). However, it is more convenient to use the average pressure gradient ∇P = ∆P L as a control parameter, whereL is the length of the system and the mean operator is defined as . = 1 WL . dxdy. The lateral boundary conditions are assumed to be periodic. By construction, ∇P has therefore a constant direction, aligned with the x-axis.
The permeability field was chosen to be distributed according to a log-normal distribution which has been a common model since the work of Gelhar and Axness [18]. It has the advantage of being consistent with the field data but also of allowing the mean and variance to be varied independently. The permeability field is obtained by generating a Gaussian field δf of zero mean and given standard deviation σ f . δf has a correlation length λ (see [35] for more details). The permeability field is then given by: where κ 0 = exp (f 0 ), is a parameter characterizing the average permeability of the medium. The threshold field u c is expected to depend on bothγ 0 and the pore structure.γ 0 is a characteristic of the fluid rheology and is therefore constant. The pore structure may however vary spatially and is related to the permeability κ. A simple relationship can be established using phenomenological arguments. For porous media with a typical pore size d, a scaling analysis leads to κ ∼ d 2 et u c ∼ γ 0 d, which gives: This scaling has for example been evaluated in [34]. It is important to stress that this argument is very crude. While this scaling law is certainly valid in the case of homothetic transformations, it is not necessarily applicable to more complex structural changes. In other words, this scaling law is most likely valid when changing the diameter of a packet of monodisperse beads (or sand). But it is probably more complicated if the nature of the medium changes radically, from a sandpile to a porous rock for example.
As mentioned in the introduction, two-dimensional Darcy's law can also be used to solve the flow in heterogeneous fractures. The scaling is then slightly different. By defining h as the opening, the local permeability leads to κ ∼ h 3 and u c ∼ γ 0 h 2 . It follows: The field u c is thus determined from the permeability field κ(x, y) by assuming a more generic scaling law of the form: where A is a prefactor and γ a parameter with γ ∈ [0; 1]. Although there are many different parameters, in this study we will focus on the disorder amplitude σ f , the flow index n and the exponent γ. As discussed in Appendix C, the numerical method is particularly efficient for specific values of n ∈ [1/3; 1/2; 2/3; 1; 3/2; 2; 3], which are the values used in this work. The parameters A, µ and κ 0 will be kept constant with A = 1, µ = 0.1 and κ 0 = 1.
Based on these parameters, a characteristic velocity u 0 = Aκ γ 0 and a characteristic pressure gradient g 0 = µu0 κ0 are defined and will be used for non-dimensionalization. Basically, u 0 and g 0 represent the velocity and pressure gradient at which the system would change behavior if the field were homogeneous (i.e σ f = 0).

Simple case: Heterogeneous permeability -homogeneous critical velocities
In order to describe the problem qualitatively, a simplified version is presented here, where the permeability is heterogeneous but the critical velocity field u c is homogeneous (i.e γ = 0). Figure 3 (left) shows the flow field for different average pressure gradients ∇P . For a very small pressure gradient, the flow field is heterogeneous but all velocities are below the (single-valued) threshold u c = u 0 (figure not shown). As the pressure gradient increases, the velocity field increases as a whole. At some point, some regions reach the threshold (Fig. 3, top). At these locations, the viscosity is changed, which disturbs the surrounding flow field. This modification therefore favor or disfavor the velocity in the close vicinity.  : For a shear-thinning truncated rheology n = 1/2 in a medium with a heterogeneous permeability field with σ = 1 with a homogeneous threshold field (γ = 0). As function of the mean pressure gradient: a: mean velocity ux . The red dashed lines represent the power law ux ∝ ∇P and ux ∝ ∇P 1/n . b: mean velocity rescaled with the Newtonian Darcy's law. c: percentage of regions above the threshold (|| u|| > uc). d: relative standart deviation of the velocity field σu ux . Vertical dashed lines represent the transitions of the flow regimes (I, II and III).
Here, the fluid is shear-thinning n < 1, so the viscosity is locally decreased. The flow is then increased downstream and upstream, while it is decreased on the lateral sides. Thus, some correlations of regions that have changed flow regime can be expected. For shear-thinning (resp. thickening), the regime change should thus be correlated along the flow direction (resp. perpendicular to it). As the pressure gradient increases, more and more regions change their behavior until the entire domain is in the non-Newtonian regime. We can also note how the velocity distribution is crucially altered as it passes through the threshold value in Fig. 3 (right column). Fig. 4 shows the evolution of various interesting quantities. Fig. 4.a displays the average velocity as a function of the average pressure gradient. As expected, for low pressure, the average flow rate varies linearly with the average pressure gradient. Above a certain value, the curve starts to slowly deviate from the linear trend to a power law type trend. Three flow regimes are thus identified. For the first one (regime I), the average flow rate increases linearly. At a high flow rate (regime III), the mean flow follows a power law with exponent u x ∝ ∇P 1/n . Between these two regimes, a transitional one (regime II) is observed. This change in behavior can be magnified by plotting the ratio u x / ∇P for example ( Fig.4.b). Another convenient quantity to characterize this transition is the percentage of regions in the nonlinear regime R( ∇P ) ∈ [0, 1] ( Fig. 4.c). A final very important quantity, shown in Fig. 4.d, is the relative standard deviation of the velocity field, σu ux , which characterizes the heterogeneities of the flow and is of great importance, for example, in the problem of species transport. Here we observe that regime I and III correspond to a plateau value of the standard deviation.
This example shows how permeability heterogeneity affect the flow regime. Since regions with high permeability have higher local velocity, they are more likely to reach the nonlinear viscosity regime. However, the more general case γ = 0 is more complex because the velocity threshold is also distributed in space. Additionally, regions of higher permeability also correspond to a higher velocity threshold. A competition between these two effects is thus expected. .a displays the average flow rate as a function of the average pressure gradient ∇P for a shear thinning fluid, n = 1/2, for γ = 1/2 and for different magnitudes of heterogeneity σ f . Similarly to the homogeneous velocity threshold case, a transient behavior is observed around u x ∼ u 0 and ∇P ∼ g 0 . In Regime I, the flow rate follows a linear behavior u x ∝ ∇P , while in regime III, it follows a power law u x ∝ ∇P 1/n . The disorder then smoothes the transition while increasing its range. For σ f = 0, the transition is abrupt, while for σ f high, the transition is smoother and extending over half a decade. This effect is more apparent after normalizing the mean flow by Darcy's Newtonian law µ ux κ0∇P ( Fig. 5.b). The second effect of heterogeneity is to significantly increase the flow rate in the high pressure gradient regime, where the flow rate is almost doubled between σ f = 0 and σ f = 2. This effect will be further detailed later as it can be predicted using a perturbation expansion approach (Appendix A). Fig. 5.c plots the relative standard deviation of the flow field, σux ux , as a function of the mean pressure gradient and for different σ f . Like previously, we observe a change of plateau when changing the flow regime. And the values of these plateaus increase with the magnitude of the permeability heterogeneity σ f . Here again, the range of the transition between these two asymptotes depends on the heterogeneities of the porous medium, abrupt at low σ f while smoother at high σ f . We can also notice that these curves are very symmetrical. In fact, these curves can be fitted remarkably well with a hyperbolic tangent function (in semi-log representation). Fig. 5.d shows the normalized standard deviation of the velocity as a function of the applied pressure gradient, for a shear-thickening fluid (n = 2). The most notable difference is that the heterogeneity of the velocity field decreases in regime III compared to regime I. Shear-thickening fluids therefore attenuate the permeability field heterogeneity while shear-thickening fluids enhance it.

Influence of the amplitude of the field heterogeneities
We have seen that the normalized standard deviation is constant in regimes I and III. It is noteworthy to mention that this is also the case for all normalized moments of the velocity distribution (Skewness, Kurtosis, etc.). In fact, the velocity distribution actually keeps a constant shape in both regimes as observed in Fig.  6. The distribution is only shifted when the mean pressure gradient varies. As a result, the field u/ u x is

Statitical properties of regime I and III
The flow field in regimes I and III is governed by the equation of a power-law fluid. For such fluid, it is however possible to determine the mean flux and standard deviation of the velocity using a perturbation expansion approach. The principle is to extend the work of [18] on Newtonian fluids to power law fluids. It consists in expanding the permeability, pressure and velocity field around the mean value and in assuming that the deviation terms are small for sufficiently small σ f . The complete calculation is provided in Appendix A for a governing equation of the form : − ∇P = c( r)|| u|| n−1 u, where we define r = (x, y) for the sake of conciseness. In the case of the truncated rheology, we have c( r) = µ κ( r)u n−1 c . Because κ and u c are distributed according to a log-normal distribution, it is also the case for c( r) with: σ log c = (1 + γ(n − 1))σ f . It follows: and with α = 1/n and D 0 = µ −α κ α+γ(1−α) 0 A 1−α . These results confirm qualitatively the previous observations that lower n increase the mean flow rate and the flow heterogeneity. As observed in pore network model [7], shear-thinning fluids are indeed more channelized and thus more heterogeneous. The limit n → 0 is interesting because it corresponds to a pure plastic flow (viscoplastic fluid without viscosity). The flow field is then expected to converge to a single flow path as for the Bingham fluid in the large Bingham number limit [35]. However, the expansion of the perturbations predicts a divergence of the standart deviation of the velocity at n = 0. Since the calculations assume small perturbations for each field, they are expected to fail for a sufficiently small value of n. On the other hand, increasing n decreases the average velocity and the heterogeneity of the flow. The expansion approach should thus be better. Fig. 7 compares the numerical mean and standard deviation of the flow field for a sufficiently small (resp. large) applied pressure gradient ∇P /g 0 = 0.01 (resp. ∇P /g 0 = 100) against the predictions of eqs. (12) and (13). Both figures show a very good agreement between the analytical predictions and the simulations, even for a heterogeneity amplitude as high as σ f ∼ 2, which is quite significant for a lognormal distribution. An expansion of order 2 seems thus sufficient to predict both the mean flux and standard deviation of the flow field in regimes I and III.

Pressure range of regime II
We now discuss the range of the transient regime II. Fig. 8 represents different phase diagrams of the system as a function of the mean pressure gradient while varying σ f , γ or n. The crosses represent the system in the intermediate regime whereas the dots represent the system in regime I or III. Fig. 8.a displays the evolution of the transient regime as a function of σ f . This figure mainly confirms the previous observation that the range increases significantly with σ f because the velocity field is more heterogeneous. Considering that the plot is on a logarithmic scale, the growth is in fact very significant as the range increases faster than a power-law. Fig. 8.b shows the variation of the pressure range of regime II as a function of the rheological index n. The transition range becomes narrower with increasing n. This result is contradictory to the previous observation that shear thickening fluids have lower velocity heterogeneity than shear thinning ones. It can also be noted that the lower pressure limit is not really affected by the value of n.
The influence of γ is also not intuitive (Fig. 8.c) because the pressure range exhibits a non-monotonic behavior. The pressure range initially decreases with γ but increases again above a certain value of γ.
The reason behind the last two observations is that the transition depends not only on the velocity distribution but also on the distribution of u c ( r). The connection between both distributions is not obvious because regions of higher velocity are more likely to have higher permeability, and thus a higher threshold. One way to highlight this phenomenon is to estimate the limits of regime II by exploiting the invariance of the velocity field in regimes I and III.
Starting at a very low flow rate, all regions are in the Newtonian regime and thus || u( r)|| = u x || ω 1 ( r)||. When increasing the average flow rate, the first location that changes its rheology occurs where || u( r)|| = u c ( r). It follows the condition of this first occurrence: max ux || ω1( r)|| uc( r) = 1. This means that the change of regime will start at the average velocity u x = u 1 satisfying: Similarly, starting with a high pressure gradient, all regions are in the power-law regime and u( r)/ u x = ω n ( r) is constant. While decreasing u x , the first point changing its behavior occurs when u x = u 2 , with: Assuming that at order zero ∇P 1 µ κ0 u 1 and ∇P 2 µκ 0 −1−γ(n−1) A 1−n u n 2 , it leads to an estimation for the pressure bounds of the transient regime: Lower bound: Upper bound: These two boundary estimates have been plotted in Fig 8 and show good correspondence with the numerical simulations. This shows also that the relevant quantity is in fact the extension of || ωm( r)|| uc( r) (with m = 1 or n) which does not necessarily follow the extension of the velocity field ω m . This is confirmed in Fig. 9 which represents the evolution of the standart deviation of || ω n ( r)|| (the normalized velocity in regime III) and the standart deviation of || ωn( r)|| uc( r) as function of γ. For low γ values, both quantities decrease with γ. However || ωn( r)|| uc( r) is non-monotonic and becomes increasing above a certain value of γ. This fact explains then the non-monotonic evolution observed in Fig. 8.b.

Statistical properties of the flow field in the transient regime II
The flow field in the transient regime has interesting statistical properties. One way to apprehend it, is to notice a similarity with the problem of percolation. Indeed, as shown in Fig. 3, when u x is increased, more and more regions satisfy the criterion ||u|| > u c . Connected regions satisfying this criterion allows to define clusters, which become larger and more numerous as the pressure gradient increases. If the field || u||/ u x = || ω( r)|| were constant, the transition would occur at which would correspond to a percolation problem. Similar behaviors are therefore expected, such as cluster fractality and the presence of criticality. However, it is important to recall that the problem is not strictly equivalent to percolation because the field u ux is not constant in the transient regime. In particular, the change in viscosity introduces correlations in the velocity field and thus changes the shape of the clusters.
Here, the clusters are identified using a Hoshen-Kopelman algorithm [36]. Their shape are then characterized by their total size S and the two dimensions of the bounding rectangle that contains it: L along the flow direction and W transversely to it.  Fig 10.a displays the size distribution P (S) for different mean pressure gradient and for the parameters n = 1/2 and σ = 1. For any applied pressure, the distribution follows a decaying power law over a wide range of sizes. However, a large-scale cutoff, S 0 , is found that varies with the pressure gradient. The plotting of the variation of S 0 as a function of pressure ( Fig. 10.b) shows that S 0 diverges at a certain value of the average pressure gradient ∇P = ∇P c according to a power law:

Size distribution
Combining both observations, it follows the scaling law for the size distribution: which is confirmed by the good overlap of the distributions rescaled according to eqs. (18) and (19) plotted in Fig. 10.c. This scaling law is thus similar to the one found in other problems with a critical transition such as percolation [37], avalanches (e.g [38,39]) or yield-stress fluid in porous media [35,40]. The sizes are distributed on many scales up to a size limit. And this size limit diverges as the control parameter, ∇P , approaches a critical value ∇P = ∇P c . The system then exhibits an infinitely broad range of scales (e.g. fractal). This scaling law is characterized by the two exponents τ S and ν S . 14 A similar scaling law could also be observed for the length L of the clusters: allowing to also identify the exponents τ L and ν L . For each set of parameters (σ f , γ, n) this scaling law is observed leading to the determination of the exponents τ S , τ L , ν S and ν L . Tables 1-3 report these exponents according to the different sets of parameters. The most remarkable result is the fact that τ S seems to be indeed independent of the parameters. This is a characteristic found in many critical systems, where some exponents are independent of the details of the disorder distribution. Such behavior is often referred to as "universal" as, for example, the exponents τ and ν in the percolation or avalanches of an elastic line in a random medium [41]. The results obtained here seem to suggest a universal behavior for the exponent τ S . The observed value τ S = 1.65 ± 0.05 is however very different from the standard percolation problem (τ perc = 2.05), which indicates that it would be of a different universality class.
The trend is less clear with the exponent ν S , which seems to vary with the rheological index n and also the heterogeneities σ f . However, it should be noted that the determination of the exponent ν S is generally more prone to errors because it requires the determination of ∇P c , which is also subject to uncertainties. The error can be estimated at about 10%. It is then difficult to conclude on the universality of this exponent.

Cluster's shape
The shape of the clusters can be characterized by the aspect ratio L/W . Fig. 11.a shows the width W of each cluster as a function of their length L for two different sets of parameters. A power-law type relationship is then observed: which is a characteristic of the self-affine fractal structure typically found in anisotropic critical systems (e.g. avalanches, directed percolation, front propagations, etc.). It characterizes the fact that, although many different cluster sizes are present, the aspect ratio is not the same at each scale. If ζ < 1, larger clusters are more elongated than the smaller ones as represented in Fig. 12.a. For ζ > 1 bigger clusters are more elongated in the direction transverse to the flow ( Fig. 12.b). Because clusters are not compact, viz. they may contain holes, another interesting quantity to analyze is the surface area (S) as a function of the enclosing box size W L, as shown in Fig. 11.b. Here, again the relationship observed is a power-law: which is also a characteristic of a fractal structure. The measured exponents for different sets of parameters, σ, n and γ, are displayed in table 1-3. The exponent of size β seems to be almost constant, within the error bar, for any parameter value: β 0.77±0.05. This exponent seems to be universal. More surprising is the evolution of ζ which takes only two values: either ζ 0.85 ± 0.05 or ζ 1.15 ± 0.05 depending on the parameters. The two cases shown in Fig. 12 are in fact the only two observable exponents. In this figure, the two cases appear to be very similar, as if they were rotated by 90 • . In fact, such rotational symmetry can be proven when the rheological parameter n is modified. As detailed in the Appendix B, rotation by 90 • of the velocity field and pressure field gradient is equivalent to solving the flow with inverse rheology (i.e n → 1/n).
We can therefore expect that, by reversing n → 1/n, the shape of the clusters remains the same but rotated by 90 • . This corresponds to the inversion of the roles of W and L. It thus leads to the relation: which seems to be satisfied by the numerical observations. The value ζ(n = 1) 1 appears therefore as a marginal value. It remains quite remarkable that this exponent is constant for any n exponent of shear thinning (or shear thickening).   More unexpected is the similar change when γ is varied ( Table 3). The value ζ is constant until a certain value γ ∼ 0.6, where it switches to its inverse value. Some remarks can be made to interpret this switch. First, this change in correlation confirms the fact that the spatial correlation of the velocity field, correlated in the direction of flow, does not necessarily follow the correlation of the field u( r)/u c ( r). Moreover, the value of γ at which the change occurs, seems to correspond to the non-monotonic change observed in Fig.  9.a. To understand this change, it is useful to consider the two extreme cases γ = 0 and γ = 1.
Before, it should be noted that in the non-linear Darcy's law eq. (8), the criterion || u|| > u c ( r) also corresponds to a criterion for local pressure gradient || ∇P ( r)|| > g c ( r), with u c = Aκ γ and g c = µAκ γ−1 .
For γ = 0, u c is a constant: the cluster, defined by || u|| > u c , are expected to follow the correlation of the velocity field. The cluster should then be more elongated in the streamwise direction. For γ = 1, the situation is different because u c ( r) is now a random variable, but g c is a constant field. The clusters, equivalently corresponding to || ∇P ( r)|| > g c , should then have a similar direction of correlation as ∇P . In general, the pressure gradient tends to be correlated in the direction transverse to the flow 1 . The clusters are then expected to have a correlation direction transverse to the flow. Depending on the γ values, the shape of the cluster is then the result of a balance between these two opposite effects.
If the cluster elongation changes with γ, it is still quite surprising that the correlation exponent ζ is constant and only switches from one value to another. This suggests that the exponent is also universal but the principal direction of correlation is determined by the value of γ.  Table 3: Measured exponents for different γ with (n, σ f ) = (0.5, 1).

Conclusion
In this paper, we have studied the flow in macroscopic heterogeneous porous media with nonlinear rheology exhibiting a change in behavior such as truncated rheology. The influence of heterogeneities in permeability and velocity threshold fields was considered.
By varying the mean pressure gradient, three flow regimes are observed. At low pressure gradient, the whole system is in the constant viscosity regime. The total flow rate then increases linearly with the mean pressure gradient. At a high pressure gradient, when the entire system is in the nonlinear viscosity regime, the total flow rate increases non-linearly with the same exponent. Transiently, the different regions of the medium change their viscous regime inhomogeneously with the mean pressure gradient, which induces a progressive change of the mean flow.
The first and last regime are relatively easier to analyse since they obey to the equation of a power-law fluid. This allow for instance a good prediction using a perturbation expansion approach.
The intermediate regime is more complex. Qualitatively, the pressure range of this regime is expected to be related to the width of the velocity distribution in regime I and III and thus to the amplitude of the permeability heterogeneities. For a wider velocity distribution, the transition starts at a lower average velocity and ends at a higher velocity. Such behavior could be used to quantify the heterogeneities of a field for example. While this trend was indeed observed, the results showed that it is more complicated to make quantitative predictions. Indeed, a very important aspect of the problem is the relationship between velocity, permeability and local threshold. In this paper, a power-law dependency between the permeability and the threshold has been assumed. The exponent γ modifies drastically the correlation of the clusters.
The reason behind this observation is that different mechanisms are at work. First, regions with high permeability have a higher velocity field but also a higher velocity threshold. If the velocity threshold varies weakly with permeability, the correlation follows that of the velocity field, in the direction of flow. Conversely, if the velocity threshold varies strongly with permeability (i.e. high γ), the transition may not correspond to the highest velocities regions. But at the same time, a similar reasoning can be made by considering the pressure field gradient: low permeability regions are more likely to have a higher pressure gradient but also a higher pressure threshold. The correlations could therefore be dominated by the pressure gradient field, transverse to the flow direction. The competition between these two effects is thus balanced by the γ exponent.
The relationship between the permeability and the two thresholds seems then to be very important. Here, a power law has been assumed but it is important to point out that in practice, the relationship is probably more complicated.
A remark can be made. In the literature, a very common approach is to model a porous medium or a fracture by a bundle of tubes or parallel layers [42][43][44][45]. This facilitates the analysis of the influence of heterogeneity and leads to a good qualitative understanding of the problem. In these models, the flow field is by construction infinitely correlated in the streamwise direction, whereas the pressure gradient is uncorrelated in the crossflow direction. These models are thus expected to be unable to capture behaviors associated with the change in correlation. These models should therefore be taken with caution when applied to a 2D or 3D medium.
The statistical properties of the flow field in the transient regime also appeared very rich. Indeed, in this regime, regions above their threshold define clusters that exhibit fractal and critical properties. 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 universal behaviors. Another interesting feature is that the shape exponent ζ is constant but only switches to its inverse depending on the value of γ and n.
This statistical feature is thus very reminiscent to related problems such as percolation [37] and yield stress fluid in porous media [35]. 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. For yield stress fluid, it was found τ S 1.15 and ζ 0.75. The size distribution exponent τ S is different while the roughness exponent is similar to the present case. As discussed previously, a slightly closer problem could be the percolation (directed or not), but the correlation of the velocity field evolves with the applied gradient of pressure. The observed exponent are indeed different: τ S 2.1 for percolation and τ S 1.26 for directed percolation. The present case seems to be intermediate and if the exponents are "universal" they fall in a different universality class.
There are many interesting directions to pursue this work. One important question is how to generalize it to the 3D permeability field. Although rotational symmetry is no longer applicable in 3D, there remains the important fact that the velocity field is correlated along the stream direction while the pressure gradient is correlated in both transverse directions. Thus, a change in correlation should still be expected depending on γ. Critical behaviors should also probably be observed but with different exponents. Another interesting study would be to investigate the problem of species transport in such a system. Indeed, the dispersion of a tracer depends on the heterogeneity of the velocity field and its correlation. It is therefore expected to observe a change in behavior due to the change in rheology. Furthermore, since the molecular diffusion coefficient is generally related to the viscosity of the fluid, it is expected to be different when the viscosity is below or above the threshold. Finally, another direction of investigation could be other porous media problems with similar behaviors. For example, the problems of two-phase flow [46][47][48][49], emulsion driven in a porous medium [50] or erosion of a granular bed [51], present a similar critical behavior with the appearance of preferential flow paths as a function of flow rate. These problems have in common the property that at a certain critical velocity, the local flow conditions are drastically modified, because the bubbles are mobilized or because the grains rearrange themselves. It would therefore be very interesting to study the similarities and differences between these problems.

Acknowledgements. -
This work is supported by "Investissement d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM). This work was partly supported by the Research Council, through its INTPART funding scheme, project number 309139. I would like to thanks D. Salin, R. Kostenko, A. Hansen and A. Rosso for fruitfull discussions.

Appendix A. Perturbation expansion for a power-law rheology
In section 3, Fig. 5, we have seen that in regimes I and III, the moments of the velocity distribution (mean, standard deviation, etc.) are constant. In these regimes, the flow is governed by Darcy's law for a power law fluid in a heterogeneous medium. In this case, the mean and standard deviation can be determined using a perturbative approach. Following the work of Gelhar and Axness [18] for Newtonian fluids, the principle is to expand the solution around the mean value and assuming sufficiently small perturbations (i.e σ small).
The flow field is assumed to be solution of a power-law rheology in heterogenous porous media in the form: − ∇P = c( r)|| u|| n−1 u, with the free divergence: The field c( r) is assumed to be distributed according to a log-normal. The principle is to decompose each field u, P , and c( r) into a mean part and the spatially fluctuating part To simplify the notations, we introduce G = ∇P and U = u x . The mean flow is assumed to be along the e x axis. It is more convenient to use the g field defined by: c( r) = exp (ḡ + δg) = C 0 exp(δg).
Expanding up to the second order, one have: Expanding and taking the spatial average . of eq. (A.1), all the first order terms vanish by definition. It yields: Along the e x axis, it follows: This expression relates the mean gradient G = ∇P to the mean velocity of the flow U = u x , if we know the different cross-correlation terms of the spatially fluctuating fields, which are determined next.
The first order of eq. (A.1) gives: Thus, Taking the curl leads to: It is now more convenient to write this equation in Fourier space. Definingˆ u(k x , k y ) andĝ(k x , k y ), respectively the Fourier transform of δ u and δg, gives: ik xûy − ik y (Uĝ + nû x ) = 0. (A.10) Using the free divergence in Fourier space: it follows the relationship between the fluctuation of δ u and δg, in Fourier space: From these expressions, it is possible to determine the different cross correlation terms using Parseval's formula: F H dxdy = FĤ * dk x dk y , for any field F ( r) and H( r). Thus, yĝĝ * dk x dk y (A.14) These equations are very general and should apply to any field distribution and correlation c( r), provided that the amplitude of the heterogeneities is small enough.
Using now the particular distribution field g( r) witĥ where B is a normalisation prefactor determined by σ g = δg 2 . The different correlation functions can be derived after some manipulation: sin 2 θ cos 2 θ (n cos 2 θ + sin 2 θ) 2 dθ = σ 2 It is important to note here that these results are independent of the correlation length λ and the shape of the correlation function. This is due to the fact that the correlation function is istropic, so that the integrals eqs. (A.17-A. 19 can be split into a function depending only on ||k|| multiplied by another depending only on the angular coordinate θ 2 . The first correlation term, Eq. (A.17), is the standard deviation of the velocity field which quantifies the heterogeneity of the velocity field: These expressions can be used in eq. (A.6) to determine the mean pressure gradient G by imposing the mean flow rate U : A similar procedure could be used 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 form u = −D|| ∇P || α−1 ∇P , with α = 1/n and D = C −α . This results to 3 : with D 0 = C −α 0 . Surprisingly, this expression is very similar to eq. (A.21). This similarity in fact originates from a symmetry property of 2D flow fields where the role of pressure and velocity can be switched. This will be demonstrated and discussed below.
The idea consists in rotating the two fields u and ∇P by 90 • . 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 e z : Defining the rotated fields q = e z × ∇P and Z = e z × u, it can be shown that ∇. q = 0 and ∇ × Z = 0. This means that q is a flux vector and Z derives from a potential field Z = ∇Ψ. The two new fields satisfy: Since g(|| u||) = || ∇P ||, it follows || ∇Ψ|| = g −1 (|| q||). Yielding to: As a result, the rotated fields q and ∇Ψ 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: is thus equivalent to: with q c = µ uc κ . Using the relation u c = Aκ γ , the parameters have thus changed to: It is remarkable that for the most natural value γ = 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 withf 0 → −f 0 and the same σ f . It follow that the study can be limited to shear thinning fluids (n < 1) without loss of generality.

Augmented Lagrangian method
This Augmented Lagragian method has been introduced to solve non-Newtonian Stokes equation and has been used by many authors (see for instance [54,55]). In this paper the method was adapted to solve the non-linear Darcy's equation: 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 (in the y−direction). As described above, the solution of such system of equations is equivalent to finding, among all admissible velocity fields the minimum of the functional: The main idea of the augmented Lagrangian method is to introduce a secondary field v in order to decouple the nonlinear problem G(|| u||) from the flow equation. The equality u = v is then guaranteed by the introduction of a Lagragian vector field ξ. An extra term α 2 ( u − v) 2 is also added to enhance the convergence, where α is a small parameter. Another Lagragian field χ is introduced to impose the free divergence, the problem can thus be recast into a saddle point determination: If we now differentiate this functional, we obtain: This set of equations is quite cumbersome. However, the main advantage of this approach lies in the fact that, for a given u and ξ, finding the saddle point for u and χ eqs. (C.6-C.7), is equivalent to solving: with the boundary conditions χ = P in and χ = P out . This is a classical linear Darcy equation with a source term, which can therefore be solved using classical methods. Here, a second order finite difference method has been used. It is important to note that in this equation the permeability is homogeneous and constant, which allows a very fast solution at each step.

Validation
This section presents the validation of the numerical method, in particular, the influence of the mesh size. In this problem, the main characteristic length is the field correlation length λ which determines the amplitude of the velocity and pressure gradients. The numerical resolution is then related to N , the number of mesh nodes per length λ. In this paper, all the simulations were performed with N = 5 and a total system size 1024 × 1024. Fig. C.13.a displays the mean flow rate as function of ∇P , for different N . Although the difference between high and low resolution is not noticeable in this figure, a relative error between low and higher resolution can be defined: Fig. C.13.b represents the evolution of this error as a function of ∇P . The error depends on the flow regime and is, surprisingly, more important in the linear regime than in the non-linear one. The maximum error, however, does never exceed 1%. Figure C.13.c shows the convergence of this error as function of the mesh resolution N and different ∇P . The convergence rate is then slightly faster than N −2 , which is in agreement with the second order finite difference scheme used.

Solutions of eq. (C.12
) For the given fields v and ξ, the alghorithm requires to find u satisfying:  It must be noted that u and B are colinear and with the same orientation because the left term in eq. (C.18) is positive. It is thus sufficient to determine the norm of || u|| = u. The equations then become: where A = µ κu n−1 c . 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.