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 [6], 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. [29], based on the work of Koltermann and Gorelick [22], who adapted the Kozeny-Carman equation to represent bimodal grain-size mixtures.
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 [5]. 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 networkinspired 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 [12].
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 [5].
Channel Node Figure 1: An illustration of a rectangular network.
Some macroscopic properties of porous media are mainly determined by the connectivity of the pore system [5]. 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,45]. 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 [7], 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. [4,14,31,35,37]).
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,40]. In Sect. 5.1, we will prove that 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,32]. Another stabilised finite element method is employed by Wan [41], Korsawe and Starke [23] and Tchonkova et al. [38], based on the Galerkin least-squares method.
Poroelasticity problems have been attracting attention from the scientific computing community [21,44] (and references therein). Another 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,28]. 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 [9,27]. 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 [7] 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 interval 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,42]: where σ and v are defined by the following equations σ = λtr(ε)I + 2µε; (2) 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 [30,39]) with θ 0 the initial porosity. Subsequently, the permeability can be determined using the Kozeny-Carman equation [43] κ(x, t) = d 2 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 [45], 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 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 Fig. 3, the histograms for k i = 0.1, k i = 0.5 and k i = 0.9 are shown. The meanμ and the standard deviation σ of the distribution of f c depend on the normalised permeability κ n as presented 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.
Channel Node Figure 4: An illustration of a triangular network.

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 = 1 − f 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   [36]. 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 depends on the percolation threshold and hence on the network topology. The normalised permeabilities, for the different network topologies used in this paper and for the Kozeny-Carman equation, are depicted in Fig. 8.

Problem formulation
The following numerical experiments are designed to study the different relations for the porosity and the permeability. In both experiments, the boundary conditions are chosen such that a decrease in porosity is realised in some parts of the computational domain.

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 Figure 9: Sketch of the setup for the two-dimensional problem with high pump pressure.
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 Figure 10: Sketch of the setup for the two-dimensional squeeze problem.
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:

Numerical method
To present the variational formulation of problem (1), we first introduce the appropriate function spaces. Let L 2 (Ω) be the Hilbert space of square integrable scalar-valued functions, with inner product (f, g) = Ω f g dΩ. And let H 1 (Ω) denote the subspace of L 2 (Ω) of functions with first derivatives in L 2 (Ω). We further introduce the function spaces Q(Ω) = {q ∈ H 1 (Ω) : q| Γ2 = p pump and q| Γ4 = 0}.
In addition, we consider the bilinear forms The variational formulation for problem (1) with boundary and initial conditions (13)- (14), consists of the following, using the notationu = ∂u ∂t : Find (u(t), p(t)) ∈ (W × Q) such that with the initial condition u(0) = 0, and where Subsequently, problem (17) 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 u } and Q k h = Q ∩ P k h with basis {ψ j ∈ Q k h : j = 1, . . . , 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 (17) is formulated as follows: while for m = 0: u 0 h = 0. The discrete Galerkin scheme for problem (1) with boundary and initial conditions (15)-(16) is derived similarly.
In the network-inspired relation (7), the permeability is equal to zero for θ ∈ [0,θ]. Hence for large percolation thresholds, is it more likely that the permeability goes to zero. In the next section, the convergence of problem (1) to a saddle point problem when κ goes to zero is proven.

Convergence to a saddle point problem
In this section, we will prove that the solution of the variational formulation (17) converges to the solution of the related saddle point problem as τ κ → 0, both in the continuous as in the discrete system. Since we are interested in the case τ κ → 0, we can assume without loss of generality that κ is bounded and that it fulfils the following relation (see Fig. 8): where κ max > 0 is constant and f (θ) ∈ [0, 1]. Consequently, define the bilinear form and define u = u m − u m−1 , which gives u m = u + u m−1 . This transformation results in the following variational formulation: For m ≥ 1, find ( u, p m ) ∈ (W × Q) such that The well-posedness of this problem has been analysed in [34,46]. Using the inequality |2ab| ≤ (a 2 + b 2 ) for a, b ∈ R, we can show that where the second inequality is a result of the definition of the H 1 -norm Furthermore, define the spaces R(Ω) = {q ∈ L 2 (Ω) : q| Γ2 = p pump and q| Γ4 = 0} and R 0 (Ω) = {q ∈ L 2 (Ω) : q| Γ2∪Γ4 = 0}. Hence, assuming that τ κ max → 0, the variational formulation becomes: For m ≥ 1, find ( u 0 , p m 0 ) ∈ (W × R) such that Brenner and Scott [11] showed, using Korn's inequality, that the bilinear form a(u, w) is coercive in W . In addition, the bilinear form b(p, w) satisfies the inf-sup condition on (W × R) as proven in [10]. The continuity of both bilinear forms follows from the Cauchy-Schwarz inequality. Hence, from Brezzi's splitting theorem [10] we can conclude that a unique solution exists for problem (P 0 ). This problem is referred to as a saddle point problem. We assume that our domain allows classical solutions to exist for problem (P 0 ), and as a result the solution of this problem is in the spaces (W × Q). Then we prove that solutions to (P τ κ ) converge to solutions of (P 0 ), which illustrates that the solutions of (P τ κ ) inherit the saddle point structure of (P 0 ). Theorem 1. The solution of (P τ κ ) converges to the solution of (P 0 ) for τ κ max → 0.
Proof. Subtracting (P 0 ) from (P τ κ ) yields Using (25b) and the bilinearity of the form c θ (p, q) gives From the coercivity of a(u, w) and the definition of c θ (p, q), we can state that where the last inequality follows from (23). Let τ κ max → 0, then follows that a( u − u 0 , u − u 0 ) → 0. Consequently, the coercivity of a(u, w) implies that u → u 0 as τ κ max → 0. Therefore, we have a( u − u 0 , w) → 0 for all w ∈ W . Subsequently, this last result together with Eq. (25a) lead to b(p m − p m 0 , w) → 0. Hence, we can conclude that p m → p m 0 as τ κ max → 0.
Proof. The convergence proof in the discrete problem is similar to the continuous problem. Hence, it is sufficient to show that a(u, w) is coercive in W k h . Since W k h ⊂ W , the coercivity property is automatically verified in the discrete case.
Since Biot's poroelasticity problem converges to the related 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 [8,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 (1). 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,32] 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 [30,33], is adopted to solve the discretised quasi-two-dimensional problem (1). 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 (1) and the boundary conditions (13) 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 onedimensional 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). p pump 50 · 10 5 / 5 · 10 5 Pa Uniform load σ 0 3 · 10 6 N/m 2 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. 11b-13b. 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. 11c-13c. 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 the network-inspired relation results in higher flow rates than the Kozeny-Carman relation. Furthermore, the flow rate changes significantly as a function of the percolation threshold. Hence the water flow depends on the topology of the network. The negative values of the flow rate Q out for the very high values of the percolation threshold from the network-inspired relation are caused by violation of the M-matrix property that gives loss of monotonicity of the numerical solution. 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 networkinspired 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 networkinspired 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 Figs. 16a-18a, the impact of the imposed vertical load on the computational domain is shown. The magnitude of the velocity is small near the boundary segments where the vertical load is applied and large in the middle near the symmetry axis y = H/2. In the outlet, the magnitude of the velocity is maximal near the upper and lower boundary segments and decreases towards the symmetry axis. In all three cases, the fluid flows mainly in the horizontal direction. In the region where the load is imposed, the domain is squeezed, resulting in a larger density of the grains. This leads to a lower porosity in this region, with minimum values 0.8809, 0.8811 and 0.8810 for the Kozeny-Carman relation, the triangular structured network-inspired relation and the rectangular networkinspired relation respectively. The abrupt transition in the boundary condition between boundary segments Γ 1 and Γ 2 (and between Γ 4 and Γ 5 ), which results in a discontinuous force, results in a small numerical artefact at the location of these transitions. As expected from Fig. 8 and the minimum values for the porosities, the decrease in permeability by the Kozeny-Carman relation, as shown in Fig. 16b, is larger than the decrease in the porosities obtained using the network-inspired relation, Figs. 17b and 18b.
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 large percolation thresholds.

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 networkinspired 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 networkinspired 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 large percolation thresholds, spurious oscillations appeared due to the violation of the M-matrix property in the discretisation matrix that resulted from the convergence of Biot's problem to the related saddle point problem, as proven in Sect. 5.1. 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.