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.

1 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, 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) a 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 modeling can help to compute permeabilities and understand the microstructures as well as flow patterns in three-dimensional pore structures. To compute fluid flow directly within 3-D 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 2-D or 3-D 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 computed tomography (CT) and nuclear magnetic resonance imaging (NMR) (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 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 efficiency of oil recovery (Wasan and Nikolov, 2003; Huang et al., 2013), whereas the bubbles and/or crystal content of magmas control 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 behavior of these suspensions is commonly described using an effective rheology, exhibiting non-Newtonian behavior in most cases. For magmas, it is not quite clear which physical process is responsible for the non-Newtonian behavior (Deubelbeiss et al., 2011) as the non-Newtonian behavior usually originates 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 (Lithosphere and Mantle Evolution Model) to model fluid flow on the pore scale with both Newtonian and 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 solutions and then used to perform the permeability computations for a digital Fontainebleau sandstone sample (Andrä et al., 2013b).

2 Fluid flow in porous media

Fluid flow in porous media can be characterized with the Reynolds number which relates inertial to viscous forces:

$$\begin{array}{}\text{(1)}& \mathit{\text{Re}}={\displaystyle \frac{\mathit{\rho }vL}{\mathit{\eta }}},\end{array}$$

where ρ is fluid density, v is velocity in direction of the flow, L is the characteristic length, and η is 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 (Bear, 1988). For geological applications, Reynolds numbers typically are around 10⁻⁹ to 10⁻¹⁰ for magmas (Glazner, 2014) and range from 10⁻⁸ to 10⁻⁵ for groundwater flow. This allows to simplify the incompressible Navier–Stokes equations to the Stokes equations (ignoring gravity):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial v_i}{\partial x_i}}& {\displaystyle }=\mathrm{0}\text{(3)}& {\displaystyle \frac{\partial }{\partial x_j}}\left[\mathit{\eta }\left({\displaystyle \frac{\partial v_i}{\partial x_j}}+{\displaystyle \frac{\partial v_j}{\partial x_i}}\right)\right]-{\displaystyle \frac{\partial P}{\partial x_i}}& {\displaystyle }=\mathrm{0},\end{array}$$

where P denotes pressure, v the velocity component and x the spatial coordinate.

If the pore structure of a porous medium is known, Eqs. (2) and (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 η:

$$\begin{array}{}\text{(4)}& k=-{\displaystyle \frac{\mathit{\eta }LQ}{\mathrm{\Delta }PA}},\end{array}$$

where A is the cross-sectional area of the porous medium. Equation (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 flow using micro-permeabilities as input parameters. Therefore, an accurate prediction of micro-permeabilities 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 is not trivial and has been subject to a significant amount of research (Kozeny, 1927; Carman, 1937, 1956; Mavko and Nur, 1997). Due to the strong dependency of the permeability not only on porosity but also on the 3-D structure of the pore space, these approaches still suffer from inaccuracies. Due to the development of pore-scale numerical models, it has 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 the z direction, the volume-averaged velocity component v_m is calculated (e.g., Osorno et al., 2015):

$$\begin{array}{}\text{(5)}& v_{\mathrm{m}}={\displaystyle \frac{\mathrm{1}}{V_{\mathrm{f}}}}\underset{V_{\mathrm{f}}}{\int }\left|v_z\right|dv,\end{array}$$

where V_f is the volume of the fluid phase. Making use of Eq. (4) and $Q=v_{\mathrm{m}}\cdot A$, the intrinsic permeability k_s of the sample can then be computed as

$$\begin{array}{}\text{(6)}& k_{\mathrm{s}}={\displaystyle \frac{\mathit{\eta }{v}_{\mathrm{m}}L}{\mathrm{\Delta }P}}.\end{array}$$

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

$$\begin{array}{}\text{(7)}& \mathit{\eta }=\left\{\begin{array}{ll}{\mathit{\eta }}_{\mathrm{1}},& \text{if}\dot{\mathit{\epsilon }}<{\dot{\mathit{\epsilon }}}_{\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{0}}{\left(\frac{\dot{\mathit{\epsilon }}}{{\dot{\mathit{\epsilon }}}_{\mathrm{0}}}\right)}^{n-\mathrm{1}}\\ {\mathit{\eta }}_{\mathrm{2}},& \text{if}\dot{\mathit{\epsilon }}>{\dot{\mathit{\epsilon }}}_{\mathrm{2}},\end{array}\right.\end{array}$$

where η₁ and η₂ are the upper and lower cutoff viscosities at the corresponding strain rates (${\dot{\mathit{\epsilon }}}_{\mathrm{1}}$ and ${\dot{\mathit{\epsilon }}}_{\mathrm{2}}$). η₀ is the fluid viscosity at the reference strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}$, and $\dot{\mathit{\epsilon }}=\sqrt{\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{\epsilon }}}_{ij}{\dot{\mathit{\epsilon }}}_{ij}}$ is 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 modeling (called n^′ here), where $n^{\prime }=n^{-\mathrm{1}}$.

In the case of non-Newtonian fluids, the definition of a permeability is not as straightforward as in the Newtonian case. 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 non-linear 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, Bird et al. (1960) 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):

$$\begin{array}{}\text{(8)}& v_{\mathrm{m}}={\left({\displaystyle \frac{k}{{\mathit{\eta }}_{\mathrm{eff}}}}{\displaystyle \frac{\mathrm{\Delta }P}{L}}\right)}^{\frac{\mathrm{1}}{n}}=K_F{\left(F\right)}^{\frac{\mathrm{1}}{n}},\end{array}$$

where k is the permeability, η_eff an effective viscosity and K_F a related model parameter. If n=1 and η_eff=η, Eq. (4) is recovered. Both the fraction k∕η_eff and K_F depend on porosity ϕ, stress exponent n, the reference viscosity η₀ 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 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∕η_eff 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.

3 Method

We solve the system of governing Eqs. (2) and (3) on a cubic lattice using the finite difference code LaMEM, which 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 (Fig. 1).

Figure 1 Staggered grid and location of variables.

Pressures are defined in the middle of the staggered grid cell, whereas velocities are defined on cell faces. Based on the data from CT scans, each cell is assigned either a fluid or a solid phase. The discretized system is then solved using 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 Portable, Extensible Toolkit for Scientific Computation (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. 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⁻⁸ and an absolute convergence tolerance of 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-hand side. Therefore, convergence at iteration k is reached for

$$\begin{array}{}\text{(9)}& \left|\right|r_k|{|}_{\mathrm{2}}<max\left(\mathrm{rtol}\cdot \left|\right|b{|}_{\mathrm{2}},\mathrm{atol}\right),\end{array}$$

where $r_k=b-C{x}_k$, with b the right-hand-side vector, x the solution vector of the current time step 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 grid, the location of the interface therefore does not correspond to the location of the interface-parallel velocity component. To illustrate this issue, the discretization stencil of a shear stress component τ_xy is shown in Fig. 2. When no interfaces are present (Fig. 2a), the finite different discretization results in the following expression (k index is omitted for brevity):

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{xy(i,j)}=\mathit{\eta }{\displaystyle \frac{v_{x(i,j)}-v_{x(i,j-\mathrm{1})}}{\mathrm{\Delta }y}}+\mathit{\eta }{\displaystyle \frac{v_{y(i,j)}-v_{y(i-\mathrm{1},j)}}{\mathrm{\Delta }x}}.\end{array}$$

When stencils contain rock cells (e.g., Fig. 2b), we can straightforwardly enforce the no-flow conditions at their boundaries:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& v_{x(i,j-\mathrm{1})}=\mathrm{0},\\ & v_{y(i-\mathrm{1},j)}=\mathrm{0},\end{array}\end{array}$$

to obtain

$$\begin{array}{}\text{(12)}& {\mathit{\tau }}_{xy(i,j)}=\mathit{\eta }{\displaystyle \frac{v_{x(i,j)}}{\mathrm{\Delta }y}}+\mathit{\eta }{\displaystyle \frac{v_{y(i,j)}}{\mathrm{\Delta }x}}.\end{array}$$

This form, however, does not enforce interface-parallel velocities to be zero at the interface locations, which results in sub-optimal convergence. Alternatively, the exact constraints can be enforced:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\displaystyle \frac{v_{x(i,j)}+v_{x(i,j-\mathrm{1})}}{\mathrm{2}}}=\mathrm{0},\\ & {\displaystyle \frac{v_{y(i,j)}+v_{y(i-\mathrm{1},j)}}{\mathrm{2}}}=\mathrm{0},\end{array}\end{array}$$

which will give

$$\begin{array}{}\text{(14)}& {\mathit{\tau }}_{xy(i,j)}=\mathrm{2}\mathit{\eta }{\displaystyle \frac{v_{x(i,j)}}{\mathrm{\Delta }y}}+\mathrm{2}\mathit{\eta }{\displaystyle \frac{v_{y(i,j)}}{\mathrm{\Delta }x}}.\end{array}$$

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”. 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 non-monotonous 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 used more than two grid cells but did not provide any information about convergence of the effective permeability.

Figure 2 Staggered-grid stencil rescaling. (a) Stencil without interfaces and (b) with solid interfaces. The sketch is based on a cross section of Fig. 1 in the x–y plane.

4 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 Eqs. (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:

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathit{\eta }& {\displaystyle }={\mathit{\eta }}_{\mathrm{c}}\cdot \stackrel{\mathrm{̃}}{\mathit{\eta }}\text{(16)}& {\displaystyle }{x}_i& {\displaystyle }=l_{\mathrm{c}}\cdot {\stackrel{\mathrm{̃}}{x}}_i\text{(17)}& {\displaystyle }\mathit{\tau }& {\displaystyle }={\mathit{\tau }}_{\mathrm{c}}\cdot \stackrel{\mathrm{̃}}{\mathit{\tau }}\text{(18)}& {\displaystyle }{v}_i& {\displaystyle }=v_{\mathrm{c}}\cdot {\stackrel{\mathrm{̃}}{v}}_i={\displaystyle \frac{l_{\mathrm{c}}{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\eta }}_{\mathrm{c}}}}\cdot {\stackrel{\mathrm{̃}}{v}}_i,\end{array}$$

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 using the L₂ norm of their relative misfit:

$$\begin{array}{}\text{(19)}& \left|\right|{\mathit{\delta }}_k|{|}_{\mathrm{2}}=\sqrt{{\left({\displaystyle \frac{k_{\mathrm{comp}}-k_{\mathrm{ana}}}{k_{\mathrm{ana}}}}\right)}^{\mathrm{2}}},\end{array}$$

where k_comp is the computed and k_ana the analytically obtained permeability.

4.1 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 and Lifshitz, 1987)

$$\begin{array}{}\text{(20)}& {\displaystyle }v={\displaystyle \frac{\mathrm{\Delta }P}{\mathrm{4}\mathit{\eta }L}}\left(R^{\mathrm{2}}-r^{\mathrm{2}}\right)\text{(21)}& {\displaystyle }Q={\displaystyle \frac{\mathit{\pi }\mathrm{\Delta }P}{\mathrm{8}\mathit{\eta }L}}{R}^{\mathrm{4}},\end{array}$$

where $\frac{\mathrm{\Delta }P}{L}$ is the pressure drop in the z direction, R the radius of the pipe and r the integration variable. The characteristic scales in this case are given by η_c=η₀, τ_c=ΔP and l_c=R, so that the pipe radius R, fluid viscosity η and pressure difference ΔP all take values of 1. The cubic model domain has an edge length of four units. Combining Eq. (21) with Eq. (4), the non-dimensional permeability is then given by

$$\begin{array}{}\text{(22)}& k={\displaystyle \frac{L\frac{\mathit{\pi }\mathrm{\Delta }P}{\mathrm{8}\mathit{\eta }L}{R}^{\mathrm{4}}}{\mathrm{\Delta }PA}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{128}}}.\end{array}$$

To assess the effect of different spatial resolutions, we conduct a resolution test where we increase resolution from 8³ to 256³ nodes with a constant grid spacing in each direction. Four sets of resolution tests were conducted. In the first two sets, 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 the x and y directions. 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.

Figure 3 Hagen–Poiseuille benchmark results. Shown is the error norm $\left|\right|{\mathit{\delta }}_k|{|}_{\mathrm{2}}$ vs. spatial resolution. The different curves show cases where the tube surface coincides with a nodal point (ON NODE) or not (OFF NODE). Blue lines represent simulations using stencil rescaling, whereas red lines denote simulations without stencil rescaling. To highlight convergence, black lines with given slopes were added.

As expected, the numerical results generally show higher accuracy when stencil rescaling is employed and when node 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.

4.2 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 several straight tubes with different radii (Fig. 4). We use four pipes with non-dimensional radii given as R₁=2, R₂=1, R₃=8, R₄=4. The viscosity of the fluid is 1 and edge length of the cubic domain is 8. The simulations are performed in a similar manner to the single-tube benchmark by increasing the number of grid points from 8³ to 256³. For each tube, the analytical solution (Eqs. 20, 21) is computed and the cumulative analytical permeability value is compared against computed values. The non-dimensional permeability in this case reads as

$$\begin{array}{}\text{(23)}& k={\displaystyle \frac{LQ\left(R_{\mathrm{1}}^{\mathrm{4}}+R_{\mathrm{2}}^{\mathrm{4}}+R_{\mathrm{3}}^{\mathrm{4}}+R_{\mathrm{4}}^{\mathrm{4}}\right)}{\mathrm{\Delta }PA}}.\end{array}$$

The individual tubes contribute to the absolute permeability as follows: P1=0.3662 %, P2=0.0229 %, P3=93.7514 %, P4=5.8595 %.

Similar to what was 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 a single tube, the results are more accurate in 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.

Figure 4 Multiple tube Hagen–Poiseuille benchmark. Lines and symbols correspond to the same cases as in Fig. 3. To highlight convergence, black lines with given slopes were added.

4.3 Newtonian flow through simple cubic (SC) sphere packs

In order to verify the code for more complex geometries such 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 in all directions.

The permeability of an SC sphere pack is given by Sangani and Acrivos (1982) and Bear (1988):

$$\begin{array}{}\text{(24)}& k={\displaystyle \frac{{\mathit{\varphi }}^{\mathrm{3}}\cdot d_{\mathrm{sp}}^{\mathrm{2}}}{\mathrm{180}\cdot (\mathrm{1}-\mathit{\varphi }{)}^{\mathrm{2}}}},\end{array}$$

where d_sp is the sphere diameter and ϕ is the porosity for simple cubic packing of $\mathrm{1}-\frac{\mathit{\pi }}{\mathrm{6}}\approx \mathrm{0.476}$, respectively.

Figure 5 shows the increase in accuracy with increasing number of grid points employed. The presented relative errors of the permeability value are computed in the same manner as shown in Eq. (19). The simulations employing stencil rescaling show a better convergence and seem to saturate against a relative error of 10⁻¹, demonstrating the influence of boundary effects through applied no-slip boundary conditions (finite size effect).

Figure 5 Computed $\left|\right|{\mathit{\delta }}_k|{|}_{\mathrm{2}}$ norm of the misfit between analytically and numerically computed permeabilities. The inset shows the discretization using four spheres in each direction (64 spheres in total). Streamlines are computed around those spheres and colorized with the computed velocity. Blue dots show results using stencil rescaling and red dots results with the standard method. To highlight convergence, black lines with given slopes were added.

4.4 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 behavior. For the single-tube configuration described in Sect. 4.1 and a power-law rheology, the velocity within the tube is given by, e.g., Turcotte and Schubert (2002):

$$\begin{array}{}\text{(25)}& v_z\left(r\right)={\displaystyle \frac{C_{\mathrm{1}}}{\frac{\mathrm{1}}{n}+\mathrm{1}}}\cdot {\left({\displaystyle \frac{\mathrm{\Delta }P}{L}}\right)}^{\frac{\mathrm{1}}{n}}\cdot \left({\left[{\displaystyle \frac{R}{\mathrm{2}}}\right]}^{\frac{\mathrm{1}}{n}+\mathrm{1}}-r^{\frac{\mathrm{1}}{n}+\mathrm{1}}\right),\end{array}$$

where $C_{\mathrm{1}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{0}}^{-\frac{\mathrm{1}}{n}}$ (see Appendix B), R is the tube radius and r the width of the tube in Cartesian coordinates. Figure 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).

Figure 6 Comparison of analytical and numerical velocities for Hagen–Poiseuille flow with a power-law fluid. Analytical velocities are represented as colored lines and numerical velocities as colored symbols.

5 Application to Fontainebleau sandstone

To verify the ability of the code to handle more complex flows through natural samples and to validate previously computed permeability values, we used the CT data for a Fontainebleau sandstone sample provided by Andrä et al. (2013b) with dimensions 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³ 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 permeability value of ≈1100 mD (Keehm, 2003).

5.1 Newtonian flow

As in previous tests, we compute the permeability of the extracted subsample using Eqs. (4) and (5). Figure 7a shows streamlines colored using computed fluid velocities and Fig. 7b the local pressure.

Figure 7 Newtonian fluid flow through the Fontainebleau sandstone sample. Streamlines colored using computed fluid velocities are shown in panel (a) and streamlines colored using fluid pressures are shown in panel (b).

For a resolution of 256³, we obtain permeabilities which are comparable to previously computed permeabilities of the same sample (Fig. 8; Andrä et al., 2013b), with the rescaled stencil 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³, 768³, 1024³). 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. Figure 8 shows 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³, the computed permeability values decreased by ≈24.6 % when using a grid resolution of 1024³. Additionally, the benefit of stencil rescaling can also be seen here, as, e.g., the simulation with a resolution of 512³ 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). Figure 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 behavior has also been observed in Lattice–Boltzmann method (LBM) simulations (e.g., Khirevich et al., 2015; Khirevich and Patzek, 2018).

Figure 8 Computed permeability values against grid resolution. Orange symbols denote simulations using Lattice–Boltzmann method (LBM) and explicit jump Stokes (EJ Stokes); both methods are used in Andrä et al. (2013b). Blue data points represent simulations using stencil rescaling, while simulations represented by red dots use the standard method. The brown dotted line symbolizes the experimental estimate from Keehm (2003).

5.2 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 power-law rheology given in Eq. (7). The edge length of the subsample was 1.92 mm, which corresponds to a CT resolution of 256³ voxels. To increase accuracy, we increased this resolution by a factor of 2 to a resolution 512³. As seen in Sect. 5.1, 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 η₀=1 Pa and η₁ and η₂ were set to 10⁻³ and 10⁶, 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 to 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 $(\mathrm{19}\mathrm{982}\pm \mathrm{9})\times {\mathrm{10}}^{-\mathrm{4}}$ and $(\mathrm{1}\mathrm{000}\mathrm{009}\pm \mathrm{3})\times {\mathrm{10}}^{-\mathrm{6}}$ are in good agreement with the imposed power-law coefficients of 0.5 and 1 (Eq. 8).

Figure 9 Computed results on the Fontainebleau sample using non-Newtonian rheology. Panel (a) shows the mean velocity against the applied pressure at the top boundary. Red and blue triangles symbolize each simulation, and the corresponding dotted black line represents the fitted curve through the obtained data with slope $\frac{\mathrm{1}}{n}$. Panel (b) illustrates computed streamlines of the Fontainebleau subsample using a power-law coefficient of 0.5. Solid material is displayed in grey and the streamlines are colored according to computed velocities.

6 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 the case that the fluid–solid interface and the numerical grid are at least partially aligned. Cases using the stencil rescaling solutions with a velocity change on a computational node produce smaller relative errors.

Similar to studies using the 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 the Fontainebleau sandstone sample with grid resolution of 1024³, calculations employing stencil rescaling give approximately the same permeability value as suggested by laboratory measurements, while simulations without employing stencil rescaling overestimate the computed permeability by ≈14.72 %. (Fig. 8). However, this behavior may also be the opposite depending on the numerical implementation of the respective numerical method (e.g., Khirevich et al., 2015; Khirevich and Patzek, 2018; Zakirov and Galeev, 2019).

The computation of permeabilities in a three-dimensional pore space using micro-CT data strongly depends on the reasonable quality of the micro-CT images followed by several steps of segmentation in order to resolve tiny fluid pathways. Although high-quality input data are 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 have 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 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).

Furthermore, it should be kept in mind that solver options like convergence criteria may influence the obtained permeability result. Figure 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 numbers of CPUs depending on the size of the computed domain. As an example, a simulation with 512³ voxels uses 1024 CPUs, 185 GB RAM and requires 5790 s computation time. Apart from LaMEM, finite difference codes like FDMSS (Gerke et al., 2018) 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.

Appendix A: Convergence criteria

To determine whether a numerical solution converges, two convergence criteria, which are absolute and relative convergence tolerance, are used. 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 $<{\mathrm{10}}^{-\mathrm{7}}$. Thus, for all further simulations, a relative convergence tolerance of 10⁻⁸ is used (Fig. A1). A change in the absolute convergence tolerance did not have any effect on the computed solution; therefore, we use an absolute convergence tolerance of 10⁻¹⁰.

Figure A1 Results of simulations for (a) Hagen–Poiseuille single tube, (b) simple cubic sphere pack and (c, d) Fontainebleau sandstone using different relative/absolute convergence tolerances.

Appendix B: Definition of C₁

C₁ is an constant arising during the derivation of Eq. (25), which is related to the non-linear rheology used in Turcotte and Schubert (2002). This rheology is written as

$$\begin{array}{}\text{(B1)}& \dot{\mathit{ϵ}}=C_{\mathrm{1}}{\mathit{\tau }}^{n^{\prime }},\end{array}$$

where n^′ is the stress exponent as used for power-law materials in geodynamics. Replacing τ with $\mathit{\tau }=\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}}$ leads to

$$\begin{array}{}\text{(B2)}& \dot{\mathit{ϵ}}=C_{\mathrm{1}}(\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}}{)}^{n^{\prime }}.\end{array}$$

Solving Eq. (B2) for η results in

$$\begin{array}{}\text{(B3)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{1}}^{-\frac{\mathrm{1}}{n^{\prime }}}{\dot{\mathit{ϵ}}}^{\frac{\mathrm{1}}{n^{\prime }}-\mathrm{1}}.\end{array}$$

We can now define a reference viscosity η₀ at a reference strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}$. This reference viscosity then reads as

$$\begin{array}{}\text{(B4)}& {\mathit{\eta }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{1}}^{-\frac{\mathrm{1}}{n^{\prime }}}{\dot{{\mathit{ϵ}}_{\mathrm{0}}}}^{\frac{\mathrm{1}}{n^{\prime }}-\mathrm{1}}.\end{array}$$

Assuming $\dot{{\mathit{ϵ}}_{\mathrm{0}}}=\mathrm{1}$ and solving for C₁ then provides us with the following expression:

$$\begin{array}{}\text{(B5)}& C_{\mathrm{1}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{0}}^{-n^{\prime }}=\mathrm{2}{\mathit{\eta }}_{\mathrm{0}}^{-\frac{\mathrm{1}}{n}}.\end{array}$$

Appendix C: Permeabilities of different Fontainebleau subsamples

In order to show numerical convergence of the given Fontainebleau sample, several subsamples were extracted and the resolution increased to 512³, 768³ and 1024³ grid points. Figure C1 displays the convergence with increasing grid resolution. The different subsamples show a variance of around 12 % for the computed permeability value.

Figure C1 Numerical convergence of different Fontainebleau subsamples with increasing grid resolution. All subsamples displayed were computed using stencil rescaling. For comparison, the computed permeabilities from Andrä et al. (2013b) are shown. The dotted brown line symbolizes the experimental estimate taken from Keehm (2003).