Permeability Estimation of Regular Porous Structures: A Benchmark for Comparison of Methods

The intrinsic permeability is a crucial parameter to characterise and quantify fluid flow through porous media. However, this parameter is typically uncertain, even if the geometry of the pore structure is available. In this paper, we perform a comparative study of experimental, semi-analytical and numerical methods to calculate the permeability of a regular porous structure. In particular, we use the Kozeny–Carman relation, different homogenisation approaches (3D, 2D, very thin porous media and pseudo 2D/3D), pore-scale simulations (lattice Boltzmann method, Smoothed Particle Hydrodynamics and finite-element method) and pore-scale experiments (microfluidics). A conceptual design of a periodic porous structure with regularly positioned solid cylinders is set up as a benchmark problem and treated with all considered methods. The results are discussed with regard to the individual strengths and limitations of the used methods. The applicable homogenisation approaches as well as all considered pore-scale models prove their ability to predict the permeability of the benchmark problem. The underestimation obtained by the microfluidic experiments is analysed in detail using the lattice Boltzmann method, which makes it possible to quantify the influence of experimental setup restrictions.

Pressure interaction force acting from particle j on particle i, [N] F V i j Viscous interaction force acting from particle j on particle i, [N]

Introduction
The characteristics of flow through porous media play an important role in a wide range of natural and industrial applications. A classical example is found for soil through which groundwater seeps. In this context, we can intuitively understand the terms porosity and permeability. Porosity is a measure of the void space in the porous medium, whereas permeability is a measure for the resistance of the porous medium itself against flow. In this regard, the famous law of Darcy (1856) provides a correlation between flow and the corresponding driving force in a porous medium via Here, v is the fluid's filter (Darcy) velocity, p is the pressure, μ is the dynamic viscosity and K is the second-order intrinsic permeability tensor. In particular, the Darcy velocity v = φ w F is the fluid's seepage velocity w F weighted by the porosity φ. Gravitational forces are neglected in (1). The computation of the total (void) porosity for porous media with a known geometrical structure is straightforward. However, this is not the case for the permeability. The permeability is crucial for physically consistent modelling and accurate numerical simulations of flow and transport processes in porous media, but associated with great uncertainty. Therefore, it deserves our particular attention.
In this paper, we use the following definition of scales to distinguish the relevant processes and parameters accordingly. The intrinsic permeability is an effective parameter that accounts on the scale of a representative elementary volume (REV scale) for the geometric configuration on the pore scale. The REV scale is the one where Darcy's law is valid; it is therefore often denoted also as the Darcy scale. Our proposed methodology involves different mathematical and numerical models that all resolve the flow through the porous medium on the pore scale. It is compared to an experimental determination of permeability, which is inherently at the REV scale. For further details on permeability and porosity as effective REV-scale parameters, we refer to e. g. Helmig (1997); Hommel et al. (2018).
Different experimental, semi-analytical and numerical techniques to compute or estimate the permeability exist in the literature. The intrinsic permeability can be determined experimentally by imposing a flux, measuring the corresponding pressure drop and applying Darcy's law. In recent years, microfluidic experiments have been increasingly used for various experimental investigations of porous media, cf. Yoon et al. (2012); Karadimitriou and Hassanizadeh (2012); Scholz et al. (2012).
Several approximations for the semi-analytical description of the permeability-porosity relationship exist, the most famous being the Kozeny-Carman equation, cf. Kozeny (1927); Carman (1997). In this context we mention also the Hazen relation (Hazen 1892), the powerlaw relation by Verma and Pruess (Verma and Pruess 1988) and the Timur and Morris-Biggs relations (Timur 1968;Morris et al. 1967). For a recent overview of porosity-permeability relationships for evolving porous media we refer the reader to Hommel et al. (2018). These relations have in common that they try to describe all variability using the porosity and one or several fixed scaling factors. In practice, these scaling factors are often fitted to the observed permeability data, if available. However, the structure and the permeability of a porous medium is not solely described by its porosity; the same porosity can yield different permeability values, e.g. Schulz et al. (2019). The shape, distribution, interfacial tension and roughness of the grains as well as the coordination number, connectivity and size of the pores significantly affect the permeability of the medium (Berg 2014;Millington and Quirk 1961;Valdes-Parada et al. 2009;. However, this information is unfortunately unknown for most settings. The Kozeny-Carman relation includes geometrical information through the tortuosity, which can be estimated for simple grain packings (Yazdchi et al. 2011). An extensive review concerning the validity of the Kozeny-Carman relation and its modifications for different porous-medium geometries is provided in Schulz et al. (2019).
Different averaging techniques can be applied to compute the permeability: the homogenisation theory based on two-scale asymptotic expansions, the volume averaging theory and the numerical upscaling (Whitaker 1999;Hornung 1997;Auriault et al. 2009;Gray and Miller 2014). The theory of homogenisation provides a useful tool to efficiently compute the permeability of regular porous media requiring solutions to flow problems in a periodic unit cell. Another advantage of the method is that it is not restricted to the scalar case. Depending on the geometrical characteristics (e.g. thin porous media, parabolic velocity profile) different assumptions can be made and, therefore, different flow problems in unit cells are solved (Fabricius et al. 2016;Chamsri and Bennethum 2015). Volume averaging (Whitaker 1999;Gray and Miller 2014) offers an alternative technique to compute permeability based on solving local flow problems. However, under periodic closure conditions, which is the case for the benchmark problem considered in this work, the resulting set of equations coincides with the one obtained from homogenisation (Valdes-Parada et al. 2009). Therefore, volume averaging is not considered in this paper.
Finally, solutions of the pore-scale resolved models can be obtained using the lattice Boltzmann method (LBM), Smoothed Particle Hydrodynamics (SPH) or the Finite Element Method (FEM) and the simulation results can then be upscaled in order to compute the intrinsic permeability (Pan et al. 2001;Sivanesapillai et al. 2014).
In reservoir simulations with a length scale of kilometres, a field scale is typically introduced. In this case, Darcy's law is considered to be valid both on the Darcy scale and the field scale. The Darcy-scale permeability tensor is usually highly heterogeneous and needs to be upscaled in order to develop efficient numerical methods (Farmer 2002;Wu et al. 2002;. However, such methods are beyond the scope of this article. In the literature, comparisons of individual methods for permeability estimation exist (Schulz et al. 2019;Chamsri and Bennethum 2015;Song et al. 2019;Sugita et al. 2012;Guibert et al. 2015;Yazdchi et al. 2011). For highly complex geometries, the numerical methods can generally only be applied to a small subsample of the domain, whereas experimental and empirical relations are applied for the domain as whole, giving inherent uncertainties and inaccuracies due to rock heterogeneities (Song et al. 2019;Guibert et al. 2015). Therefore, we perform a comparative study of a broad variety of methods with a designed benchmark experiment. The goal of this paper is to estimate the intrinsic permeability of an artificially produced regular porous medium using different techniques, to investigate the applicability of each method and to validate the considered approaches for computing effective properties of porous media. The manuscript can serve as a benchmark for researchers working on modelling flow and transport processes in porous media, where the geometric structure is available but the effective parameters are not known (artificially produced composite materials, geometry reconstructed by imaging techniques). The results of the paper will also be of use for scientists working on the development of advanced averaging methods for computing effective properties of porous media.
We would like to draw the reader's attention to the limitations of this work in terms of scales. Since one of the objectives of this work is to compare the numerical schemes with physical experiments, the experimental potential was the limiting factor in terms of the length scales under consideration. Our experimental infrastructure would not allow for formations with ultra small features (smaller than a few microns). Additionally, even if such features were possible to achieve, in terms of the fabrication method for the artificial porous medium employed, this would induce pressures which would lead to the deformation of the material, leading to increased conditional inaccuracies. Consequently, this experimental approach in general cannot account for very low permeabilities (nano-or micro-Darcy), which would also demand for the use of gases instead of fluids, and a fundamentally different approach to estimate the corresponding permeability (Klinkenberg effect). In accordance to this limitation, the corresponding numerical schemes under investigation and comparison, also do not bear these characteristics and assumptions.
The paper is organised as follows. In Sect. 2, the setup of the benchmark problem, which serves as the basis for all considered methods, is described. Section 3 presents calculations using the Kozeny-Carman equation providing a semi-analytical porosity-permeability relation. Further, four mathematical homogenisation approaches are given in Sect. 4. In particular, a three-dimensional (3D) approach is described in Sect. 4.1, a classical two-dimensional (2D) approach in Sect. 4.2, a very thin porous medium (VTPM) approach in Sect. 4.3 and a pseudo 2D/3D approach in Sect. 4.4. Numerical methods to compute the intrinsic permeability of the benchmark problem via upscaling the pore-scale simulations are discussed in Sect. 5, namely FEM in Sect. 5.1, SPH in Sect. 5.2 and LBM in Sect. 5.3. The experimental setup and the experimental results of the benchmark problem are presented in Sect. 6. Finally, a comparison of the results is given in Sect. 7 and a concluding discussion is provided in Sect. 8.

Benchmark Problem
As a common basis for the discussion of different permeability estimates, it is crucial to precisely describe the chosen geometry of the porous structure and the related modelling assumptions. As a benchmark example, a simple regular and homogeneous porous medium is chosen. In a 3D domain of dimensions 14 mm × 10 mm × 0.091 mm (length l × width b × thickness h), equidistantly aligned cylinders with the same radius r are embedded, cf. Fig. 1. The radius for the manufactured sample investigated in the experiments is r = 0.4 mm, cf. Sect. 6. In addition, the other methods consider also radii of 0.35 mm, 0.45 mm, 0.47 mm and 0.49 mm. The smallest repeating unit consists of a square of 1 mm edge length. We chose the REV having one cylinder in the centre, which corresponds to the unit cell in this example. The computation of the porosity φ is straightforward and given in Table 1 for different radii. In terms of basic modelling assumptions, we consider no-slip boundary conditions at the internal boundaries (interfaces) of the cylinders as well as at the top and bottom surface of the domain. Furthermore, we assume creeping flow conditions of a viscous liquid at very low Reynolds numbers and, thus, neglect inertia effects. Finally, we assume the solid structure as impermeable and rigid, i. e. no solid deformations.

Semi-Analytical Approach
The well-known Kozeny-Carman equation is a semi-analytical, semi-empirical relation for estimating the permeability of porous media (Kozeny 1927;Carman 1997): where k is the intrinsic permeability, φ is the porosity, σ is the specific surface area and c kc is the Kozeny-Carman constant. The porosities φ of the considered geometries are given in Table 1. The specific surface area for cylinders is σ = 2/r . In general, it is difficult to estimate the value of c kc since it depends on many factors including the flow tortuosity and roughness of the grains. Carman (1997) originally proposed the value c kc = 5 or granular porous media, which has also been suggested for cylindrical grains by Schulz et al. (2019). We will hence also use c kc = 5 in the current study. Table 3 lists the values of permeability k obtained from relation (2) for different radii of solid cylinders. A typical procedure in practical applications is to use the Kozeny-Carman constant as a fitting parameter. This is not done here as our goal is to find the permeability values given by the Kozeny-Carman equation without using extra information to fit parameters. However, fitting of c kc for cylindrical obstacles does not give one unique value of c kc applicable for several cylinder radii (Yazdchi et al. 2011). Moreover, the Kozeny-Carman equation (2) provides a permeability estimation for porous media, which are extended over infinite space. Thus, the influence of the walls in the benchmark problem (cf. Fig. 1) is not considered in the Kozeny-Carman approximation, cf. the discussion in Sect. 7. For modifications of the Kozeny-Carman relationship (2) we refer the reader to e.g. Schulz et al. (2019).

Mathematical Homogenisation
In this section, we compute the permeability using the theory of homogenisation (Auriault et al. 2009;Hornung 1997). We consider fluid flow through the three-dimensional porous with a periodic arrangement of pores, cf. Fig. 1a. We denote by ε = δ/ the ratio of the characteristic pore size δ to the length of the domain of interest. At the pore scale, flow in the pore space ε is described by the steady Stokes equations completed with the no-slip condition on the boundary of solid inclusions and appropriate boundary conditions on the external boundary ∂ . We denote by v ε and p ε the non-dimensional velocity and pressure of the fluid. This problem formulation ensures non-trivial flow in the range of Darcy's law (Hornung 1997).
As it is common in the theory of homogenisation, we define the porous microstructure by periodic repetition of the scaled unit cell Y , which consists of the solid part Y s and the fluid part Y f , cf. Fig. 2. To obtain the permeability of the porous medium we follow the classical procedure of homogenisation (Hornung 1997). We assume two-scale asymptotic expansions for the pore-scale velocity and pressure where a, b ∈ N 0 depend on the homogenisation approach applied and v i , p i are y-periodic functions. Computing the derivatives ∇ = ∇ x + ε −1 ∇ y , substituting expansions (4) into the pore-scale problem (3) and combining terms with the same degree of ε, we obtain Darcy's law (1) valid in the porous-medium domain .
We consider different assumptions for the geometric structure of the medium (full 3D, 2D, very thin porous medium and pseudo 2D/3D) with the appropriate unit cells (circular and cylindrical solid inclusions) and apply four different homogenisation approaches to compute the intrinsic permeability. We use the software package FreeFem++ (Hecht 2012) for numerical simulations of all homogenisation approaches considered in this section.

Classical Three-Dimensional Approach
We consider solid cylinders with height h > 0. Thus the corresponding unit cell Y phys = Y phys s ∪ Y phys f is also characterised by the macroscopic height h, cf. Fig. 2b. To compute the permeability by means of homogenisation we solve the following threedimensional cell problems for j = 1, 2, 3: taking into account Y phys f π j dy = 0. The non-dimensional intrinsic permeability is given by and π j are the solutions to the cell problems (5) for j = 1, 2, 3. For the considered geometry (cylindrical solid inclusions) we obtain the diagonal permeability The intrinsic permeability k is computed numerically, the units are taken into account and the physical permeability values are presented in Table 3 for different cylinder radii. Taylor-Hood P2/P1 elements are applied for the velocity and the pressure, respectively. For the solid inclusions with radius r = 0.4 the 3D unit cell Y phys is partitioned into twelve elements in x 3 -direction and approx. 29 000 elements in the whole fluid part of the unit cell Y phys f .

Classical Two-Dimensional Approach
The classical two-dimensional approach corresponds to the case where no effect from the top and bottom of the porous domain is included, and hence corresponds to the height h of the medium being infinite. This is a very strong simplification which is not valid for our benchmark in general (see discussion in Sect. 7). However, in comparison to the 3D homogenisation approach the 2D approach is computationally much cheaper. In this case, we define a two-dimensional unit cell Y (Fig. 2a) and obtain the cell problems in two space dimensions ( j = 1, 2): where Y f π j dy = 0. With this approach we obtain the non-dimensional permeability and π j are the solutions to the cell problems (6) for j = 1, 2. Again, we get a diagonal tensor K = (k i j ) i, j=1,2 = diag{k, k} and present the physical permeability values k in Table 3 for different radii. For discretisation of velocity and pressure the Taylor-Hood finite element pair is used with approx. 20 500 elements for the radius r = 0.4.
We observe that for highly porous structures (r = 0.35, r = 0.4) the permeability is one order of magnitude higher than the one obtained by the 3D approach. This is due to the fact that the 2D approach neglects the wall friction at the top and bottom, compared to the 3D approach, where the no-slip condition is valid at these boundaries.

Very Thin Porous Media
We apply the homogenisation approach proposed in Fabricius et al. (2016) for Very Thin Porous Media (VTPM) where the cylinder height h is much smaller than the distance δ − 2r between the cylinders. In this case, the relation h = δ 2 is used for the derivation of the permeability tensor. The main idea is to pass the limit h/ε → 0 in the 3D approach. In this case we get the following cell problems for j = 1, 2: taking into account Y f w j dy = 0. The non-dimensional permeability is given by where w j is the solution to the cell problems (7) for j = 1, 2. As in Sects. 4.1 and 4.2 , we obtain a diagonal permeability tensor K = diag{k, k} and present the physical permeability values k in Table 3. For the numerical simulations P1 finite elements are used, and for inclusions of the radius r = 0.4 the fluid part of the unit cell is partitioned into approx. 20 500 elements. The cell problems (7) are computationally very cheap compared to problems (5) and (6). However, applying the VTPM approach we assume that the height h is much smaller than the interspatial distance between two solid obstacles (Fabricius et al. 2016) and the fluid flow in x 3 -direction is negligible. Therefore, this approach fails for porous media where the distance between the cylinders is smaller than the height. For r = 0.45 the interspatial distance is 0.1 > h. The interspatial distance for r = 0.47 and r = 0.49 is 0.06 < h and 0.02 < h, accordingly. Therefore, the computational results for the two last radii in Table 3 (homogenisation VTPM) are not relevant.

Pseudo Two-/Three-Dimensional Approach
Another simplification of the 3D homogenisation approach is obtained assuming that the velocity is horizontal and its vertical variability is resolved through a parabola as no-slip boundary conditions are assumed and the viscous boundary layer is larger than half of the cell height. Such a simplification is applicable when the height of the domain is dominating the flow profile. The procedure is in accordance to Flekkøy et al. (1995), where the Stokes flow between two parallel plates is considered. In particular, the effect from the top and bottom boundaries is incorporated into the model through a viscous drag force. Hence, this pseudo 2D/3D homogenisation approach also relies on the porous medium being thin, where the thickness is still comparable to the typical pore size. Therefore, only the flow in the middle of the domain (x 3 = h/2) needs to be considered, and the following 2D-cell problems are obtained ( j = 1, 2): The non-dimensional permeability in this case is given by where the scaling factor 2/3 comes from the integral of the parabolic profile. As in Sects. 4.1-4.3, we obtain a diagonal permeability tensor K = diag{k, k} and present the physical permeability values k in Table 3. For the numerical simulations we applied Taylor-Hood elements and for the setting with radius r = 0.4 mm approx. 20 500 elements are used.
Computational costs for the cell problems (8) are the same as for the classical 2D approach but the permeability estimation complies with the 3D approach for moderate cylinder radii (r < 0.47 mm). However, the difference between the 3D approach and the pseudo 2D/3D approach increases for larger radii (r = 0.47 mm and r = 0.49 mm). This is due to decreasing gaps between the cylinders leading to a violation of the assumption on the parabolic flow profile.

Pore-Scale Models
In this section, three different numerical methods for pore-scale resolved computations are applied to the benchmark problem, namely the FEM, SPH and LBM. In these models, implementations of the Navier-Stokes equations are available and, thus, used in the numerical computations. However, the solution of this more general formulation reduces to the solution of the stationary Stokes problem in the limit case of Reynolds numbers tending towards zero, as it is postulated in the benchmark problem.

Finite-Element Method for the Navier-Stokes Equations
The simulation is realised using the commercial finite-element (FE) software tool Abaqus /CFD 1 . Thereof, a FE implementation of the incompressible Navier-Stokes equations is used for the performed computational fluid dynamics (CFD) simulations. For the flow domain, a constant flow rate of 10 μL/min is prescribed at the inlet surface, while the outlet surface is assigned with a constant pressure of 0 Pa, cf. Fig. 3a.
The remaining surfaces, i. e. the top, bottom, lateral and internal surfaces of the cylinders, are assigned with no-slip boundary conditions. For the spatial discretisation, a finite-element mesh is generated with linear hexahedral elements since quadratic elements are not available in Abaqus for fluids. Linear elements are not optimal to approximate the parabolic velocity profile but numerically much cheaper. Therefore, the geometry is meshed in total with 327 440 elements using ten elements along the thickness of the structure in order to adequately resolve the arising parabolic velocity profile. For the material properties of the liquid, a density of 1 000 kg/m 3 and a viscosity of 8.9 × 10 −4 Pa s are chosen according to water at 25 • C. For the considered example, a linear pressure drop from the left to the right side (i. e. a constant pressure gradient in the steady state) is obtained, cf. Fig. 3b. The volumetric flow rate through the pore structure is provided as an output variable within Abaqus. Thus, the Darcy filter velocity can be computed by dividing it with the permeable cross-sectional area. Based on these quantities, the intrinsic permeability is obtained via the Darcy equation (1), applied in the flow direction. The estimated permeability values are collected in Table 3.

Smoothed Particle Hydrodynamics for the Navier-Stokes Equations
In Smoothed Particle Hydrodynamics (SPH), the discretisation of the governing weakly compressible Navier-Stokes equations spans a set of integration points x i , also referred to as particles, cf. Monaghan (1988); Gingold and Monaghan (1977). At each of these interacting collocation points, scalar fields , cf. Nomenclature of physical properties are interpolated employing convolution with a continuously differentiable kernel function W (x, R). The kernel radius R determines a sphere of influence and likewise declares neighbouring particles  j. The approximation of integrals converts continuous field functions into particle properties cf. Nomenclature, kernel representation into spatial discretisation W (x, R) = W i j , and turns differential operators into short-range interaction forces. Hence, the motion equation to be solved for each fluid particle is found as where F V i j denotes viscous interaction forces, F P i j pressure interaction forces and F B i body forces. In this context, the implemented SPH formulation is derived on top of the tool HOOMD-blue (Anderson et al. 2008;Glaser et al. 2015;Sivanesapillai et al. 2016) providing state of the art no-slip no-penetration fluid-solid boundaries (Adami et al. 2012), ghost-particle methods for periodic boundaries (Sivanesapillai et al. 2016), a Verlet time integration scheme and stabilisation via an artificial viscosity term (Monaghan and Gingold 1983). Internal surfaces and boundaries in x 3 -direction are assigned with no-slip conditions. The domain is periodic in x 1 -and x 2 -directions, cf. Fig. 1a. For the liquid, material parameters of water are chosen. A body force equivalent to a pressure gradient is applied as the driving force.
In terms of the considered benchmark problem, cf. Fig. 1, the arising flow is computed for a single REV, cf. Fig. 1b, as in homogenisation. The REV structure is discretised using a resolution of 200×24×200 particles for all radii configurations. The discrete structure of the REV and the velocity field in steady state are shown in Fig. 4. Based on the resulting velocity profiles of the SPH simulation, the permeability is computed using the rearranged version (9) of Darcy's law (1). As in Sect. 5.1, the mean particle velocity w F (seepage velocity) needs to be multiplied by the porosity φ to obtain the filter velocity v used in (1). The computed permeabilities are collected in Table 3.
As a basic requisite of the SPH approach the fluid domain is discretised with fluid particles. For a meaningful approximation of the velocity profile in the pore-space passages, at least 10-15 particles are required. For example, with the chosen resolution, a radius r = 0.35 mm results in 60 particles, while r = 0.47 mm allows only 12 particles over the passage. This leads to an enormous increase of the computational effort for larger radii (r > 0.47 mm) due to the need to increase the resolution. Therefore, the SPH is not suited to simulate the benchmark problem for a radius of 0.49 mm with acceptable computational costs.

Lattice Boltzmann Method
The lattice Boltzmann method (LBM) is based on the mesoscopic representation of fictional particle movements (Succi et al. 1991) and provides the third pore-scale resolved method discussed in this paper. Particles stream on a lattice with discrete velocities and collide to relax to an equilibrium. The probabilities f i (x, t) of finding a particle with velocity c i at lattice node x at time t evolve according to the lattice Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) collision operator BGK , the time step t and the number Q of discrete velocities. The collision operator governs the rate at which the probabilities f i relax towards the equilibrium distribution where w i are the weights of each discrete velocity and c s is the speed of sound. The discretisation of velocities limits the applicability of the lattice Boltzmann method to flows with a maximum Mach (Ma) number of 0.3 (Krüger et al. 2017), which is the case here (Ma ≈ 0.001). It is commonly known that the continuum flow approximation holds up to a Knudsen number of Kn ≈ 0.001. For the application and flow rates chosen here Kn 0.001, taking the cylinder diameter (millimetre range) as reference length. Further details on the Knudsen number limits of the LBM can be found in Silva and Semiao (2017); Silva (2018), whereas He and Luo (1997) provides a detailed derivation of lattice Boltzmann equation from the Boltzmann equation. The LBM is not directly applicable for high-speed flows (Re > 10000) as the spatial and temporal resolutions required for such flows become computationally prohibitive (Jain et al. 2016). For the benchmark problem, Re < 1.0 (Fig. 8b). Therefore, the lattice Boltzmann equation is valid within the limits for the application considered here. The macroscopic velocity v and pressure p are obtained from the particle distribution functions, cf. Succi et al. (1991). For this benchmark problem, a D3Q19 stencil is chosen, i. e. 19 discrete velocities (Q = 19) in three dimensions (D = 3) are used (Krüger et al. 2017). No-slip boundary conditions are prescribed at the obstacles as well as the top and bottom walls of the domain. A flow rate guaranteeing a low Reynolds number is prescribed in the form of velocity boundary conditions at the inlet and outlet. The simulations are performed using the LBM solver of the open-source software ESPResSo, cf. Arnold et al. (2013); Roehm and Arnold (2012). In order to justify the resolution for the LB simulation, a grid convergence study (using a quarter of the benchmark setup) is performed prior to the calculation of the permeability of the benchmark geometry, cf. Fig. 5.
In accordance to Sect. 5.1, at least 10 lattice cells are required to resolve the quadratic velocity profile in the pore space. Therefore, a height resolution of 10 lattice cells is chosen to minimise the numerical effort. For the radius r = 0.49 mm, the lattice resolution was increased to 40 cells in x 3 -direction and only a subset of 4×1 cylinders, periodically repeated in x 2 -direction, is considered. The permeability is estimated using the computed velocity and pressure values according to (9). The values are collected in Table 3.

Micromodel Experimental Setup
The experimental setup is comprised of the micromodel, a syringe pump, the pressure sensors and a computer for data acquisition. The micromodel is manufactured out of Poly-Di-Methy-Siloxane (PDMS), following the principles of soft-lithography (Xia and Whitesides 1998;Karadimitriou et al. 2013;Karadimitriou and Hassanizadeh 2012;Auset and Keller 2004). Using this technique, micromodels with a very low surface roughness, in the order of a few tens of nanometers, are produced. Therefore, the influence of the surface roughness on the permeability estimation of the benchmark problem is negligible. The micromodel's flow network is presented in Fig. 6. A CETONI neMESYS 1000N mid-pressure syringe pump 2 is used for the introduction of the fluid, in combination with 2.5 ml glass syringes. The flowinduced pressure difference between the inlet and the outlet of the flow network p, cf. Fig. 6, is measured with two Elveflow MPS2 pressure sensors 3 , with a measurement range of 0-1 bar and acquired with a 16 bit data acquisition system (CETONI QMix I/O). The flow network is connected to the syringe pump and the pressure sensors via Teflon (Polytetrafluoroethylene, PTFE), 1/16 OD, 0.5 mm ID tubing. The connected computer is able to control the syringe pump and communicate with the pressure sensors via a CETONI BASE 120 module 2 , cf. Fig. 7.

Permeability Estimation
The intrinsic permeability of the flow network is estimated by imposing a flux and measuring the arising pressure drop p = p 2 − p 1 , cf. Fig. 6. For this purpose Darcy's law (1) is employed. In order to increase the accuracy of the results and get a better estimation of the measurement error, three runs of series of different flow ratesxxx ranging from 1 to 7 L/s are applied for one minute each, cf. Fig. 8, while the induced pressure is measured at the inlet and outlet, cf. Fig. 6.
The absolute pressure drop across the length of the flow network is directly proportional to the flow rate. Therefore, a linear regression was applied. Rearranging Darcy's law (1), we obtain the following expression for the intrinsic permeability where A = bh is the cross-sectional area of the porous domain, l is the length of the porous domain and μ = 10 −3 Pa s is the dynamic viscosity of water. The permeability is calculated by equation (13) using the slope s pq (pressure over flow rate) of the linear regression, cf. Table 2. 1.4 × 10 −2 9.1 × 10 −7 3.9 × 10 +11 3.9 ± 0.2 × 10 −11 The corresponding Reynolds numbers range from 0.1 to 0.7, cf. Fig. 8b. The error bars in Fig. 8b are related to the pressure measurement. Note that apart from the quantifiable error in the pressure measurement, there are additional uncertainties in the procedure towards obtaining the permeability of the porous domain. The pressure drop caused by the triangular inlet and outlet domains is neglected in this case, cf. Fig. 6. Also, during the photo-lithography step of the manufacturing process, the photo-resist-covered silicon wafer is exposed to ultra violet light under a mask. If the illumination source is not collimated, this will create a slope in the resulting photo-resist walls, which will later be transferred to the walls of the actual micromodel via the soft lithography process. This bias in the process will eventually create channels in the flow network with a smaller cross-sectional area than the intended one, leading to a potentially smaller intrinsic permeability of the network in total. As expected, this is found in comparison to the homogenisation and pore-scale methods, cf. Fig. 11. Therefore, further studies of the full microfluidic setup (Fig. 6) are performed using LBM.

Numerical Study of the Microfluidic Setup
For the simulation of the entire microfluidic device with the inlet and outlet, the computer aided design (CAD) file was discretized with about 176 million lattice cells. This resolution was chosen based on our previous comprehensive mesh convergence study (Jain et al. 2016). For this simulation the Musubi LBM solver (Klimach et al. 2014) is chosen since this feature is currently not implemented in the ESPResSo framework. Simulations are executed using 512 CPU cores of the compute cluster installed at the Institute for Computational Physics at the University of Stuttgart. The boundary conditions are prescribed according to the experiments with a flow rate of 1 L/s, which translates to a velocity boundary condition at the inlet and a zero-pressure condition at the outlet.
The pressure drop across the full device ( p 2 − p 1 ) was calculated as p LBM = 3.96 mbar from LBM simulations, which is slightly different from the experimental measurement p exp = 4.2 ± 1.2 mbar, but well within the error margin. This indicates that the inclusion of the inlet and outlet in the numerical simulations results in a better agreement with the experimental values, compared to the LBM simulations without the inlet and outlet (Fig. 5).
To further quantify the role of the inlet and outlet, we analysed the flow physics in these regions, and its effect on the flow field within the porous structure thereof. As can be observed from Fig. 9c, the inlet geometry leads to an uneven flow distribution not only in the inlet itself, but also ranging into the array of cylinders. These observations are reflected in Fig. 10a which shows the normalised pressure across the micromodel geometry. The pressure values are obtained by averaging across various x 2 − x 3 cross sectional planes placed between every column of cylinders in the porous domain, and at a distance of 1mm in the inlet and outlet domains. The dotted vertical lines indicate the beginning and end of the porous domain corresponding to M1 and M2 in Fig. 9a.
It can be clearly seen from Fig. 10a that the disturbed flow from the tortuous inlet causes significant pressure drop even after the first array of cylinders, and the outlet has same effect up to second-to-last array of cylinders (highlighted with grey vertical bars). To further elucidate the role that these irregularities play in the permeability estimates we calculate the stepwise pressure drop p i at the cross sectional plane i as a central derivative, i. e. p i = ( p i−1 − p i+1 )/2. The results are shown in Fig. 10b and confirm that the effects of the inlet and outlet channels carry over into the regular domain about one column deep before diminishing. Quantitatively, the permeability computed from the pressure drop between M1 and M2 is 17.4 × 10 −11 m 2 , while the permeability computed from the linear pressure drop region between the second and second to last column of cylinders is 15.9 × 10 −11 m 2 . This suggests that even if in the experiment the pressure is measured directly in front of the regular domain, there would be still a systematic difference to models assuming a periodic repetition of cylinders in x 1 -direction. This difference, however, is of the order of the variation of the results provided by the numerical and semi-analytical methods, cf. Table 3 and can therefore be considered minor. Due to the current design of the micromodel, the exact pressure drop of the main porous domain cannot be measured in the experiment, making a more detailed quantitative comparison with simulations of the benchmark problem impossible. Furthermore, despite a detailed grid convergence study there are non-quantifiable errors in any numerical methods that can contribute to the discrepancies with experiments. It is remarked that numerical dissipation in LBM, even at the scales of grid spacing, and the numerical dispersive effects are smaller compared to other second-order accurate methods (Marié et al. 2009). It can thus be inferred that the LBM simulation, which captures the whole experimental setup, achieves a marvellous agreement with the experimental measurements.

Comparison of Methods
The results from all applied approaches are collected in Table 3 and visualised in Fig. 11. It stands out that for r = 0.4 mm the 3D and pseudo 2D/3D homogenisation approaches as well as all pore-scale models provide similar permeability estimates.
An analysis of the homogenisation and the pore-scale approaches shows that the noslip boundary conditions in x 2 -direction do not influence the results, thus the width of the benchmark problem is chosen to be large enough. In contrast, the no-slip boundary condition in x 3 -direction has an influence on the results of the estimated permeability. This can be seen in the difference between the permeability values obtained from the Kozeny-Carman relation and the classical 2D homogenisation with the ones obtained by the other methods. The 2D homogenisation assumes no frictional effect from the top and bottom boundaries, which is obviously a crude assumption for the considered case. To evaluate the influence of the frictional effects, a corresponding FEM simulation for cylinders (r = 0.4 mm), periodically repeated in all directions (infinite height), yields a permeability of k = 1.9 × 10 −9 m 2 . This value is much closer to the value k = 3.9 × 10 −9 m 2 obtained from the Kozeny-Carman relation (2), cf. Sect. 3, and the value k = 1.8 × 10 −9 m 2 obtained from the 2D homogenisation, cf. Sect. 4.2, which do not account for the frictional effects of the surrounding walls.  . 11 Visualisation of the results collected in Table 3 As mentioned in Sect. 3, the Kozeny-Carman constant c kc can often be fitted in practical applications. If a fit of the parameter c kc towards the 3D homogenisation permeability values is carried out for this benchmark problem, then c kc = 167 gives a least-squares fit, or c kc = 133 minimises the 2-norm of the logarithm of the permeability values. However, even with fitted c kc one would find that the curve still has the wrong shape. Therefore, relationships like Kozeny-Carman fail to describe the evolution of the permeability even for simple porous media. The 3D homogenisation and the pseudo 2D/3D homogenisation approaches consider the frictional influence in x 3 -direction, while the other homogenisation methods do not. In particular, for the 3D approach no additional simplifications are made, making it the most accurate homogenisation approach which can be used as reference solution for the other homogenisation approaches. For decreasing pore-space passages (r → 0.5 mm), the frictional influence is reduced relatively and, thus, the permeability estimate of the classical 2D homogenisation approach tends towards the 3D one. Homogenisation approaches are computationally more efficient than pore-scale simulations but they are solely applicable for periodic microstructures. Pore-scale models provide accurate, detailed local information, such as 3D velocity and pressure fields also in the case of arbitrary geometries. We compare simulation results obtained by the pore-scale models and present in Fig. 12 the velocity profiles of FEM, SPH and LBM along the pore throat in the center of the domain (x 1 = 7.5 mm, x 2 = [4.9, 5.1] mm). The considered pore-scale methods deliver very similar velocity profiles. All models work with different driving forces, therefore the normalised and not the absolute velocities are compared.
For the considered benchmark problem, the use of a 1D Darcy's law is reasonable to estimate the scalar-valued intrinsic permeability. However, if certain 3D effects may play a role they can be extracted from pore-scale simulations in further studies, e. g. to estimate an anisotropic permeability tensor. As seen in other comparison studies using 3D CT scans of real rock samples, further uncertainties in terms of rock heterogeneities are expected to affect the estimated permeability results (Guibert et al. 2015;Zhang et al. 2019;Song et al. 2019).

Conclusions
The intrinsic permeability k of the benchmark problem (r = 0.4 mm) was estimated by the applicable approaches -classical 3D and pseudo 2D/3D homogenisation as well as all porescale models -in the range between 16.4 × 10 −11 m 2 and 17.7 × 10 −11 m 2 , cf. Table 3 and Fig. 11. In addition to the benchmark problem, further pore structures with different radii were studied and the results were discussed with regard to the individual strengths and limitations of the used methods. The experimentally determined permeability is underestimated in comparison with the computational approaches that consider the regular domain only and not the inlet and outlet channels between the location of pressure measurement and the porous domain, due to the design of the PDMS micromodel. However, it is consistent with our expectations and confirmed by the performed numerical LBM study of the microfluidic setup. The Kozeny-Carman relation does not provide a good estimate for the permeability. The 3D homogenisation gives accurate information about the permeability, while simplified 2D approaches have a limited range of applicability. Classical 2D homogenisation yields an overestimation of the permeability due to the wall effects on the top and bottom boundaries of the porous domain. Pore-scale resolved numerical simulations provide accurate estimates about the permeability, but are limited by computational costs.

Conflicts of interest There are no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.