Network-inspired versus Kozeny–Carman based permeability-porosity relations applied to Biot’s poroelasticity model

Water injection in the aquifer induces deformations in the soil. These mechanical deformations give rise to a change in porosity and permeability, which results in non-linearity of the mathematical problem. Assuming that the deformations are very small, the model provided by Biot’s theory of linear poroelasticity is used to determine the local displacement of the skeleton of a porous medium, as well as the fluid flow through the pores. In this continuum scale model, the Kozeny–Carman equation is commonly used to determine the permeability of the porous medium from the porosity. The Kozeny–Carman relation states that flow through the pores is possible at a certain location as long as the porosity is larger than zero at this location in the aquifer. However, from network models it is known that percolation thresholds exist, indicating that the permeability will be equal to zero if the porosity becomes smaller than these thresholds. In this paper, the relationship between permeability and porosity is investigated. A new permeability-porosity relation, based on the percolation theory, is derived and compared with the Kozeny–Carman relation. The strongest feature of the new approach is related to its capability to give a good description of the permeability in case of low porosities. However, with this network-inspired approach small values of the permeability are more likely to occur. Since we show that the solution of Biot’s model converges to the solution of a saddle point problem for small time steps and low permeability, we need stabilisation in the finite element approximation.


Introduction
For the description of different physical processes, such as consolidation, it is of a pivotal importance to have a valid estimation of permeability. The permeability of porous media is usually expressed as a function of some physical properties of the interconnected pore system such as porosity. Although it is natural to assume that the permeability depends on the porosity, it is not simple to formulate satisfactory theoretical models for the relation between them, mainly due to the complexity of the connected pore space. One of the most widely used permeability-porosity relationships is the Kozeny-Carman relation [13,24]. This relation assumes that the connectivity of a porous space does not vary in time, either by assuming a pore space that stays fully connected or by taking the effective porosity initially and assuming no loss of connectivity. Therefore, the Kozeny-Carman relation assumes that flow through the porous medium is possible as long as the porosity is nonzero. Hence, this relation is not capable of predicting blocking of the flow if the porosity is too low in some parts of the porous medium leading to two or more disconnected regions. Moreover, it is empirically proven that the permeability decreases dramatically with decreasing porosity [7], indicating that the Kozeny-Carman relation is less accurate at low porosities. To improve the behaviour of this relation for small values of the porosity, Mavko and Nur [26] incorporated a simple porosity adjustment into the Kozeny-Carman relation, by taking the percolation threshold into account. This approach resulted in a better prediction of the permeability for low porosities. However, in the work of Mavko and Nur the percolation threshold was chosen empirically to give a good fit for the experimental data. Another approach to take into account the percolation threshold was presented by Porter et al. [30], based on the work of Koltermann and Gorelick [22], who adapted the Kozeny-Carman equation to represent bimodal grain-size mixtures.
In the literature, the relationship between the permeability and the porosity is also derived using upscaling-based approaches, in case the underlying pore geometry of a representative elementary volume is prescribed [34]. Application of these approaches provides a link from the micro-to macroscopic behaviour of the material and allows the derivation of effective properties. Starting from mathematical models at the pore scale, such as the Stokes equation for fluid flow, an averaging procedure is performed in order to derive effective models. Possible averaging methods that may be applied to these equations are volume averaging [48] and two-scale asymptotic expansion [4]. The effective permeability tensor is given as the integral over solutions of auxiliary cell problems which are defined on the representative elementary volume. Upscaling methods directly enable calculating the full, potentially anisotropic, permeability tensor and one only needs the geometric information in terms of a representative elementary volume as input [34]. In contrast, wellestablished relations such as the Kozeny-Carman equation relate the porosity to a scalar permeability coefficient. Another upscaling method is the fast marching method (FMM), that is employed by Sharifi et al. [36] as an efficient tool for permeability upscaling. Their approach is based on the minimisation of the difference between dynamic behaviours of fine-scale and up-scaled reservoir models that is provided by the FMM. Upscaling approaches are also used to derive Biot's poroelasticity theory. Mikelić and Wheeler [27] derived the Biot-Allard equations via the use of homogenisation and asymptotic expansions in the space and time variables, the quasi-static Biot's system is then obtained by taking the singular limit of the more general dynamic Biot-Allard's system. In [42], the model for linear thermo-poroelasticity with nonlinear convection is derived using homogenisation. While the nonlinear quasi-static thermo-poroelastic equations are obtained using the two-scale asymptotic expansion method [12].
The Kozeny-Carman equation is based on only having spherical grains in the porous medium, whereas these grains can have various shapes. In this sense, the Kozeny-Carman relation represents a limit case. Another limit case is the assumption that the voids between the grains are represented by straight channels. In this study, we briefly introduce a new approach for the permeability that is derived on a micro-scale network model. At the pore-scale, this network model offers a detailed description of the porous medium [6]. The current modelling framework deals with the latter limit case, and can be used for general network topologies (such as rectangular, triangular and cubic arrangements of the channels). In contrast with the Kozeny-Carman equation, the new network-inspired approach yields that the permeability is nonzero only if the porosity is larger than a specific percolation threshold, that depends on the topology of the network. This means that the new approach may give a better prediction of the permeability in physics problems where abrupt changes in the porous medium occur resulting in a loss of connectivity of the connected pore space. In addition, the current model provides a computational framework to determine the percolation threshold for any given network with straight channels, which implies that there is no need for fitting on the basis of experiments if the network topology is known. This numerically determined threshold value agrees very well with the known values of the percolation thresholds from the literature. The new permeability-porosity approach is derived from percolation theory, which is a branch of probability that describes the effects of randomness arising from the porous media structure [11].
In percolation theory, the nodes of the network are called sites and the edges are called bonds. This leads to two approaches: site percolation and bond percolation. In this paper, we consider the bond percolation approach, whose basic idea can be explained as follows. Consider a network such as shown in Fig. 1. Whether a bond is open or closed depends on a certain probability and it is independent of the neighbouring bonds. If one bond is open and the nearest neighbours are closed, it is said that a 1-cluster is formed. If two adjacent bonds are open, they form a 2-cluster and so on. If the probability p of a bond being open increases, then larger clusters are created. There exists a critical probability at which the first cluster that spans the entire network from the inlet to the outlet is obtained. For infinite networks, this critical probability is well defined and it is called the percolation threshold p c [6].
Some macroscopic properties of porous media are mainly determined by the connectivity of the pore system [6]. Hence, it is possible to find a relation between the permeability of a porous medium and the percolation properties using the network model. For instance, if each of the open bonds conducts a fluid and there exists a cluster that spans the network from the inlet to the outlet, then a volumetric flow is possible through the network.
By adding more open bonds to the cluster, the volumetric flow through the network may increase. Since the number of open bonds represents the porosity of the network, it is possible to find a relation between the volumetric flow and the porosity [2,49]. In addition, Darcy's law gives a relationship between the permeability of the porous medium and the volumetric flow. As a result, a relation between the permeability and the porosity can be derived.
The applicability of these permeability-porosity relations will be demonstrated using two illustrative numerical examples. In these examples, flow of an incompressible fluid through a linearly elastic, saturated porous medium is modelled, using the classical theory of poroelasticity. Originally developed by Biot [8], poroelasticity theory assumes a superposition of solid and fluid components and couples the mechanical deformation of a porous solid with fluid flow through its internal structure. This approach is valid in the infinitesimal deformation range. Poroelasticity problems exist widely in the real world, making this theory of a great interest due to its applicability in various branches of science and engineering (see e.g. [5,14,32,37,39]).
In order to investigate the difference between the Kozeny-Carman approach and the equation based on the percolation theory, the boundary conditions in both academic poroelasticity examples are chosen such that a decrease in porosity is realised in some parts of the computational domain. This decrease in porosity will lead in both relations to a decrease in the permeability of the porous medium. When using the finite element method to solve the poroelastic equations, it is well known that the numerical solution of these equations may exhibit nonphysical oscillations in the pressure field for low permeabilities and short time steps [18,43]. Under these conditions, the resulting finite element discretisation approaches saddle point problems. Hence, a considerate numerical methodology in terms of possible spurious oscillations is needed. Therefore, a Galerkin finite element method based on Taylor-Hood elements has been developed, combined with the stabilisation technique proposed in [1,33]. Another stabilised finite element method is employed by Wan [44], Korsawe and Starke [23] and Tchonkova et al. [40], based on the Galerkin least-squares method.
Poroelasticity problems have been attracting attention from the scientific computing community [21,47] (and references therein). Another class of numerical methods that are used to solve the poroelasticity equations are the finite volume method combined with a nonlinear multigrid method as adopted by Luo et al. [25]. In addition, stabilised finite difference methods using central differences on staggered grids are used by Gaspar et al. [16,17] to solve Biot's model. However, the numerical solution of poroelasticity equations has been traditionally associated with finite element methods [3,29]. Furthermore, a monolithic approach for solving the quasi-static two-field poroelasticity equations is employed in this paper, which involves solving the coupled governing equations of flow and geomechanics simultaneously at every time step. Another approach that is widely used in coupling the flow and the mechanics in porous media is the fixed stress split method [10,28]. Hong et al. [19,20] applied the fixed stress splitting scheme to multiple-network poroelasticity systems.

Governing equations
The model provided by Biot's theory of linear poroelasticity with single-phase flow [8] is used in this study to determine the local displacement of the grains of a porous medium and the fluid flow through the pores, assuming that the deformations are very small. The fluid-saturated porous medium has a linearly elastic solid matrix and is saturated by an incompressible Newtonian fluid. Let Ω ⊂ R 2 denote the computational domain with boundary Γ , and x = (x, y) ∈ Ω. Furthermore, t denotes time, belonging to a half-open time inter-val I = (0, T], with T > 0. The initial boundary value problem for the consolidation process of an incompressible fluid flow in a deformable porous medium is stated as follows [1,45]: where σ and v are defined by the following equations In the above relations, σ and ε denote the effective stress and strain tensors, p the pore pressure, u the displacement vector, v Darcy's velocity, λ and μ the Lamé coefficients; κ the permeability of the porous medium and η the fluid viscosity. The parameter values used in this paper are given in Table 5. In addition, appropriate boundary and initial conditions are specified in Sect. 4.

The permeability-porosity relations
In this study, we consider the spatial dependency of the porosity and the permeability of the porous medium. The porosity θ is computed from the displacement vector using the porosity-dilatation relation (see [31,41]) with θ 0 the initial porosity. Subsequently, the permeability can be determined using the Kozeny-Carman equation [46] κ( where d s is the mean grain size of the soil. The Kozeny-Carman relation assumes that the permeability becomes zero if and only if the porosity also becomes zero. A new approach for the relation between the porosity and the permeability is inspired by the fluid flow through a network, where the fluid flows into the edges (channels) of the network. The network-inspired relation takes into account that a certain number of the channels are removed randomly, therefore the fluid can not flow through those particular channels causing an alteration in the permeability and the porosity of the network. The number of removed channels increases from 1% of the total number of channels to 100%. For a certain number of removed channels, there are no connected paths anymore between the inlet and the outlet. In this case, the fluid will stop flowing and the permeability will be expected to become zero. For an arbitrary network topology, the network-inspired relation states: where κ 0 is the initial permeability. The percolation threshold porosityθ = p c θ 0 , represents the minimal porosity needed to have connection via voids or channels from one end to the other, and is dependent on the topology of the network. Forθ = 0.5θ 0 , the normalised permeabilities (κ/κ 0 ) for the network-inspired relation and the Kozeny-Carman relation are depicted in Fig. 2 as function of the normalised porosity (θ /θ 0 ). In the coming section, we will show how the permeability-porosity relation (7) is derived using a network model.

The network-inspired permeability-porosity relation
In this section, we explain the procedure followed to obtain the network-inspired permeability-porosity relation. According to Balberg [2] and Wong [49], it is possible to obtain a relation between the number of open bonds (or channels) in the network N o and the flow rate through the network Q. This relation is determined by a power law as follows: in whichN is the number of bonds corresponding with the percolation threshold porosity and b is an exponent that can be determined by theory or via computer simulations. A similar relation between the permeability κ and the porosity θ will be obtained in this section. First, we consider a network with all channels open and with a pressure gradient p imposed over the horizontal direction of the network. Second, we assume that mass conservation holds in each node of the network n i , hence where S i = {j | node n j is adjacent to node n i }, and q ij is the flow rate in the channels connected to node n i . This flow rate is given by the Poiseuille flow where r and l are the radius and the length of a channel between two neighbouring nodes, respectively, and p ij is the pressure drop.
A linear system for the pressure p i in each node arises from substituting Eq. (10) in Eq. (9). This system is solved via a direct method and, subsequently, the flow rate in each channel is computed via Eq. (10). Finally, the flow rate through the network Q is computed by the summation of the flow rates in the channels that are connected with the outlet of the network. After this, we randomly close the channels, starting with 1% of the total number of channels in the network until that 100% of the channels is closed. In each stage, 500 simulations are performed and in each simulation, the linear system for the pressure is solved and the flow rate through the network is computed. From the computed flow rate Q, the permeability of the network can be determined using Darcy's law where A is the cross-sectional area of the network and L is the length over which the pressure gradient is taking place. In addition, there is a direct relation between the porosity of the network and the volume of the open channels V o in which V t is the total volume of the channels. This procedure yields a relation between the permeability and the porosity of the network. In the coming sections, this procedure will be demonstrated for three two-dimensional networks: rectangular, triangular and triangular unstructured.

Rectangular network
We start by describing the results obtained for a rectangular network (see Fig. 1). In this case, we use a network with N x = 100 horizontal nodes and N y = 60 vertical nodes. We computed the fraction of closed channels f c = V c /V t , that is needed to obtain a normalised permeability κ n = κ/κ 0 such that k i -0.05 < κ n < k i + 0.05, where k i ∈ {0.1, 0.2, . . . , 0.9}. In  Table 1.

Triangular network
In this section, we present the results obtained with a triangular network as shown in Fig. 4. The number of nodes in the horizontal direction N x = 100 and the number of nodes in the vertical direction N y = 60. The coordination number of the interior nodes is eight or four (see Fig. 4). In Fig. 5, the histograms for k i = 0.1, k i = 0.5 and k i = 0.9 are depicted. In addition, the values of the meanμ and the standard deviation σ of the distribution of f c are presented in Table 2.

Triangular unstructured network
Finally, we present the results obtained with a triangular unstructured network such as shown in Fig. 6. The total number of nodes in this network is 6921. The histograms of the fraction of closed channels for some values of the normalised permeability are depicted in Fig. 7. Moreover, in Table 3, the values of the meanμ and the standard deviation σ of the distribution of f c are given. This method is also used to compute f c for κ n = 0. The meanμ of the smallest value of f c for which holds κ n = 0 is depicted in Table 4. Since θ /θ 0 = 1f c , we observe from the results obtained with the different networks that the permeability becomes zero for porosities below a certain threshold. These percolation thresholds are different for the different networks used in this study, as can be seen in Table 4. The values are in good agreement with the literature [38]. We also observe that, for all three networks, the relation between the permeability and the porosity exhibits an almost linear increase for values of the porosity larger than the percolation threshold. Moreover, the slope of these linear lines  Table 2 The meanμ and the standard deviation σ of the distribution of f c for the triangular network

Problem with high pump pressure
In this numerical experiment, the infiltration of an incompressible fluid through a filter into a two-dimensional area is considered, as shown in Fig. 9. During the infiltration, a high pump pressure is used to inject water into the porous medium. Leading to a compression of the material against the rigid right boundary Γ 4 .
The computational domain Ω is a two-dimensional rectangular surface with Cartesian coordinates x = (x, y). In order to solve this problem, Biot's consolidation model, as described in Sect. 2, is applied on the computational domain Ω with width L and height H.
The fluid is injected into the soil through a filter placed on boundary segment Γ 2 . More precisely, the boundary conditions for this problem are given as follows: where t is the unit tangent vector at the boundary, n the outward unit normal vector and p pump is a prescribed high pump pressure due to the injection of the fluid. Figure 9 shows the definition of the boundary segments. Initially, the following condition is imposed:

Squeeze problem
The infiltration of a fluid through a filter into a rectangular two-dimensional area is shown in Fig. 10. In this numerical experiment, the porous medium is squeezed by applying a vertical load on the middle of the top and bottom edges of the domain. This problem is solved using Biot's consolidation model, that is applied on the computational domain Ω with width L and height H. The fluid is injected into the soil through a filter placed on boundary segment Γ 3 . The boundary conditions for this problem are given as follows: where t is the unit tangent vector at the boundary, n the outward unit normal vector, p pump is a prescribed pump pressure due to the injection of the fluid and σ 0 is the intensity of a uniform vertical load. Figure 10 shows the definition of the boundary segments. Initially, the following condition is imposed:
In addition, we consider the bilinear forms where the notation (17) is used for the bilinear form c(p, q).The variational formulation for problem (1a)-(1b) with boundary and initial conditions (13a)-(13f)- (14), consists of the following, using the no- with the initial condition u(0) = 0, and where Subsequently, problem (18a)-(18b) is solved by applying the finite element method. Let P k h ⊂ H 1 (Ω) be a function space of piecewise polynomials on Ω of degree k. Hence, we define finite element approximations for W and Q as . . , n p }, respectively [1]. Afterwards, we approximate the functions u(t) and p(t) with functions u h (t) ∈ W k h and p h (t) ∈ Q k h , defined as in which the Dirichlet boundary conditions are imposed. Simultaneously, discretisation in time is applied using the backward Euler method. Let τ be the time step size and define a time grid {t m = mτ : m ∈ N}, then the discrete Galerkin scheme of (18a)-(18b) is formulated as follows: while for m = 0: u 0 h = 0. The discrete Galerkin scheme for problem (1a)-(1b) with boundary and initial conditions (15a)-(15g)-(16) is derived similarly.
In the network-inspired relation (7), the permeability is equal to zero for θ ∈ [0,θ ]. Hence for high percolation thresholds, is it more likely that the permeability goes to zero. From an intuitive point of view, it can be seen that, as τ κ → 0, the solution to the semi-discrete problem inherits the properties of a Stokes-like equation, for which it becomes favourable to use stabilisation techniques (possibly in combination with Taylor-Hood, mini-elements, or Crouzeix-Raviart elements).
Since Biot's poroelasticity problem converges to a saddle point problem as τ κ → 0, the Taylor-Hood elements can be used to solve this problem numerically. These elements represent finite element pairs that satisfy the inf-sup condition for the saddle point problem [9,15]. Although the inf-sup condition and the coercivity and boundedness of bilinear form a(u, w) warrant existence and uniqueness of the finite element solution, the inf-sup condition is not sufficient for reliable numerical solutions of Biot's problem (1a)-(1b). Since for low permeabilities and/or small time steps, the approximation to the pressure exhibits nonphysical oscillations due to loss of the M-matrix property, as shown in [1]. In order to improve the monotonicity properties of the finite element scheme and to obtain oscillation free approximations of the pressure, the stabilisation procedure outlined in [1,33] is applied in this study. Therefore, a stabilisation term is added to the continuity equation in Biot's model with stabilisation parameter β = √ x 2 + y 2 4(λ+2μ) .

Numerical results
The Galerkin finite element method, with triangular Taylor-Hood elements [31,35], is adopted to solve the discretised quasi-two-dimensional problem (1a)-(1b). The displacements are spatially approximated by quadratic basis functions, whereas a continuous piecewise linear approximation is used for the pressure field. From the system of equations (1a)-(1b) and the boundary conditions (13a)-(13f) it is obvious that this two-dimensional problem is symmetrical in the y-direction and that it can be reduced to a one-dimensional problem. Aguilar et al. [1] solved this one-dimensional problem analytically and showed that the finite element method with Taylor-Hood elements gives accurate numerical results. For the time integration, the backward Euler method is applied. The numerical investigations are carried out using the matrix-based software package MATLAB (version R2015a).
The computational domain is a rectangular surface with width L = 2.0 m and height H = 1.0 m. The domain is discretised using a regular triangular grid, with x = y = 0.02. In addition, values for some model parameters have been chosen based on literature (see Table 5).
Furthermore, the Lamé coefficients λ and μ are related to Young's modulus E and Poisson's ratio ν by  The impact of the permeability-porosity relation on the water flow is defined in this study as the impact on the time average of the volumetric flow rate Q out at a distance L from the injection filter. The initial permeability κ 0 used in Eq. (7) is computed by the Kozeny-Carman relation (6), in order to have a reliable comparison between the two relations. In the generations of the simulation results, the time step size is chosen to be τ = 0.5.

Numerical results for the problem with high pump pressure
In order to obtain some insight into the impact of a high pump pressure on the water flow, we present an overview of the simulation results in Figs. 11-13. In these simulations, water is injected into the soil at a constant pump pressure of 50 bar. The simulated fluid velocity, permeability and porosity profiles that have been obtained using the Kozeny-Carman relation are provided in Fig. 11, while the simulated results that have been obtained using the network-inspired relation with p c = 0.3232, corresponding with a triangular structured network, are provided in Fig. 12. In Fig. 13, the simulated results that have been obtained using the network-inspired relation with p c = 0.4935, corresponding with a rectangular network, are depicted. As shown in these figures, the injected water flows in the horizontal direction through the domain from the inlet to the outlet. Furthermore, the magnitude of the velocity |v| in the inlet is larger than in the outlet. Note that due to the continuity equation (1b), the fluxes at the inlet and the outlet are not necessarily equal. The change in |v| can be explained by the permeability profiles shown in Figs. 11(b)-13(b). In these figures we observe that the permeability decreases almost linearly from the inlet to the outlet. In addition, the permeability obtained using the Kozeny-Carman relation exhibits a larger decrease than the permeabilities obtained with the network-inspired model. The normalised permeability using the Kozeny-Carman relation decreases from 1 to 0.4659, while the normalised permeabilities for the network-inspired relation decrease from 1 to 0.7519 and from 1 to 0.6684 for the triangular structured and the rectangular network respectively. This behaviour is clarified by Fig. 8 and Figs. 11(c)-13(c). Due to boundary condition (13e), the values of the normalised porosity in all three cases are almost the same in the outlet and are equal to 0.8321. In Fig. 8, we see that for this value of the porosity, the normalised permeability is the lowest for the Kozeny-Carman relation and the highest for the network-inspired relation derived from the triangular structured network. This explains the difference in decrease in the permeability profiles. Note that the values of the material properties in Table 5 belong to a porous medium consisting of sand grains, hence a high pump pressure will apparently not result in hydraulic fracturing. This study therefore neglects this phenomenon.
In Fig. 14, the time average of the volumetric flow rate Q out is depicted for different values of the percolation threshold. As expected from Fig. 8, for low percolation thresholds To demonstrate this, the problem is solved for a coarser grid x = y = 0.04, as shown in Fig. 15. In this figure, it becomes even worse, and we observe spurious oscillations for the network-inspired relation for percolation thresholds larger than 0.85 approximately. The reason for this behaviour is that for large values of the percolation threshold, the permeability goes to zero very soon. Therefore, even the used stabilisation did not diminish the nonphysical oscillations. Probably a stronger stabilisation will alleviate these spurious oscillations. This was behind the scope of the current paper.

Numerical results for the squeeze problem
The impact of the imposed vertical load on boundary segments Γ 1 and Γ 5 is shown in Figs. 16-18, using the Kozeny-Carman relation and the network-inspired relation respectively. In these simulations, water is injected into the porous medium at a constant pump pressure equal to 5.0 bar. The simulated fluid velocity, permeability and porosity profiles that are obtained using the Kozeny-Carman relation are provided in Fig. 16, while the simulated results that are obtained using the network-inspired relation with p c = 0.3232, corresponding with a triangular structured network, are provided in Fig. 17. In Fig. 18, the simulated results that are obtained using the network-inspired relation with p c = 0.4935, corresponding with a rectangular network, are depicted. In Fig. 19, the time average of the volumetric flow rate Q out is depicted for different values of the percolation threshold. Similarly to the high pump pressure problem, the flow rates for small values of the percolation threshold using the network-inspired relation are higher than the flow rates obtained using the Kozeny-Carman relation. In addition, we observe that the flow rate depends significantly on the percolation threshold and hence on the topology of the network for high percolation thresholds.

Figure 19
The time average of the volumetric flow rate Q out as a function of the percolation threshold p c

Discussion and conclusions
In this paper, the network-inspired permeability-porosity relation is applied on two poroelasticity problems. This numerical experiment is designed in order to analyse the applicability of this microscopic relation on the macro-scale. Furthermore, we compare the results obtained with the network-inspired relation to the Kozeny-Carman relation which is often used in these physical problems. In the first problem, a high pump pressure is imposed in the inlet of a porous medium package. This high pressure forces the grains to move towards the outlet. In the second problem the package is squeezed by applying a load on the middle of the top and bottom edges of the domain. The purpose of considering these poroelasticity problems is to create a large density of the grains in the computational domain which results in a decrease of the porosity. In these problems, Biot's model for poroelasticity is used to determine the water pressure and the displacements of the grains that are needed to compute the porosity. From the porosity the permeability is determined either by the network-inspired relation or by the Kozeny-Carman relation.
Depending on the topology, three different percolation thresholds, corresponding with a rectangular network (p c = 0.4935), triangular structured network (p c = 0.3232) and triangular unstructured network (p c = 0.3438), are distinguished. However, since the topology of macro-scale porous media is not known, computations are also performed with percolation thresholds in the interval [0, 0.975] to investigate the influence of the percolation threshold (and hence the topology of the porous medium) on the flow rate. First, the problems are solved with the Kozeny-Carman relation, the network-inspired relation based on the triangular structured network and the relation based on the rectangular network. From the numerical results we conclude that the permeability obtained using the Kozeny-Carman relation exhibits a larger decrease than the permeabilities obtained with the network-inspired relations, which is clarified by Fig. 8. In contrast, the porosity profile is not affected significantly by the selected permeability-porosity relation. Second, the time average of the volumetric flow rate was computed for percolation thresholds in the interval [0, 0.975]. For low percolation thresholds the network-inspired relation results in higher flow rates than the Kozeny-Carman relation, as expected from Fig. 8. In addition, it is shown that the flow rate changes significantly as a function of the percolation threshold which means that the water flow depends on the topology of the network. For high percolation thresholds, spurious oscillations appeared due to the violation of the Mmatrix property in the discretisation matrix that resulted from the convergence of Biot's problem to the related saddle point problem. The results for these percolation thresholds could be improved by using a finer grid.
For the studied problems and the set of parameters chosen, we noticed that the applied permeability-porosity relations result in small changes in the porosity while a major change is realised in the permeability profiles. A possible explanation for this behaviour is that the relation between the velocity field and the change of the displacements in time as stated in Eq. (1b), is not strong enough to lead to significant changes in the porosity profile.