Pore-scale permeability prediction for Newtonian and non-Newtonian fluids

The flow of fluids through porous media such as groundwater flow or magma migration is a key process in geological sciences. Flow is controlled by the permeability of the rock; thus, an accurate determination and prediction of its value is of crucial importance. For this reason, permeability has been measured across different scales. As laboratory measurements exhibit a range of limitations, the numerical prediction of permeability at conditions where laboratory experiments struggle has become an important method to complement laboratory approaches. At high resolutions, this prediction becomes computationally very expensive, which makes it crucial to develop methods that maximize accuracy. In recent years, the flow of non-Newtonian fluids through porous media has gained additional importance due to, e.g., the use of nanofluids for enhanced oil recovery. Numerical methods to predict fluid flow in these cases are therefore required. Here, we employ the open-source finite difference solver LaMEM (Lithosphere and Mantle Evolution Model) to numerically predict the permeability of porous media at low Reynolds numbers for both Newtonian and non-Newtonian fluids. We employ a stencil rescaling method to better describe the solid–fluid interface. The accuracy of the code is verified by comparing numerical solutions to analytical ones for a set of simplified model setups. Results show that stencil rescaling significantly increases the accuracy at no additional computational cost. Finally, we use our modeling framework to predict the permeability of a Fontainebleau sandstone and demonstrate numerical convergence. Results show very good agreement with experimental estimates as well as with previous studies. We also demonstrate the ability of the code to simulate the flow of power-law fluids through porous media. As in the Newtonian case, results show good agreement with analytical solutions.


Introduction
Fluid flow within rocks is of interest for several Earth Science disciplines including petrology, hydrogeology and petroleum geoscience, as fluid flow is relevant to the understanding of magma flow, groundwater flow, and oil flow respectively (Manwart et al., 2002). Permeability estimates can be inferred on several scales ranging from macroscale (crust) (Fehn and Cathles,20 1979; Norton and Taylor Jr, 1979) over mesoscale (e.g. bore hole) (Brace, 1984) to pore scale (e.g. laboratory) (Brace, 1980). Permeability at crustal scale is of great importance as crustal scale permeability is a function of its complex microstructure, therefore an accurate prediction of permeability on the pore scale is necessary (Mostaghimi et al., 2013). Typical limitations for laboratory measurements on pore scale are: (i) change of the sample's microstructure and therefore its physical properties through cracking and self-filtration (Zeinijahromi et al., 2016;Dikinya et al., 2008) (ii) pressure changes due to the influence of wall effects (Ferland et al., 1996) and finally (iii) difficulties to measure irregular grain shapes and small grain sizes of the porous medium (Cui et al., 2009;Gerke et al., 2015).
At this point numerical modelling can help to compute permeabilities and understand the microstructures as well as flow 5 patterns in three dimensional pore structures. To compute fluid flow directly within 3D pore structures it is necessary to determine the morphology of the investigated sample. This can be achieved by digital rock physics (DRP). It is a powerful tool which allows to improve the understanding of both pore scale processes and rock properties. DRP approaches use 2D or 3D microstructural images to compute fluid flows (Fredrich et al., 1993;Ferreol and Rothman, 1995;Keehm, 2003;Bosl et al., 1998), which are obtained using modern techniques including x-ray computer tomography and magnetic resonance imaging 10 (Dvorkin et al., 2011;Arns et al., 2001;Arns, 2004). In a first step the obtained microstructural images undergo several stages of segmentation (binarization, smoothing etc.) necessary to create a three dimensional pore space. The subsequent computation of fluid flow through the reconstructed three dimensional pore space is tackled with either Lattice-Boltzmann (Bosl et al., 1998;Pan et al., 2004;Guo and Zhao, 2002) , Finite Difference (Manwart et al., 2002;Shabro et al., 2014;Gerke et al., 2018) or Finite Element methods (Garcia et al., 2009;Akanji and Matthai, 2010;Bird et al., 2014). The computed velocity field is then 15 used to estimate permeability (Keehm, 2003;Saxena et al., 2017) and other physical properties (Saxena and Mavko, 2016;Knackstedt et al., 2009).
In recent years, the flow of non-Newtonian fluids has gained significant interest due to their use in a wide range of applications including geology, medicine and other industrial processes (e.g. Johnston et al., 2004;Choi, 2009;Suleimanov et al., 2011;Mader et al., 2013). Nanofluids contain nanometer-sized particles and have been shown to significantly enhance the effi-20 ciency of oil recovery (Wasan and Nikolov, 2003;Huang et al., 2013), whereas the bubbles and/or crystal content of magmas controls their rheology and thus ultimately their eruption style (Mader et al., 2013;Cassidy et al., 2018). If the suspended particles are much smaller than the system to be modeled, the behaviour of these suspensions is commonly described using an effective rheology, exhibiting non-Newtonian behaviour in most cases. For magmas it it not quite clear which physical process is responsible for the non-Newtonian behaviour (Deubelbeiss et al., 2011) as the non-Newtonian behaviour usually originates 25 from the interaction of suspended particles with each other and the surrounding fluid. Therefore, it is necessary to develop numerical models that can simulate non-Newtonian flow through porous media.
In this paper we enhance the open-source finite difference solver LaMEM to model fluid flow on pore-scale with both Newtonian as well as non-Newtonian rheologies. We show that rescaling the staggered grid stencil to better describe velocity components parallel to the fluid-solid interface significantly improves the accuracy. The code is verified using analytical 30 solutions and then used to perform the permeability computations for a digital Fontainebleau sandstone sample (Andrä et al., 2013b).
Fluid flow in porous media can be characterized with the Reynolds number which relates inertial to viscous forces: where ρ is fluid density, v is velocity in direction of the flow, L is the characteristic length and η the fluid viscosity. Due to the small pore size, flows in porous media commonly exhibit small Reynolds numbers and are thus considered to be laminar 5 (Bear, 1988). For geological applications, Reynolds numbers typically are around 10 −9 − 10 −10 for magmas (Glazner, 2014) and range from 10 −8 to 10 −5 for ground water flow. This allows to simplify the incompressible Navier-Stokes equations to the Stokes equations (ignoring gravity): where P denotes pressure, v the velocity component and x the spatial coordinate.
If the pore structure of a porous medium is known, eq.(2) and eq.(3) can be used to directly model laminar fluid flow within this medium. However, at larger scales direct numerical simulation of porous flow is not feasible. In the case of Newtonian fluids, it is common to define a permeability k which relates the flow rate Q to the applied pressure gradient ∆P/L as well as fluid viscosity η: where A is the cross-sectional area of the porous medium. Eq.(4) is also known as Darcy's law and forms the basis of an effective description of Newtonian fluid flow in porous media (Andrä et al., 2013b;Saxena et al., 2017;Bosl et al., 1998). As stated above, this permeability is commonly determined by experimental methods on all scales. With the advent of numerical models for subsurface fluid flow (e.g. FEFLOW (Diersch, 2013)), it has become possible to predict large scale subsurface fluid 20 flow using micro permeabilities as input parameter. Therefore an accurate prediction of micro permeabilies is necessary.
One possibility to do this is to relate the porosity φ of the medium to its permeability k. Deriving the exact nature of this relationship it not trivial and has been subject to a significant amount of research (Kozeny, 1927;Carman, 1937Carman, , 1956Mavko and Nur, 1997). Due to the strong dependency of the permeability not only on porosity, but also on the 3D structure of the pore space, these approaches still suffer from inaccuracies. Due to the development of pore-scale numerical models, it has 25 become possible to determine and refine the porosity-permeability relationship using direct numerical simulation on the basis of computed tomography (CT). These simulations typically provide solutions for fluid velocity v and pressure P for a given pressure gradient across the sample. From the velocity field in zdirection the volume-averaged velocity component v m is calculated (e.g. Osorno et al., 2015): where V f is the volume of the fluid phase. Making use of eq.(4) and Q = v m · A, the intrinsic permeability k s of the sample can then be computed as: As described above, the flow of non-Newtonian fluids through porous media has gained considerable attention in recent years. Here, we use a power law rheology given by: where η 1 and η 2 are the upper and lower cutoff viscosities at the corresponding strain-ratesε 1 andε 2 . η 0 is the fluid viscosity 10 at the reference strain-rateε 0 andε = 1 2ε ijεij the effective strain-rate. n is the power law exponent. With the definition adopted here, fluids with n < 1 are called shear-thinning, while fluids with n = 1 behave as Newtonian fluids and n > 1 are considered shear-thickening fluids. Note that this definition of n differs from the common definition used in geodynamical modelling (called n here), where n = n −1 .
In the case of non-Newtonian fluids, the definition of a permeability is not as straightforward as in the Newtonian case. 15 Several studies have attempted to describe porous media permeability for non-Newtonian fluid rheologies. Until now a general description could not be found as used approaches differ. To develop a nonlinear variant of Darcy's law, Bird et al. (1960) assumed that porous media can be represented by parallel pipes and scaled up these capillary models to general porous media.
By doing so, he suggested that the average velocity v m scales as a function of the driving force F or the pressure gradient ∆P/L (Bird et al., 1960;Larson, 1981): where k is the permeability, η ef f an effective viscosity and K F a related model parameter. If n = 1 and η ef f = η, eq.(4) is recovered. Both the fraction k/η ef f as well as K F depend on porosity φ, stress exponent n, the reference viscosity η 0 and the pore scale geometry of the medium. Consequently, a simple expression for the permeability k has not been found yet.
Attempts to generalize Darcy's law based on eq.(8) include effective medium theories (Sahimi et al., 1990), pore network 25 models (Shah and Yortsos, 1995) and pore-scale numerical simulations (Aharonov and Rothman, 1993;Vakilha and Manzari, 2008). Irrespective of the chosen approach and the exact form of either k/η ef f or K F , eq.(8) implies that a logarithmic plot of v m vs. either ∆P/L or F should produce a straight line with slope 1/n.

Method
We solve the system of governing equations (2) and (3) on a cubic lattice using the finite difference code LaMEM, which 5 has originally been developed to simulate large scale deformation of the Earth's lithosphere and mantle (Kaus et al., 2016).
Here, we will focus on modeling the flow of a fluid with both linear and non-linear viscosity η through a rigid porous matrix.
LaMEM employs a staggered grid finite difference scheme (Harlow and Welch, 1965) to discretize the governing equations  on the data from CT-scans, each cell is assigned either a fluid or a solid phase. The discretized system is then solved using 10 an iterative multigrid scheme to obtain values for velocities v and pressure P . To this end, we employ multigrid solvers which are part of the PETSc library (Balay et al., 2010). As only cells belonging to the fluid phase exhibit non-zero values for the velocity, the velocity components belonging to solid cells are directly set to zero and only considered as boundary conditions. This greatly reduces the degrees of freedom of the system to be solved and hence also the computational cost.
Pressures are fixed on the top and bottom boundaries and free slip boundary conditions are employed on the side boundaries. 15 As described above, no slip boundary conditions apply at the solid-fluid interface. To solve the linear system of equations a V-cycle geometric multiplicative multigrid solver is used (Fedorenko, 1964;Wesseling, 1995). The multigrid solver operates on up to five multigrid levels depending on the given input model. Convergence criteria are given by a relative convergence tolerance of 10 −8 and an absolute convergence tolerance of 10 −10 (see appendix A1). The absolute convergence tolerance atol is defined as the absolute size of the residual norm and rtol the decrease of the residual norm relative to the norm of the right 5 hand side. Therefore convergence at iteration k is reached for: where r k = b − Cx k with b is the right-hand-side vector, x the solution vector of the current timestep k and C the matrix representation of a linear operator (Balay et al., 2010).
Assigning solid and fluid phases to different cells defines the location of the fluid-solid interface. In the case of a staggered When stencils contain rock cells (e.g. fig. 2b) we can straightforwardly enforce the no-flow conditions at their boundaries: to obtain: This form, however, does not enforce interface-parallel velocities to be zero at the interface locations, which results in suboptimal convergence. Alternatively, the exact constraints can be enforced: which will give: The specific expression will depend on the exact subset of cells occupied by rock. The discretization of the other components is performed in a similar manner. The above modification of the shear stress discretization stencil is called here "stencil rescaling". 25 Similar approaches have already been presented in the literature (Vasilyev et al., 2016;Mostaghimi et al., 2013;Manwart et al., 2002). Both Manwart et al. (2002) and Mostaghimi et al. (2013) presented tests to validate their method. The test performed in Manwart et al. (2002) (permeability of a cubic array of spheres) exhibits nonmonotonous convergence of the numerical solution. Mostaghimi et al. (2013) validated their method by comparing the numerical solution to the analytical solution of flow between two parallel plates. They found that they were able to compute the velocity "to within machine accuracy" if they 5 used more than two grid cells, but did not provide any information about convergence of the effective permeability.

Comparison with analytical solutions
To verify the method presented above, we performed a series of benchmark tests where we compared numerical solutions of simplified model setups to their respective analytical solutions. For simplicity, we non-dimensionalized the governing equations (2) and (3) as well as the rheology given in eq.(7) with characteristic values for viscosity η c , length l c , stress τ c and velocity v c : 5 where the characteristic value for v c can be derived from the other characteristic values. Non-dimensional values are denoted with a˜. For the remainder of this section, we will only use non-dimensional values and drop the˜for simplicity. Benchmark tests are organized as follows: first, we will present three benchmark tests for the flow of a Newtonian fluid through i) a single tube, ii) multiple tubes and iii) through a simple cubic sphere pack, which is followed by a benchmark test of power law fluid flow through a single tube. The difference between numerically and analytically computed permeabilities is then expressed 10 using the L 2 norm of their relative misfit: where k comp is the computed and k ana the analytically obtained permeability.

Newtonian flow through a single vertical tube
For a single vertical tube, the analytical solutions for both velocity v and flow rate Q are given as (e.g. Poiseuille, 1846;Landau 15 and Lifshitz, 1987): where ∆P L is the pressure drop in z-direction, R being the radius of the pipe and r the integration variable. The characteristic scales in this case are given by η c = η 0 , τ c = ∆P and l c = R so that the pipe radius R, fluid viscosity η and pressure difference 20 ∆P all take values of 1. The cubic model domain has an edge length of 4 units. Combining eq.(21) with eq.(4), the nondimensional permeability is then given by: permeability was computed using the standard finite difference approach without stencil rescaling. The two sets then differ due to the exact location of the pipe. In set 1, the location of the pipe was chosen in such a way that the pipe surface aligned with the numerical grid (standard, ON NODE) so that computational nodes were directly located on the fluid/solid interface in in 5 xand y-direction. In set 2 (standard, OFF NODE), the location of the pipe was shifted so that the fluid/solid interface was located between the respective computational nodes. The same procedure was applied to sets 3 and 4 where stencil rescaling was employed. The reason to do that was to determine the effect of well-aligned computational nodes, as this is often not the case in more complex geometries.
As expected, the numerical results generally show higher accuracy when stencil rescaling is employed and when node 10 locations and interfaces of the tube are aligned (see Fig.3). The order of convergence is linear for cases without rescaling or when the tube interface does not coincide with grid nodes, but superlinear if both rescaling is employed and interface and node location coincide.

Newtonian flow through multiple vertical tubes
In natural rocks larger channels tend to dominate the overall permeability. To assess this effect, we compute the flow through 15 several straight tubes with different radii (Fig.4). We use four pipes with non-dimensional radii given as R 1 = 2, R 2 = 1, Similar as observed for the single tube setup, the results show a lower relative error for calculations employing the stencil rescaling compared to those without. Furthermore as shown for the setups with single tube the results are more accurate in 25 cases where the numerical grid aligns with the tube surface. As expected, the overall permeability is dominated by the largest tube, as we do not see any significant changes within the relative error of the computed permeability.

Newtonian flow through simple cubic (SC) sphere packs
In order to verify the code for more complex geometries as the vertical tube, we here consider simple cubic (SC) sphere packs.
Sphere packs provide a geometry for different packings as the porous medium is homogeneous. The setup has dimensions of 2 30 in all directions. The permeability of an SC sphere pack is given by (Sangani and Acrivos, 1982;Bear, 1988): where d sp is the sphere diameter and φ is the porosity for simple cubic packing of 1 − π 6 ≈ 0.476, respectively. Fig.5 shows the increase in accuracy with increasing number of grid points employed. The presented relative errors of the permeability value is computed in the same manner as shown in eq.(19). The simulations employing stencil rescaling show a better convergence and seem to saturate against an relative error of 10 −1 , demonstrating the influence of boundary effects through applied no-slip boundary conditions (finite size effect).

Power law fluid flow through a single vertical tube
In order to verify the computed value we compare this setup against an analytical solution of Hagen-Poiseuille flow with power law fluid behaviour. For the single tube configuration described in sec.4.1 and a power law rheology, the velocity within the tube is given by (e.g. Turcotte and Schubert, 2002): where C 1 = 2η − 1 n 0 (see Appendix B), R is the tube radius and r the width of the tube in Cartesian coordinates. Fig.6 shows a good agreement between the numerical and analytical velocities for non-Newtonian fluids when using 0.5 and 2 as values for the power law exponents, covering most fluids used for enhanced oil recovery (e.g. Najafi et al., 2017;Xie et al., 2018).

Application to Fontainebleau sandstone
To verify the ability of the code to handle more complex flows through natural samples and to validate previously computed 5 permeability values we used the CT data for a Fontainebleau sandstone sample provided by Andrä et al. (2013b) with dimen-sions 2.16 mm × 2.16 mm × 2.25 mm (resolved with 288×288×300 grid points). In order to optimize the computation and reduce computational resources a subsample with dimensions of 256 3 is used for further computations. The sample mainly consists of monodisperse quartz sand grains and is therefore a very popular sample for numerical and experimental permeability measurements. Furthermore sandstone is known to be an ideal reservoir rock and is of certain interest for several geological fields, especially in exploration geology. Laboratory measurements of the given sample with porosity ≈ 15.2 % result in a 5 permeability value of ≈ 1100 mD (Keehm, 2003).

Newtonian flow
As in previous tests we compute the permeability of the extracted subsample using eq.(4) and (5). Fig.7a shows streamlines colored using computed fluid velocities and Fig.7b the local pressure. For a resolution of 256 3 , we obtain permeabilities which are comparable to previously computed permeabilities of the same sample ( fig.8, Andrä et al., 2013b), with the rescaled stencil 10 method yielding significantly lower values at higher resolutions. As previous tests show, permeabilities may be overestimated at lower resolutions. To test this effect, we increased the resolution of the Fontainebleau subsample by a factor of 2, 3 and 4 (512 3 ,768 3 ,1024 3 ). The resolution increase is achieved by subdividing a voxel into 2,3 or 4 voxels. We do not apply any interpolation or stochastic reconstructions to conserve spacial statistics as discussed by Karsanina and Gerke (2018). Determining the effects of these more sophisticated methods on computed permeabilities will require further work in the future. Fig.8 shows 15 a comparison between the computed and measured values for the given Fontainebleau dataset. With increasing resolution of the subsample the computed permeability value converges to the laboratory value. In comparison to the initial resolution of 256 3 the computed permeability values decreased by ≈ 24.6 % when using a grid resolution of 1024 3 . Additionally the benefit of stencil rescaling can also be seen here, as e.g. the simulation with a resolution of 512 3 and stencil rescaling predicts nearly the same permeability as the case with doubled resolution and no stencil rescaling. Clearly, the models converge to a value that is close to the measured value. The numerical convergence is computed for several subsamples (see appendix C). Fig. 8 shows the convergence of a single subsample. Previous studies have also observed this convergence with increasing resolution, albeit not always from above (e.g. Zakirov and Galeev, 2019). Similar behaviour has also been observed in LBM simulations (e.g Khirevich et al., 2015;Khirevich and Patzek, 2018).

Power law fluid flow
To demonstrate the capability of the code to compute the flow of non-Newtonian fluids through porous media, we computed the average flow velocity v m for a square subsample of the Fontainebleau sandstone sample described above using the powerlaw rheology given in eq.(7). The edge length of the subsample was 1.92 mm, which corresponds to a CT resolution of 256 3 voxels. To increase accuracy, we increased this resolution by a factor of 2 to a resolution 512 3 . As seen in the section 5.1, 5 results at this resolution may overestimate the actual permeability value. The chosen resolution thus represents a compromise between accuracy and computational cost. The reference viscosity was set to η 0 = 1 Pas and η 1 and η 2 were set to 10 −3 and 10 6 respectively. Two sets of simulations using a power law exponent of 0.5 and 1 were performed. In each set the applied pressure at the top boundary is changed from 1 -16 Pa. In fig.9 we plot the applied pressure at the top boundary against the computed average velocity. For both sets of simulations the computed slopes of (19982 ± 9) × 10 −4 and (1000009 ± 3) × 10 −6 10 are in good agreement with the imposed power law coefficients of 0.5 and 1 (eq.(8)).

Discussion and Conclusion
In this paper, we presented the capability of the open-source finite difference solver LaMEM to compute the permeability of given porous media. The code was verified using a set of benchmark problems with given analytical solutions ranging from Hagen-Poiseuille flow through vertical tubes to more complex flow through simple cubic sphere packs. Using CT Data of a Fontainebleau sandstone, we then demonstrated that the code is able to predict the permeability of natural porous media. In both benchmarks and application tests, the benefits of the stencil rescaling method can be observed, as this method provides significantly more accurate results at no additional computational cost.
Benchmarking results for single and multiple tubes demonstrate that the permeability calculation improves considerably in case the fluid-solid interface and the numerical grid are at least partially aligned. Cases using the stencil rescaling solutions 5 with a velocity change on a computational node produce smaller relative errors.
Similar to studies using the Lattice-Boltzmann method (LBM) (Knackstedt and Zhang, 1994;Zhang et al., 2000;Keehm, 2003) our resolution test for the Fontainebleau subsample shows that the computed permeability value also decreases with increasing grid resolution. For instance, computing the permeability of Fontainebleau sandstone sample with grid resolution of 1024 3 , calculations employing stencil rescaling give approximately the same permeability value as suggested by laboratory 10 measurements, while simulations without employing stencil rescaling overestimate the computed permeability by ≈ 14.72 %.
The computation of permeabilities in a three dimensional pore space using micro-CT data strongly depends on a reasonable quality of the micro-CT images followed by several steps of segmentation in order to resolve tiny fluid pathways. Although 15 high quality input data is required in most cases it is usually computationally expensive to use the entire micro-CT scan with full resolution, thus representative subvolumes or a reduced numerical resolution has to be used as computational resources are limited.
Additionally the segmentation of the CT data has a considerable effect on the computed permeability as discussed in Andrä et al. (2013a), since segmentation of the acquired micro-CT data has a major effect on the three dimensional pore space and 20 therefore on the obtained value. In two phase systems (fluid/solid), segmentation is straightforward whereas it may become more difficult in multiphase systems. All of the above points are a source of uncertainty and need to be considered when comparing numerical calculations to laboratory measurements for rock samples. Furthermore we showed that LaMEM is able to compute non-Newtonian fluid flow in porous media, which is not only relevant for geosciences but also of importance for industrial applications (Saidur et al., 2011). 25 Furthermore it should be kept in mind that solver options like convergence criteria may influence the obtained permeability result. Fig.A1 (see Appendix A) highlights the effect of different relative tolerances on the computed permeability value. In order to avoid spurious results, we recommend to test the influence of the relative and absolute tolerance on the model outcome.
The simulations were performed on the clusters of University of Bayreuth and University of Mainz using different amount of CPUs depending on the size of the computed domain. As an example a simulation with 512 3 voxels uses 1024 CPUs, 30 185 GB RAM and requires 5790 s computation time. Apart from LaMEM finite difference codes like FDMSS  compute permeability of porous media more efficient, but these codes mostly are not able to compute fluid flow using non-linear viscosity.
In conclusion the capability of the open-source finite difference solver LaMEM to accurately simulate Newtonian and non-Newtonian fluid flow in porous media is successfully demonstrated for different setups with an increasing geometric complexity including pipe flow, ordered sphere packs and a micro-CT dataset of Fontainebleau sandstone.
Code availability. https://bitbucket.org/bkaus/lamem/src/master/ ; commit: 9c06e4077439b5492d49d03c27d3a1a5f9b65d32 Appendix A: Convergence criteria 5 To determine whether a numerical solution converges, two convergence criteria are used, which are absolute and relative convergence tolerance. To test the effect on the numerical solution we varied both while computing permeability of three different setups. Our results show that the obtained permeability value saturates for relative convergence tolerances < 10 −7 .
Thus for all further simulations a relative convergence tolerance of 10 −8 is used (Fig.A1). A change in the absolute convergence tolerance did not have any effect on the computed solution, therefore we use a absolute convergence tolerance of 10 −10 .

10
Appendix B: Definition of C 1 C 1 is an constant arising during the derivation of eq. (25), which is related to the nonlinear rheology used in Turcotte and Schubert (2002). This rheology is written as: where n is the stress exponent as used for power law materials in geodynamics. Replacing τ with τ = 2η˙ leads to: Solving eq.B2 for η results in: We can now define a reference viscosity η 0 at a reference strain rateε 0 . This reference viscosity then reads as:  were performed on the btrzx2 cluster, University of Bayreuth and the Mogon II cluster, Johannes Gutenberg University, Mainz. We would like to thank Kirill Gerke and Stéphane Beaussier for their constructive reviews that helped to improve the manuscript.
Fedorenko, R. P.: The speed of convergence of one iterative process, Journal of Computational Mathematics and Mathematical Physics, 4, 227-35, 1964.