Determination of the Effective Viscosity of Non-newtonian Fluids Flowing Through Porous Media

When non-Newtonian fluids flow through porous media, the topology of the pore space leads to a broad range of flow velocities and shear rates. Consequently, the local viscosity of the fluid also varies in space with a non-linear dependence on the Darcy velocity. Therefore, an effective viscosity is usually used to describe the flow at the Darcy scale. For most non-Newtonian flows the rheology of the fluid can be described by a (non linear) function of the shear rate. Current approaches estimate the effective viscosity by first calculating an effective shear rate mainly by adopting a power-law model for the rheology and including an empirical correction factor. In a second step this averaged shear rate is used together with the rheology of the fluid to calculate. In this work, we derive a semi-analytical expression for the local viscosity profile using a Carreau type fluid, which is a more broadly applicable model than the power-law model. By solving the flow equations in a circular cross section of a capillary we are able to calculate the average viscous resistance directly as a spatial average of the local viscosity. This approach circumvents the use of classical capillary bundle models and allows to upscale the viscosity distribution in a pore with a mean pore size to the Darcy scale. Different from commonly used capillary bundle models, the presented approach does neither require tortuosity nor permeability as input parameters. Consequently, our model only uses the characteristic length scale of the porous media and does not require empirical coefficients. The comparison of the proposed model with flow cell experiments conducted in a packed bed of monodisperse spherical beads shows, that our approach performs well by only using the physical rheology of the fluid, the porosity and the estimated mean pore size, without the need to determine an effective shear rate. The good agreement of our model with flow experiments and existing models suggests that mean viscosity is a good estimate for the effective Darcy viscosity providing physical insight into upscaling of non-Newtonian flows in porous media.

When non-Newtonian fluids flow through porous media, the topology of the pore space leads to a broad range of flow velocities and shear rates. Consequently, the local viscosity of the fluid also varies in space with a non-linear dependence on the Darcy velocity. Therefore, an effective viscosity µ eff is usually used to describe the flow at the Darcy scale. For most non-Newtonian flows the rheology of the fluid can be described by a (non linear) function of the shear rate. Current approaches estimate the effective viscosity by first calculating an effective shear rate mainly by adopting a power-law model for the rheology and including an empirical correction factor. In a second step this averaged shear rate is used together with the real rheology of the fluid to calculate µ eff . In this work, we derive a semi-analytical expression for the local viscosity profile using a Carreau type fluid, which is a more broadly applicable model than the power-law model. By solving the flow equations in a circular cross section of a capillary we are able to calculate the average viscous resistance µ directly as a spatial average of the local viscosity. This approach circumvents the use of classical capillary bundle models and allows to upscale the viscosity distribution in a pore with a mean pore size to the Darcy scale. Different from commonly used capillary bundle models, the presented approach does neither require tortuosity nor permeability as input parameters. Consequently, our model only uses the characteristic length scale of the porous media and does not require empirical coefficients. The comparison of the proposed model with flow cell experiments conducted in a packed bed of monodisperse spherical beads shows, that our approach performs well by only using the physical rheology of the fluid, the porosity and the estimated mean pore size, without the need to determine an effective shear rate. The good agreement of our model with flow experiments and existing models suggests that the mean viscosity µ is a good estimate for the effective Darcy viscosity µ eff providing physical insight into upscaling of non-Newtonian flows in porous media.

INTRODUCTION
Flow through porous media is ubiquitous in many natural and industrial systems. Examples include flow through biological tissues, blood vessels and bones [1][2][3] or through soils, sediments and rocks, with long-standing interest in hydrology [4,5], petroleum [6], and chemical engineering [7][8][9]. At low Reynolds numbers (Re≪1) the bulk flow of a Newtonian fluid flowing through porous media is described by Darcy's law where q is the mean flow rate per unit area, also called Darcy velocity and µ is the dynamic viscosity. The variable κ is the permeability and p/L is the pressure drop over the distance L. The proportionality constant K = κ/µ is called hydraulic conductivity and can be derived from Stokes' equation assuming a linear relation between the viscous forces and the flow velocity [10]. While Darcy's law is a sound description for the bulk behavior of a fluid whose viscosity µ is constant, many relevant fluids in e.g., in food [11][12][13] and petroleum [14,15] industry, show a much more complex constitutive law. For most of these socalled non-Newtonian fluids, the viscosity can be described by a nonlinear function of the stress-strain rate tensor E or more specifically its first principal invariantγ = 1 2 √ E : E [16]. Due to the heterogeneity of the flow velocities in the interstitial pore space, shear rates vary considerably inside the porous media. For non-Newtonian flows the coupling of the constitutive equations with the flow field leads to a spatial variable viscous resistance. Consequently, the relation between Darcy velocity and pressure drop cannot be described by a linear function anymore as in the case of Newtonian fluids. In order to obtain a bulk equation for the flow that is linear in the pressure drop, an effective viscosity µ eff -which itself depends on the flow variables-must be used in order to account for the non-linear effects i.e., Here we assumed that the permeability κ is a characteristic constant representing the complexity of the pore space alone. Several empirical and semi-empirical models have been proposed to estimate µ eff [17][18][19][20][21][22][23]. Most of these models start from a capillary bundle representation of the different flow paths through a porous medium and estimate an effective shear ratė γ eff by comparing the flow rate of a power-law fluid with that of a Newtonian Poiseuille flow [24] (see also SI). Although analytical solutions can be derived to determineγ eff for powerlaw rheologies, previous studies proposed various empirical correction factors [19,20] to relate Darcy velocity to the effective shear rate. The effective shear rateγ eff is then inserted into the constitutive law of the fluid of interest µ(γ) to obtain an effective viscosity µ eff . This approach requires an empirical factor to relate q toγ eff , which can vary over several orders of magnitude [25,26], depending on the properties of the fluid, the tortuosity and the permeability. This suggests that the above assumptions are questionable. Additionally most of these models predict a linear relationship between the effective shear rate and the Darcy velocity.
In this manuscript we show that for a Carreau fluid [27], the local viscosity can be derived directly from the fluid's constitutive law and the velocity profile in a mean pore size, using a circular capillary to mimic the flow at pore scale. Contrary to commonly used capillary bundle models, our approach does neither require the knowledge of the tortuosity nor of the permeability. The capillary is only used to calculate a fully developed average flow profile. Finally, we calculate the mean viscous resistance by spatially averaging the local viscosity µ . Comparisons of our results with flow cell experiments and existing models show that µ is a good estimate for µ eff .

Characterization of the Fluid
In order to model the flow of a non-Newtonian fluid, we first need to characterize its constitutive behavior. For most non-Newtonian fluids the constitutive relation between the deviatoric stress tensor T and the applied strain rates E can be described by a time independent scalar function µ = µ(γ), such that T = 2µ(γ)E. Here µ is a generalized viscosity which depends only on the first principal invariantγ = 1 2 √ E : E of the stressstrain rate tensor E [16]. In the case of simple shear flowγ reduces to the shear rate. Many functional forms for µ(γ) have been proposed, where the most common ones are the power-law model (Figure 1A), the Carreau model (Figure 1B), the Cross model or the Herschel-Bulkley model [18,28].
The power-law model is described by where K is the viscosity at the shear rateγ = 1 s −1 and n is the power-law index defining the steepness of the shearthinning decay for n < 1 (see Figure 1). Due to its simplicity, the power-law model is the most commonly used rheology to derive analytical expressions. However, the unbounded powerlaw model has two drawbacks: first the model does not capture the linear shear-strain relation for very low and very high shear rates, that are prevalent in most natural systems, and second the viscosity curve becomes singular in the limit of vanishing shear. Consequently constitutive models which "blend" a powerlaw regime between a Newtonian behavior at low and high shear rates-such as the Carreau model-have been proposed for real world applications. The constitutive law of a Carreau fluid is parametrized by where n is the power-law exponent, µ 0 and µ ∞ are the limits of the viscosity at zero and infinite shear and λ is the reciprocal of the critical shear rate, which describes the onset of the shear thinning regime. The parameter α describes how smoothly the Newtonian regime blends into the power-law.

Current Models
Most commonly applied models to estimate µ eff can be derived by equating the flow rate of a Poiseuille flow [29] with the flow rate of a power-law fluid [30] Q Poiseuille = Q power−law .
For a circular capillary with radius R one obtains π 8µ where we used which is also known as the Rabinowitsch equation [30] to describe the flow rate of a power-law fluid in a capillary. Solving Equation (6) for µ, a power-law viscosity µ power−law can be defined as Equation (8) can be simplified to This power-law viscosity corresponds to the viscosity of a Newtonian fluid which would have given the same pressure drop p/L along a capillary. For a power-law constitutive relation µ = Kγ n−1 Equation (9) can be inverted to obtain an effective shear rateγ ef Using the Rabinowitsch equation, we can express the term where q cap is the mean capillary velocity defined as Q power−law /(πR 2 ). Furthermore, the mean capillary velocity q cap can be defined as the Darcy velocity divided by the porosity, q cap = q . Following Savins [31], the radius R eq can be expressed by where ζ is the tortuosity, R eq is the radius of a capillary in the capillary bundle model (see detailed derivation in Supplementary Information). Inserting R eq into Equation (10) then yieldsγ Empirically it has been found by Cannella et al. [19] that the factor 1/ √ ζ does not fit realistic data and replaced the term 1/ √ ζ by a constant C, i.e., Hirasaki and Pope [20] proposed to use C = 1/ √ 25/12 ≈ 0.69 by using the tortuosity ζ of packed spheres, which has been widely reported to be 25/12 [32,33]. Ignoring the tortuosity ζ , Cannella et al. [19] found a factor of C = 6 to be suitable to describe many flows in different settings. Additionally, Cannella et al. accounted for unsaturated and multiphase flows by correcting the permeability κ to κ r,w κ and the porosity to S w . Here, κ r,w is the relative permeability and S w is the saturation. Consequently, the effective shear rate according to Cannella et al. [19,26] is given bẏ Cannella et al. then used this effective shear rate together with a constitutive law µ(γ) to calculate an effective viscosity. For this last step, mostly the Carreau model has been used due to its ability to fit a wide range of different rheologies.
Other models, that have been developed, are using more complex rheological descriptions of the fluid. Nevertheless, they generally require to correct the analytical solution with empirical factors to achieve reasonable agreement with experimental data [17].

Average Viscosity Approach
Here we present a new approach to estimate µ eff by solving directly for the viscosity profile of a fully developed Carreau flow inside a single capillary of radius R that mimics a mean pore with a mean flow rate of q . This approach assumes that a single constitutive law can be used both on the pore as well as on Darcy scale. This assumption is supported by the observation that even at pore scale, the viscosity distribution covers the whole range of viscosities given by the Carreau model equation. Consequently a single power-law is insufficient to describe the transitional behavior at the low and high shear limits.The onset of the powerlaw regime occurs at a characteristic combination of the Darcy velocity and the pore size. Therefore, it is important to determine at which Darcy velocity the onset of the non-linearities of the fluid starts to matter for a given pore size. The Carreau model allows to obtain an average viscosity profile in a pore without invoking an effective shear rate and an intermediate power-law rheology. Note that the Carreau model includes the critical shear rate 1/λ, which defines the onset of the power-law regime, see Figure 1. To model the flow of a Carreau fluid we perform the following steps: (i) We estimate the characteristic pore size in order to set the diameter of the average pore for which we examine the flow profile. This characteristic pore size can be readily obtained from an pore-or grain-size distribution. (ii) We then compute the velocity profile as a function of the flow rate. (iii) The velocity profile obtained in the previous step can subsequently be used to determine local shear ratesγ(r) = du(r)/dr. (iv) Combining the shear and the Carreau rheology (Equation 4) we obtain the local viscosity distribution µ(γ(r)) in a cross section of the capillary.
(v) Finally, we use the local viscosity, to estimate the effective viscosity µ eff by averaging the viscosity profile over the cross section of the capillary.
In order to apply this concept of a mean profile in a pore, we use a capillary with a circular cross section and assume a fully developed flow profile. The steady state Navier-Stokes equation at low Reynolds numbers in a circular capillary can be written as where the pressure gradient along the capillary is constant dp dx = const. Integrating with respect to r yields µ du dr du dr = r 2 dp dx Based on the symmetry of the flow profile, the velocity is maximal along the center line of the capillary at r = 0. By definition of a maximum, the shear rateγ = du dr vanishes, du dr r=0 = 0, which results in K 1 = 0. Thus, Equation (17) simplifies to µ(γ)γ = r 2 dp dx .
Since the pressure drop dp dx along the capillary is assumed to be constant, we can replace it by a reference pressure p ref = 0 divided by a reference length scale. Choosing the capillary radius R as characteristic length, we define Then, we insert the reference pressure from Equation (19) into Equation (18) and obtain Inserting the constitutive law of a Carreau fluid (Equation 4) into (Equation 20) and solving for r yields This expression can be rewritten by using the boundary condition for the shear rate (γ| r=R = −γ w ), whereγ w is the shear rate at the wall of the capillary. Consequently, the reference pressure is given by the following equation: In order to obtain an expression for the flow profile u(r) andγ w , we integrate r(γ ′ ) given by Equation (21) radially with respect toγ ′ . As the shear rateγ ′ (r) is an odd function, the relation −γ 0 dγ ′ = − γ 0 dγ ′ holds for all r. Consequently, the shear rateγ will be our free parameter ranging from 0 to −γ w .
The resulting integral can be expressed as Here 2 F 1 (a, b; c; z) is a hypergeometric function with parameters a, b and c. Further details about the hypergeometric function can be found in the Handbook of mathematical functions by M. Abramowitz and I.A. Stegun [34]. Using the chain rule on d(rγ ′ ), we can rewrite r(γ ′ )dγ ′ as d(rγ ′ ) −γ ′ dr, which yields Frontiers in Physics | www.frontiersin.org where we already substitutedγ ′ dr with du. Performing now the integration with respect toγ ′ gives Setting Equation (23) equal to Equation (25) allows to solve for the velocity u inside the capillary (Equation (26)).
Applying the non-slip boundary condition at the wall for velocity ( u| r=R = 0) and shear (γ| r=R = −γ w ) finally yields u max as defined in Equation (27).
Therefore, for any given maximum velocity u max at the center of a capillary, the two Equations (27) and (22) can be solved for the two unknown parameters p ref andγ w using a non-linear root finding algorithm. These two parameters can then be inserted into Equations (26) and (21) to obtain the velocity field u and the corresponding radial coordinate r. The remaining free parameter is the shear rateγ, which varies between −γ w and 0. Consequently, one can calculate the flow profile u(r) in a capillary of radius R for any u max . The mean capillary velocity q profile is then readily calculated by integrating u(r) over the capillary cross section , i.e., To obtain the Darcy velocity, the mean capillary velocity q profile has to be multiplied by the porosity , namely, q = q profile . Further, we can determine the shear rate profile from u(r) by differentiation usingγ(r) = du dr . In combination with the Carreau model (Equation 4), the shear rate can be used to calculate the local viscosity µ(r) in a cross section of the capillary from which we can infer the spatial average of µ as Based on this formalism, we now propose that the effective viscosity can be appropriately estimated directly from the average viscosity µ without using an effective shear rateγ eff . Our approach is purely based on the physical constitutive law of the fluid (here represented by a Carreau rheology) and the solution of the momentum equation of such a non-Newtonian fluid in a circular capillary. In order to test the hypothesis that µ eff is appropriately described by µ , we benchmark our approach with flow experiments and compare our model's prediction with that of Cannella et al. [19] and Hirasaki & Pope [20].

Experimental Setup
To measure the effective viscosity of a non-Newtonian fluid flowing through porous media we set up a Darcy experiment in a column (405 mm height; 50.3 mm diameter) of monodisperse glass beads with a diameter of 8 mm (Figure 2). A constant static pressure was applied using a Boyle-Mariotte bottle connected to a pressure controller (Fluigent MFCS TM -EZ). The flow rate was determined by measuring the weight of the outflowing fluid over time, while the pressure drop was measured with two pressure sensors (Keller, PAA-36XW) connected to the Darcy column. The porosity of the packed bed was = 0.4. Before determining the effective viscosity of the non-Newtonian solution we measured the permeability κ of the porous medium using a Newtonian solution (66.7 vol% Glycerol + 33.3 vol% D.I. water) with a well specified viscosity µ Newt = 0.0266 Pa · s [35]. We then measured the pressure drop over the length L = 300 mm at different flow rates, obtaining five measurements for q( P/L). Fitting a straight line to the measured data yields κ = 4.4 · 10 −8 m 2 (Figure 3).
As a non-Newtonian fluid, we used a solution of 0.05 wt% xanthan gum produced by Sigma Aldrich. The rheology curve of the xanthan gum solution was measured with an Anton Paar MCR 702 rheometer with a double gap DG 26.7 geometry. Figure 4 shows the dynamic viscosity of our solution as a function of the applied shear rate together with the fitted Carreau model with parameters µ 0 = 0.085 Pa · s, µ ∞ = 0.001 Pa · s, λ = 2 s, n = 0.48 and α = 0.8 (Equation 4). Very high shear rates have not been measured, but literature values [36] indicate that the viscosity of xanthan gum solutions approach the viscosity of water at high shear rates.
After characterizing the non-Newtonian fluid's constitutive behavior, we perform the same Darcy experiment for the non-Newtonian solution as in the Newtonian case, measuring the pressure drop at 8 different flow rates. The effective viscosity µ eff of the non-Newtonian fluid was then computed from the measured Darcy velocity q, the measured pressure drop p/L and the previously measured permeability κ according to

RESULTS AND DISCUSSION
Before discussing the upscaling of our capillary model and the averaging of the viscosity approach we analyze how the profiles of velocity, shear and ultimately viscosity behave inside a single capillary. Figure 5 shows the velocity, the shear rate and the viscosity profiles, normalized by the respective means, for a capillary of radius 1 mm for two flow rates q profile = 0.003 mm/s and q profile = 3.7 mm/s. For the rheology of the fluid we use the Careau model with parameters µ 0 = 0.085 Pa · s, µ ∞ = 0.001 Pa · s, λ = 2 s, n = 0.48 and α = 0.8. The profiles have been calculated numerically by solving Equations (22) and (27) to determine p ref andγ w together with Equations (21) and (26) as described above.
For the low flow rate q profile = 0.003 mm/s, the shear rate does not exceed the critical shear rateγ crit = 1/λ at any radial position in the capillary. Consequently the viscosity is almost constant at µ 0 resulting in a Newtonian flow behavior and a parabolic velocity profile. Here the shear rate increases linearly fromγ = 0 at the center of the capillary toγ =γ w at the wall. As soon as the flow rate exceeds a certain threshold, nonlinear effects become important and the velocity profile flattens at the center (c.f. purple line in Figure 5A). This nonlinear behavior can also be observed in the shear ratesγ = du/dr, as well as in the local viscosity µ(r) which starts to develop a very distinct maximum at the center of the capillary, Figure 5C. As a result, first averaging the shear rate and then calculating the viscosity leads to a significantly different results than calculating the average viscosity itself, namely µ( γ ) = µ(γ) . The models of Cannella et al. and Hirasaki & Pope assume a pure powerlaw rheology for the calculations of the effective shear rateγ eff at all flow rates (Equations 7, 11). This does not agree with the bahavior shown in Figures 5A-C, which indicates that the profile changes from a Newtonian to a nonlinear behavior depending on the applied flow rate.
After discussing the problem of averaging shear and viscosity in a single capillary, we now benchmark the average viscosity approach. We compare our model estimates with our experimental results and the predictions of the models by Cannella et al. and Hirasaki & Pope. Figure 6 shows the measurement of the effective viscosity as a function of the Darcy velocity using Equation (30)    fits the experimental data equally well as the average viscosity model. Nevertheless, all three models agree in terms of slope in the shear thinning power-law regime. This is consistent with the observation of Teeuw & Hesselink [22], who found that the power-law exponent of the microscopic rheology equals the exponent of the effective viscosity µ eff as function of the Darcy velocity q [19]. The major difference in the behavior of our model compared to the approach of Hirasaki & Pope [20] is found at the transition from the low shear Newtonian regime to the power-law regime. Here the former model predicts a slightly smoother transition than the average viscosity formulation. However, the error bars shown in the inset of Figure 6 reveal that this region also has the highest uncertainty due to the noise in the pressure measurements. Consequently, we cannot conclude that the proposed model, using an average viscosity, captures the transition between the low shear Newtonian regime FIGURE 7 | Characteristic pore size with radius R between three glass beads with radius d/2 (blue circles). and the power-law regime better than the previous models. Nevertheless, in favor of the proposed model, the linear relation in Equation (13) assumed forγ eff ∼ q used by both Cannella et al. and Hirasaki & Pope does not reflect the non-linear behavior of the profiles shown in Figure 5 for a single capillary.
A major challenge of models based on effective shear rates is the need of an empirical factor to computeγ eff . This, so-called C-factor, can vary over several orders of magnitude [25,26]. This problem is especially evident considering the fact that the shape factor varies for different fluids and pore geometries [26]. While the shape-factor C is used to fit the onset of the power-law regime in µ eff against Darcy velocity plot, our model describes this transition intrinsically through the Carreau parameter 1/λ.
Comparing the models of Cannella et al. and Hirasaki & Pope reveals that the empirical shape factor C of Cannella et al.'s model corresponds to C = 1/ √ ζ in the model of Hirasaki & Pope. Using the tortuosity of ζ = 25/12 for a packed bed of uniform spheres [32,33] Hirasaki & Pope arrive at a shape factor of C = 0.69. However, Cannella et al. later found empirically that this tortuosity dependence does not hold in many different settings and proposed C = 6 as a generally good value. However if we calculate the tortuosity associated with this empirical value we obtain a tortuosity ζ = 0.03, which is inconsistent with the geometric interpretation of ζ as the elongation of flow paths in a capillary bundle model. Another problematic point of the models of Hirasaki & Pope and Cannella et al. is that the tortuosity should not matter at all for the determination of the effective viscosity as the viscosity profile does not change along the capillary. Hence the viscosity, respectively the effective shear rate, should not change with the elongation of the flow paths contrary to Equation (13).
In contrast to most currently used capillary bundle models which try to mimic the flow paths inside porous media, our average viscosity approach uses a capillary only to calculate a fully developed flow profile in a single pore. Generally, we find good agreement with experimental measurements and with previous models, provided that the free parameters have been determined from experimental data [16,19,32,33,[37][38][39]. Averaging the viscosity over the mean flow profile has several advantages. First of all, the model only includes the physical rheology of the fluid and does not need a singular power-law model as an intermediate step to calculate an effective shear rateγ eff . This eliminates the problem that arises from µ( γ ) = µ(γ) . Using the mean viscosity to describe the flow resistance of the fluid is physically more intuitive than using an effective quantity derived from comparisons of non-Newtonian flows with Newtonian behavior. Additionally our model only requires a characteristic geometric length factor rather than an empirical shape factor which depends on geometric properties as well as on the fluid. Our characteristic length scale is a property of the porous medium alone and can be estimated independently from permeability, porosity or tortuosity using e.g., pore or grain size distributions. Furthermore, our approach includes the parameter λ, which describes the onset of the shear thinning regime and is given by the rheology. Consequently, the average viscosity approach provides a consistent nonlinear upscaling of a flow profile in a pore with mean pore size and does not require tortuosity [40].

CONCLUSIONS
In summary, we presented a new approach to extend Darcy's law to Carreau fluids using the mean viscosity over a representative capillary as the effective flow resistance. The major advantage of the new model is that it does not require an intermediate effective shear rate and calculates the average viscosity directly using the microscopic constitutive law of the fluid. This approach allows us to upscale the average viscosity of a single pore with a mean pore size to a Darcy scale. This procedure also avoids an empirical shape parameter, which has been replaced by a characteristic length scale derived from the physical property of the porous medium itself. Furthermore, the proposed approach does not require commonly made capillary bundle assumptions like the elongation of a capillary by using the tortuosity and it does not require the permeability [40] to estimate an equivalent capillary diameter. Experimental measurements for a flow through a packed bed of monodisperse beads and comparison with other non-Newtonian capillary models reveal that the average viscosity provides a robust estimate for the effective Darcy viscosity µ eff .

AUTHOR CONTRIBUTIONS
UE designed the research, carried out the experiment, performed the calculations, and analyzed the data. HS helped with the mathematical derivations and the analysis. MF assisted with the experimental setup and PB measured the rheology. JJ-M and JA supported by discussions and consulting. UE wrote the paper with input from all co-authors. MH supervised the research.