Shear Rate Determination from Pore-Scale Flow Fields

Aqueous solutions with polymer additives often used to improve the macroscopic sweep efficiency in oil recovery typically exhibit non-Newtonian rheology. In order to predict the Darcy-scale effective viscosity μeff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{\mathrm{eff}} $$\end{document} required for practical applications often, semi-empirical correlations such as the Cannella or Blake–Kozeny correlation are employed. These correlations employ an empirical constant (“C-factor”) that varies over three orders of magnitude with explicit dependency on porosity, permeability, fluid rheology and other parameters. The exact reasons for this dependency are not very well understood. The semi-empirical correlations are derived under the assumption that the porous media can be approximated by a capillary bundle for which exact analytical solutions exist. The effective viscosity μeff(vDarcy)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{\mathrm{eff}} (v_{\mathrm{Darcy}} )$$\end{document} as a function of flow velocity is then approximated by a cross-sectional average of the local flow field resulting in a linear relationship between shear rate γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} and flow velocity. Only with such a linear relationship, the effective viscosity can be expressed as a function of an average flow rate instead of an average shear rate. The local flow field, however, does in general not exhibit such a linear relationship. Particularly for capillary tubes, the velocity is maximum at the center, while the shear rate is maximum at the tube wall indicating that shear rate and flow velocity are rather anti-correlated. The local flow field for a sphere pack is somewhat more compatible with a linear relationship. However, as hydrodynamic flow simulations (using Newtonian fluids for simplicity) performed directly on pore-scale resolved digital images suggest, flow fields for sandstone rock fall between the two limiting cases of capillary tubes and sphere packs and do in general not exhibit a linear relationship between shear rate and flow velocity. This indicates that some of the shortcomings of the semi-empirical correlations originate from the approximation of the shear rate by a linear relationship with the flow velocity which is not very well compatible with flow fields from direct hydrodynamic calculations. The study also indicates that flow fields in 3D rock are not very well represented by capillary tubes.


Introduction
Conventional water floods have an overall recovery efficiency up to 35-40%. Significant amount of oil is left behind in the reservoir partially as residual oil is trapped in the pore space because of capillarity, but also because of bypassing caused by well placement, heterogeneity or viscous fingering instabilities. Improved oil recovery (IOR) and enhanced oil recovery (EOR) approaches aim at improving the overall recovery efficiency (Lake 1989). Many different approaches exist. To give an example, surfactant flooding targets the residual oil by enhancing the microscopic displacement efficiency. Polymer flooding on the other hand mainly targets the bypassed oil by improving the macroscopic sweep efficiency (Reuvers and Golombok 2009;van der Plas and Golombok 2015). There is also an effect on the microscopic displacement efficiency which can have also significant economic benefit. However, the main purpose of the polymer is reduce the mobility of the displacing phase. That is achieved mainly by increasing the viscosity of the aqueous phase, which improves both the mobility ratio but also the sweep because of changed pressure gradients and viscous cross-flow. Furthermore, it suppresses the viscous fingering instability (Saffman and Taylor 1958) which is largely controlled by the shock front mobility ratio (Berg and Ott 2012).
In order to increase the mobility of the displacing aqueous phase, a variety of different additives are considered where hydrosoluble polymer is the most prominent class. More specifically, hydrolyzed polyacrylamide (HPAM) is one of the most suitable polymers considered, and substantial amount of work on that matter has been conducted over the past years. However, there are also biopolymers like Xanthan (Cannella et al. 1988) and other polysaccharides considered. In the framework of Darcy-scale description of multiphase flow in porous media which is typically used for reservoir engineering, the mobility of the aqueous phase is the ratio of the permeability/viscosity. The predominant effect of the polymer to mobility reduction in the aqueous phase is the increase in viscosity, but there is also a reduction in permeability through polymer adsorption and entrapment and a polymer depletion layer causing a slip boundary condition.
The permeability of the porous medium can be significantly reduced by adsorption (Aghabozorgi and Rostami 2016) and entrapment of polymer in the pore space which can lead to plugging of smaller pores and a diversion of the pore-scale flow pattern (Sorbie 1991). A polymer-depleted layer close to the pore wall (Chauvetau 1982;Rodriguez et al. 2014) can cause a slip boundary condition (Berg et al. 2008) which in a Darcy-scale formulation would diminish an effective viscosity. It can also lead to hydrodynamic acceleration of the polymer component. Depending on the polymer formulation, these effects sometimes play a significant role in controlling the mobility reduction factor (Sorbie 1991;Rodriguez et al. 2014). The increase in in situ viscosity is a combination of rheological effects such as shear-dependent fluid viscosity, extensional viscosity (Koroteev et al. 2013;van der Plas and Golombok 2015) and visco-elasticity (Wang et al. 2010). Polymer in general and more specifically aqueous solutions of hydrosoluble polymer including HPAM exhibit a shear rate-dependent bulk viscosity (Delshad et al. 2008). The exact dependency of the shear viscosity μ on the shear rate γ , i.e., μ = μ(γ ) depends on various parameters such as polymer type, concentration, salinity, molecular weight and predominantly the molecular weight distribution, degree of hydrolysis and cross-linking state. For practical applications, where reservoir simulation is performed in order to model the chemical EOR floods in advance, these dependencies need to be implemented into the respective simulator (Delshad et al. 2008).
Several attempts have been made to model polymer rheology. None are fundamentally satisfactory. Therefore, in most practical cases, the μ = μ(γ ) relationship is measured in  1 Simplified illustration showing how shear rate (γ )-dependent shear viscosity (μ = μ(γ )) on the pore scale (left) is upscaled to the effective shear viscosity as a function of Darcy velocity μ eff = μ eff (v Darcy ) on the Darcy scale (middle) which is used in reservoir simulations in a layered formation (right). Note that most geologies are more complex and interconnectivity plays an important role which makes the transition from Darcy to reservoir scale also more complex laboratory tests by a shear rheometer. In porous media flows which exhibit a complex flow field on the pore scale, the shear rate has a spatial variation, and therefore, also the shear rate-dependent bulk viscosity μ = μ(γ ) varies from pore to pore. However, for practical reservoir engineering purpose, usually only the Darcy-scale effective viscosity is considered. Since the shear rate is a pore scale property, the effective viscosity (often termed also in situ rheology or in situ viscosity) is seen as a property depending on the Darcy velocity, i.e., μ eff = μ eff (v Darcy ). This is sketched in Fig. 1.
For most practical purpose, i.e., for practical field applications, typically a specific μ eff = μ eff (v Darcy ) relationship is required for the field of interest. For optimization of the field application, typically polymer type, concentration and molecular weight are varied which results in a large number of polymer formulations which have to be tested. For typical application-relevant polymers including biopolymers such as Xanthan and HPAM, the shear rheology effects dominate below a critical Deborah number (ratio of the stress relaxation time of the polymer and the time scale of the flow) showing shear-thinning behavior over a large range of flow velocities. After initial assessment of the bulk rheology μ = μ(γ ), eventually core flooding experiments are be conducted where the μ eff = μ eff (v Darcy ) relationship is measured in laboratory experiments for the most promising formulations by conducting core floods by varying the flow rate. Since core flooding experiments with hydrosoluble polymer are usually destructive, i.e., polymer remains inside the rock that cannot be removed (unless aggressive cleaning procedures are used that likely also alter the rock itself), and rock samples from the field of interest are often very scarce, the first set of core flooding experiments is conducted on analog rock samples. These are often taken from outcrop rock which has in most cases different porosity and permeability than the rock from the field. In order to obtain the response μ eff = μ eff (v Darcy ) of the respective formulations in the field rock, one would need to either translate the bulk rheology μ = μ(γ ) directly to μ eff = μ eff (v Darcy ) for the rock material of interest or at least translate the μ eff = μ eff (v Darcy ) determined from the analog rock material (which is available in unlimited quantity) to the scarce field cores.
Therefore, the main objective of a large body of research is to establish a relationship between bulk rheology μ = μ(γ ) and effective viscosity μ eff = μ eff (v Darcy ) to be used in Darcy-type flow simulators.
A direct link between the bulk viscosity μ = μ(γ ) and Darcy scale μ eff = μ eff (v Darcy ) can be made by using approximate analytical models (Fayed et al. 2016) or hydrodynamic simulation of the pore-scale flow field for the non-Newtonian rheology (Balhoff 2000;Afshar-poor et al. 2012;Clemens et al. 2012;Afsharpoor and Balhoff 2013;Tosco et al. 2013). Even though these approaches become more and more available, the most commonly used approach are semi-empirical correlations for the Darcy-scale effective shear rate (Delshad et al. 2008). These semi-empirical correlations are based on the Cannella or Blake-Kozeny equation (Meter and Bird 1964;Cannella et al. 1988;Sorbie 1991) which was originally developed for capillary bundles. They are semi-empirical and contain empirical factors such as the shear rate coefficient, often also termed the "C factor", which in the framework of the correlation is a constant. When comparing with experimental measurements, the "Cfactor" depends on parameters like permeability and varies case by case over three orders of magnitude (Wreath et al. 1990). That is not supposed to be because the Cannella equation already contains porosity and permeability as explicit parameters such that C should be independent of permeability. However, the experimental evidence suggests otherwise which makes it difficult to translate the effective viscosity determined experimentally, for instance, in a more readily available model or outcrop rock to the reservoir rock of interest for a field development.
The question is now what causes this explicit dependency of the "C-factor" on porosity and permeability. In general discrepancies between the correlation and experimental data on the Darcy scale can be caused by multiple effects. Polymer-related effects such as adsorption and mechanical entrapment can cause a permeability reduction of the rock. The in situ rheology can be more complex than shear thinning or a shear rheology alone. Also possibly the structure of the flow field and the relation between flow velocity and shear rate could be different from what is assumed in the semi-empirical correlation. For polymer systems, it is very difficult to distinguish the different effects. Therefore, in this work, polymer adsorption and rheology effects are excluded by considering only Newtonian fluids. Pore-scale simulations of Stokes flow on pore structures of sandstone rock and model geometries are used to investigate the relation between flow velocity and shear rate field and assess how that impacts the ability of the Blake-Kozeny correlation to predict shear rate.

The Cannella or Blake-Kozeny Correlation
A central element of this work is to assess the correlation (at same position x) between flow velocity u( x) and shear rate γ ( x) fields, where velocity is a vector and shear rate is derived from a tensor. The main motivation for looking at shear rates originates from polymer flooding. Even though all calculations in this work are performed with Newtonian fluids, for the motivation we have to go back to non-Newtonian fluids as relevant for polymer flooding. Their rheological behavior is typically parameterized as a function of shear rate. Hydrosoluble polymer systems considered for polymer flooding often exhibit rheological behavior with a Newtonian plateau up to a certain shear rate and then a shear-thinning behavior which is often parameterized using for instance the Carreau model (Delshad et al. 2008;Fayed et al. 2016). The shear viscosity μ as a function of shear rate γ is expressed as where μ ∞ the shear viscosity at infinite shear rate, μ 0 the shear viscosity at zero shear share and λ, α and n are material parameters (Delshad et al. 2008). n is the power law exponent in the shear-thinning regime, i.e., n ≤ 1. Equation (1) describes bulk rheology, i.e., continuum mechanics rheology in bulk outside the porous medium. The effective rheology inside a porous medium μ eff on the Darcy scale is interpreted via Darcy's law which relates the Darcy velocity u Darcy (which is essentially a cross section A averaged flux Q of the aqueous phase, i.e., u Darcy = u w = Q/A) to the applied pressure gradient ∇ p via where K is the permeability of the porous medium. In this Darcy-scale picture, the effective shear viscosity μ eff is then related to the (Darcy scale) effective shear rate γ eff using the same functional form of the Carreau model from Eq. (1) but replacing the bulk shear viscosity μ and γ with the effective viscosity μ eff and effective shear rate γ eff . The effective shear rate γ eff is then related to the Darcy flow rate by using semi-empirical correlations (Delshad et al. 2008) based on the Blake-Kozeny or Cannella equation (Sorbie 1991;Cannella et al. 1988) which was originally developed for capillary bundles n is again the power-law exponent of the fluid, u w the Darcy velocity of the water phase, k r,w the water relative permeability, K the absolute permeability, S w the water saturation and φ the porosity. Note that the computation of the effective shear rate in the capillary bundle model from Eq. (4) assumes a no-slip boundary condition (Berg et al. 2008). For each capillary, the shear rate is computed analytically from the parabolic (Poiseuille) velocity profile taking the derivative of the flow velocity perpendicular to the main flow direction [for details see later Eq. (13)]. Slip effects caused by the polymer depletion layer (Rodriguez et al. 2014) are not included. C is an empirical constant which in the original derivation of the capillary-tube-based model is C = 6.0 or "around 6" for Xanthan polymer (Cannella et al. 1988) but can assume other values for different systems. Taking the larger body of literature into account, C is reported to vary between 10 −1 and 10 3 (Wreath et al. 1990). Often a major unknown is the fluid rheology and, contrary to the capillary bundle model, C is not a constant but shows to depend on permeability K and porosity φ (Delshad et al. 2008) and other parameters. That raises the question where the dependency on porosity and permeability originates from and to which extent it is caused by the underlying assumptions in the correlation from Eq. (4).
One potential source is rock heterogeneity. A homogenization approach applied to heterogeneous rock on the Darcy scale showed an explicit dependency on tortuosity T , porosity φ, power law fluid index n, and more general the pore size /permeability distribution and correlations structure of the heterogeneity (Fadili et al. 2002). However, a large variation of the C-factor is also observed for relatively homogeneous rock which still raises the question about its origin.
In order to exclude potential polymer-specific effects such as permeability reduction by polymer adsorption, and more complex rheology on the pore scale such as extensional viscosity (Koroteev et al. 2013), in the following only Newtonian fluids are considered, i.e., the power law index n = 1. The situation is further simplified to single-phase flow, i.e., S w = 1, k r,w = 1. As a consequence, the term lim n→1 3n + 1 4n Combining the factor of 4/ √ 8 = √ 2 = 1.41 and the factor of 0.78 from Eq. (5), we arrive at a total pre-factor around 1.41 × 0.78 ≈ 1.1 which given the range of C-factors from 10 −1 to 10 3 sufficiently close to one (given that C can range over four orders of magnitude). That simplifies Eq. (4) to Estimating u w from Darcy's law, for typical field velocities of 1 ft/day, we obtain for Berea sandstone in the next section effective shear rates γ eff between 2 and 6 s −1 . This effective shear rate is compared with the mean shear rate mean computed from the flow field as explained in the next section.

Pore-Scale Flow Simulation
Single-phase pore-scale Stokes flow simulations are performed using OpenFOAM (Raeini et al. 2014) and the SimpleFFT solver of the GeoDICT package (version 2014; Fraunhofer ITWM, Math2Market GmbH, Kaiserslautern, Germany) (Becker et al. 2008;Berg et al. 2016). In both cases, the simulations were performed for water using a constant viscosity μ = 1 mPa s as a Newtonian fluid in order to exclude uncertainties originating from rheology effects. While the specific solvers and details of the numerical approaches differ, both simulation packages effectively simulate steady-state (fully developed) Stokes flow directly on the pore space of the porous medium, i.e., solve the flow field u for an applied pressure drop p over the porous domain for Newtonian incompressible fluids which can be expressed (neglecting gravity) as where the first equation represents the momentum balance and the second equation incompressibility.
In the momentum balance, the term μ∇ 2 u is already a simplification of the divergence of the viscous stress tensor ∇ · τ The viscous stress tensor τ is a symmetric, second rank tensor which can be expressed for an incompressible Newtonian fluid as product of the shear viscosity μ (which for a Newtonian fluid is constant) and the rate of strain tensor as (Deen 1998) The scalar shear rate γ is represented by the magnitude of the rate of strain tensor which is the key parameter for simple non-Newtonian rheology models such as the Carreau model for non-Newtonian fluids from Eq. (1). The computation of has been performed in Avizo (FEI) computing first the rate of strain tensor using the gradient functionality on each velocity vector component and then its magnitude .

Pore-Scale Flow and Shear Rate Fields Sandstone Rock
In the following, pore-scale flow fields are computed for two different sandstone rocks. The first rock is a Berea sandstone with a porosity of 19.6% and a permeability of 1193 mD in zdirection. The digital image with a domain size of (400) 3 voxels and a resolution of 5.345 µm was taken from the Imperial College website (Raeini et al. 2014). Flow simulations have been performed with OpenFOAM applying a pressure drop of 1 Pa in the z-direction (Raeini et al. 2014).
In order to check how representative the examples from the Berea sandstone are, a second rock sample "RS1" was used which is a sandstone reservoir rock with a porosity of 14.8% and a permeability of 50-200 mD. The digital image has been obtained from X-ray computed micro-tomography at 2.05 µm resolution. The image has been first filtered with a nonlocal means filter and then segmented using a watershed segmentation method using Avizo (Leu et al. 2014). The image was then down sampled by factor of 2 (obtaining an effective resolution of 4.1 µm), and the flow simulations were performed on a 600 × 600 × 950 subdomain with the SimpleFFT solver of GeoDICT by applying a pressure drop of 1 Pa in z-direction.
The properties of the rock and resulting quantities are listed in Table 1. The results of the Berea sandstone and the reservoir rock are qualitatively very similar with very minor quantitative differences originating from the differences in porosity and permeability.
The results of the pore-scale Stokes flow calculations for a (200) 3 voxel subdomain of the Berea sample are displayed in Fig. 2. Panel (A) shows the flow field (i.e., the magnitude u = | u|) and panel (B) the pressure field which has been obtained by numerically solving Eq. (7). The pressure field shows (when averaged over the pore space in the x-y plane) an overall linear decrease from inlet to outlet in z-direction (not shown). Panel (C) shows the shear rate or magnitude of rate of strain tensor computed from the flow field by using Eq. (10). Flow field and shear rate field show a much larger variation with the largest flow velocities and shear rates encountered in pore throats.
From the flow velocity and shear rate fields, respective histograms are computed using the histogram functionality of Avizo. The results are shown in Fig. 3 where on the left hand  side the histogram of the flow velocity magnitude u = | u| is displayed and on the right hand side the histogram of the shear rate magnitude = 1 2 ( : ). In the velocity magnitude histogram in Fig. 3a, the average velocity is close to (the interstitial velocity corresponding to a Darcy velocity of) 1 ft/day indicating that the simulated flow field is relevant for field conditions. In the shear rate histogram in Fig. 3b, the estimate from the Blake-Kozeny correlation from Eq. (6) for C = 6 (Delshad et al. 2008) is at the upper end of the range of shear rates and about a factor of 15 larger than the average shear rate.
The results of the reservoir rock sandstone sample "RS1" are qualitatively very similar to Berea and therefore not displayed. All results are summarized in Table 1. The mean share rates mean are typically a factor 15-32 smaller than the estimates from the Blake-Kozeny correlation in Eq. (6) for C = 6. That in itself is not really surprising because it is well understood from the literature that the "C-factor" can vary by three orders of magnitude (Wreath et al. 1990). By adjusting the "C-factor" from Eq. (6), these can be brought to a match within the range of values reported in the literature. The respective modified C-factor is then termed C' and is listed in the last column of Table 1. The magnitude is close to the value of C = 0.69 as reported by Hirasaki and Pope (1974). For all considered structures, i.e., the two sandstone rocks, capillary bundle and sphere pack, C' is actually very similar, ranging between 0.2 and 0.4. It does not exhibit the large variation over several orders of magnitude as observed in experiments with polymer (Wreath et al. 1990). Possible reasons are that in this study polymer-related and rheology-specific effects and uncertainties are excluded which  Fig. 3 Histograms of the velocity field u = | u| (a) and shear rate = 1 2 ( : ) (b) of the Berea sandstone. In the velocity histogram, the vertical red line represents the interstitial velocity for a 1 ft/day Darcy velocity indicating that the simulated flow field is in the same range of field-relevant flow rates. The vertical gray line represents the average velocity of the flow field. In the shear rate histogram, the gray line represents the average from the histogram, and the vertical red line the estimate from the Blake-Kozeny (or Cannella) correlation from Eq. (6) using C = 6 are present in the experimental studies (Wreath et al. 1990). Another possible reason is that for the considered structures the permeability varies within one order of magnitude, while in natural rock usually several orders of magnitude are encountered. Nevertheless, the data in Table 1 show the smallest C' for the smallest permeability indicating a possible trend of C' with permeability compatible with reports in the literature (Wreath et al. 1990).

Relationship Between Flow Velocity and Shear Rate in Model Geometries
The fundamental question is why the C-factor adjustment is required and whether there is a more fundamental mismatch between the correlation and its underlying assumptions. This is a long standing question that has been articulated already in Wreath et al. (1990). Teeuw and Hesselink (1980) argued that C depends on the ratios of the radii and lengths of the contractions and dilations in the underlying porous media which would in practice have a direct impact on the pore-scale flow field. Therefore, in the following, we inspect the pore-scale flow fields in more detail in particular with respect to the relation between shear rate and flow velocity. For that purpose, 2D correlation histograms between the flow velocity magnitude u = | u| from the flow simulation and the shear rate = 1 2 ( : ) are computed using the 2D histogram functionality of Avizo. We begin with considering the model geometries of sphere packs and capillary bundles which are often used as model systems for porous media. The flow fields have been computed in a similar way as the flow fields for the sandstone rocks. The respective flow and shear rate fields and the shear rate histograms are displayed in Fig. 4. Effective shear rates and the comparison with the correlation are listed in Table 1. The respective 2D shear rate versus flow velocity correlation histograms are displayed in Fig. 5.
The shear rate histogram of the capillary tubes case shown in Fig. 4e extends over a much narrower range and is overall very different from the sandstone case in Fig. 3. Therefore, for comparison, also the model structure of a sphere pack was considered. The respective flow   Table 1. e and f show the respective histograms of the shear rate γ field, shear rate field and shear rate histogram are displayed in Fig. 4b, d, f. For the sphere pack, the shear rate histogram from Fig. 4f is somewhat closer to the sandstone case from Fig. 3. The respective 2D correlation histogram is displayed in Fig. 5b. The 2D correlation histograms in Fig. 5 clearly show that there is no unique linear relationship between flow velocity magnitude and shear rate in particular not for the capillary bundle model. Strictly speaking, the correlation from Eq. (6) actually does not suggest such a linear relationship on the basis of the pore-scale flow field but only on the tube cross-sectional averaged flux and shear rate in a capillary tubes model. Let us therefore revisit the flow field in a single capillary tube with radius R as displayed in Fig. 6. For a single capillary of radius R, the velocity in (z-) flow direction is (Berg et al. 2008) and the flux with the permeability of the tube K = R 2 /8. The local shear rate [which is equivalent to the shear rate magnitude , see Eq. (10)] is then which can be expressed in terms of flux Q from Eq. (12) by eliminating the pressure gradient ∂ p/∂z obtaining which is still a local quantity of position r . The cross-sectional averaged shear rate is then . 6 Geometry of a capillary tube with radius R and constant pressure gradient ∂ p/∂z which is in essence (up to constant pre-factors) the Blake-Kozeny relationship for Newtonian fluids from Eq. (6). However, this also shows that the linear relationship is only obtained for capillary tube cross-sectional averaged flux and shear rate (magnitude). For the local flow field and shear rate field, such a unique linear relationship does not exist. In order to obtain the relationship between local shear rate γ (r ) and local flow velocity u z (r ), one can express γ in terms of u z by eliminating r using Eq. (11) and substitute in Eq. (13), obtaining The structure of Eq. (16) is significantly more complicated than the linear relationship between cross-sectional averaged shear rate and flux from Eq. (15). The 2D correlation histogram for the capillary bundle model displayed in Fig. 5a follows the functional form of Eq. (16) which is superimposed as blue line. Note that deviations of the 2D correlation histogram from the analytical model (blue line) originate from discretization issues in numerical computations and when computing derivatives during postprocessing. In the future, perhaps more robust methods can be developed. The reason why we do not observe a linear relationship between shear rate and flow velocity is that Eq. (16) contains two terms where the first is constant and the second varies with u z (r ) which has a square dependency with the radial position r in the tube, but also contains implicit dependencies on position, viscosity and permeability. An order of magnitude estimate indicates that for the typical range of pore sizes R and flow velocities u z , the ratio of the two terms under the square root can vary between 10 −2 and 10 2 . This means that a large velocity does not necessarily cause a large shear rate and vice versa. According to Eq. (16), the maximum flow velocity and minimum shear rate γ are encountered for r = 0, while according to Eq. (13) the maximum shear rate is obtained for r = R where u z is minimal (i.e., zero). That can be also clearly seen in the flow simulation for a pack of capillary tubes (with radius R = 10 µm) displayed in Fig. 4a, c. Note that as Eq. (16) is derived by substituting Eq. (11) into (13), it leads to the same cross-sectional averaged shear rate as Eq. (15).
The 2D correlation histogram for the sphere pack displayed in Fig. 5b shows clearly a stronger proportionality between flow velocity u and shear rate γ than for the capillary tubes case. For the case of capillary tubes, the root cause for the "anti-"correlation between shear rate and flow field is the no-slip boundary condition at the wall which causes maximum shear rate ∂u z /∂r at zero velocity u z = 0. For the flow field of the sphere pack, the situation is different for the regions between the spheres as indicated in Fig. 7. Figure 7 represents a cross section along the velocity field u z in a plane outside of where the spheres touch. In a (red) line perpendicular to the main flow direction (i.e., in y-direction) at the location where the flow field touches the sphere, we encounter a situation where the flow velocity u z = 0 and the shear rate ∂u z /∂ y is maximum in point (1), similar to the situation in the capillary tube. However, in the regions between spheres (blue lines in Fig. 7), the flow field is stagnant with the flow velocity u z < 3 · 10 −7 m/s ≈ 0 i.e., practically zero, and for symmetry reasons ∂u z /∂ y = 0 in point (2), i.e., velocity and shear rate converge both to practically zero. In the vicinity, a proportionality between shear rate and flow velocity is expected. That is reflected in the 2D correlation histogram in Fig. 5b for up to about half to two thirds of the velocity range. For increasing velocities further away toward point 3 where velocities initially keep increasing, eventually the velocity reaches a maximum value and the shear rate decreases again to zero because point (3) is again a symmetry point. That is also reflected in the 2D correlation histogram in Fig. 5b for the larger velocity range. It appears as if one could decompose the sphere pack's flow field into a stagnant part where flow velocity Fig. 7 Cross section through the flow field in the sphere pack from Fig. 4b (showing a plane with a detail of 3 × 3 spheres). Superimposed in yellow are the streamlines. Where the flow field is in contact with the spheres (red lines) the flow velocity is zero and the shear rate is maximum in point (1), similar to the situation in the capillary tube. However, in the regions between spheres (blue line) the flow velocity u z ≈ 0 and shear rate ∂u z /∂r = 0 at the symmetry point (2) between spheres. In the vicinity, a proportionality of flow velocity and shear rate is expected. Following the blue line, eventually the velocity reaches a maximum value and the shear rate decreases again to zero because point (3) is again a symmetry point and shear rate are proportional, and a channel-like part where flow velocity and shear rate are anti-correlated as in Fig. 5a. One can regard the sphere pack as a network of channels with cross-flow between them (while the capillary bundle has no cross-flow). Depending on the flow direction, some regions in the connections between channels that allow for the crossflow are stagnant. These stagnant regions cause the proportionality between flow velocity and shear rate.
Note that even though the flow field and shear rate-velocity correlation are different, for the average shear rate for the sphere pack, a linear relationship with the Darcy velocity is expected similar as for the capillary tube bundle.

Relationship Between Flow Velocity and Shear Rate in Sandstone Rock
The 2D shear rate-flow velocity correlation histograms for the Berea and reservoir sandstone "RS1" case are displayed in Figs. 8 and 9, respectively.
For both cases, the 2D correlation histograms show a wide spread in the relationship between flow velocity magnitude u = | u| on the horizontal axis and the shear rate = 1 2 ( : ) on the vertical axis, but no unique linear correlation, meaning that for the same Since there is no simple linear relationship between shear rate and flow velocity, this means that computing an effective viscosity from an average shear rate is not the same as computing it from an average flow velocity, i.e., μ eff (γ (v)) = μ eff (γ (v)). An empirical constant can certainly compensate for the mismatch but that constant is porous medium specific and depends on to which degree the relationship between γ and v is linear. In order to assess the sensitivity to pore geometry, in the following we consider two model geometries, i.e., capillary tubes and a sphere pack.
It appears that on the basis of the 2D correlation histograms, the flow fields of capillary tubes and sphere packs from Fig. 5 are limiting cases at opposite ends of the spectrum. On the basis of the 2D correlation histograms from Fig. 5 compared with the flow field in sandstone rock from Figs. 8 and 9, it becomes clear that in terms of flow field, capillaries are not a good representation of rock. The flow field of rock seems to lie between that of the capillary tubes case and the sphere pack, or is a combination thereof. In the discussion of the sphere pack results at the end of the previous section, we already rationalized that the flow field of a 3D structure can be decomposed into more channel-like regions with capillary tube-like flow fields and cross-connections between these channels that are stagnant causing a proportionality between flow velocity and shear rate. The flow field in any 3D structure is then a combination of those two contributions but the exact partition depends on properties like the length of pore throats and the aspect ratio of their diameter compared to pore bodies, and the coordination number.
Elongated pore throats which resemble capillary tubes are expected to have a more channel-like flow field which resembles more that of capillary tubes in Fig. 5a. That is ulti- Fig. 9 2D correlation histogram between velocity magnitude v = | v| on the horizontal axis and shear rate magnitude on the vertical axis for the reservoir sandstone rock "RS1" mately the reason why the degree to which μ eff (γ (v)) = μ eff (γ (v)) depends on the actual pore space morphology, or in other words the "C-factor" varies with the pore structure. In the Berea sandstone case in Fig. 8, it appears as if the flow field is dominated by stagnant flow regions because the shear rate converges to zero for decreasing flow velocities. For the RS1 sandstone, there is a predominant contribution from stagnant flow regions as well but it has also channel-like elements because we observer a significant shear rate at vanishing flow velocity. Therefore, the Berea case is somewhat more compatible with a linear relationship between flow velocity and shear rate than RS1 which is also reflected by the C' in Table 1 being closer to one for the Berea case compared to RS1.
Note that in this study the observed C-factors range between 0.2 and 0.4 (Table 1). The data suggest a weak dependency of the C-factor on permeability compatible with reports in the literature (Wreath et al. 1990). However, in this study the permeability range considered is only roughly one order of magnitude. In order to draw firmer conclusions on the permeability dependency of the C-factor, a variation of permeability by several orders of magnitude (Wreath et al. 1990) which is in various aspects beyond the scope of this study.
This work is limited to Newtonian rheology. When non-Newtonian rheology effects are included in the computation of the flow field, it is expected that the flow field is changing and that consequently the correlations between flow and shear rate are also changing. Compared to the Newtonian case, for a shear-thinning fluid regions with low shear rates (and consequently higher viscosity) are expected to have lower flow velocities and regions with high shear rate (and consequently low viscosity) are expected to have higher flow velocities. Dependent on the extent of the shear thinning, i.e., the relationship between viscosity and shear rate, the stagnant regions may become larger and the higher velocity flow paths more localized compared to the Newtonian case. When including extensional rheology effects, it is expected that regions with large variation in pore diameter such as pore throats would be significantly impacted by higher flow velocities, while the stagnant regions become even more stagnant. It might be useful to consider the case of the capillary tube bundle and sphere packs as reference cases for which exact analytical solutions can be computed for certain rheological models like Carreau fluids using the Weissenberg-Rabinowitsch-Mooney-Schofield integral method (Sochi 2015). But also a 3D porous rock should be considered as this work clearly showed that capillary bundles and sphere packs are only limiting cases for flow fields of sandstone rock. Resulting flow fields can be visualized and categorized by extending the methodology from this work, i.e., visualizing the flow field, shear rate and strain rate in flow direction superimposed over a representative volume of the porous rock, and in addition use the 2D correlation histograms (shear rate vs. flow velocity and strain rate in flow direction vs. flow velocity).

Summary and Conclusions
For practical applications in improved oil recovery, often semi-empirical correlations based on the Cannella or Blake-Kozeny correlation are used. They relate the in situ viscosity (interpreted from the pressure drop at given flow rate from Darcy's law) to the fluid's shear viscosity of, for instance, a hydrosoluble polymer solution through the use of an effective shear rate. This semi-empirical correlation, which is based on a capillary tubes model, is known to require tuning of the phenomenological "C-factor" which is explicitly dependent on porosity, permeability and polymer-specific properties and effects. The explicit dependency of the "C factor" on permeability means that a correlation established for a specific polymer formulation is valid only for one rock type and cannot be generally applied to another rock type which is a major limitation for practical applications. In order to shed more light on the underlying cause, the connection between effective viscosity, shear rate and flow velocity was studied systematically on the basis of pore-scale flow fields from which all quantities of interest can be computed. In order to exclude rheology and polymer-specific effects, only Newtonian fluids were considered.
Pore-scale resolved shear rates directly computed from the flow field ranged from 10 −2 to 10 1 s −1 for a 1 ft/day flow velocity typically used in field applications. The Blake-Kozeny / Cannella correlation predicted a shear rate that can be brought to an agreement with the average obtained from pore-scale flow simulation by adjusting the "C-factor" within the 3-4 orders of magnitude range reported in the literature. That is unsatisfactory when trying to make predictions for an unknown rock type without prior tuning.
The underlying reason is that the relationship between shear rate and velocity is generally more complicated than the linear relationship employed in the Blake-Kozeny relationship. In terms of the characteristics of the flow field, capillary tubes and sphere packs represent limiting cases at opposite ends of the spectrum. For sphere packs, there is a moderate linear correlation between flow velocity and shear rate but for capillary tubes the correlation is particularly weak. In capillary tubes where the flow field can be computed analytically, a linear relationship between flux and shear rate exists only on the basis of cross-sectional averages but not on the basis of the local flow field. In addition, the shear rate is related not only to the flow velocity but also the position and hence parameters like the local pore geometry. That implies that an effective viscosity computed from an average flow velocity (as used in the Cannella equation) is not the same as the average the viscosity computed on the basis of the local shear rate. Sandstone rock which lies between two limiting cases also does not exhibit a particularly strong correlation between local flow velocity and shear rate. As a consequence, an effective viscosity estimated from an average flow velocity is not the same as from an average shear rate. The degree of the discrepancy depends directly on the extent to which the relationship between shear rate and flow velocity is linear. Pore-scale flow simulations conducted on sandstone rock, a sphere pack and a bundle of capillary tubes showed that the relationship between shear rate and flow velocity strongly depends on the morphology of the pore space. That is ultimately the reason why the "C-factor" varies from rock to rock and requires individual tuning.