    # 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.

## Authors

09/13/2021

### Adaptive iterative linearised finite element methods for implicitly constituted incompressible fluid flow problems and its application to Bingham fluids

In this work, we introduce an iterative linearised finite element method...
06/16/2022

### Well-posedness and variational numerical scheme for an adaptive model in highly heterogeneous porous media

Mathematical modeling of fluid flow in a porous medium is usually descri...
12/02/2020

### Homogenization of large deforming fluid-saturated porous structures

The two-scale computational homogenization method is proposed for modell...
09/01/2020

### A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems

In this paper, we develop a simplified hybrid weighted essentially non-o...
11/05/2020

### The Darcy problem with porosity depending exponentially on the pressure

We consider the flow of a viscous incompressible fluid through a porous ...
06/21/2022

### Continuous Data Assimilation for Displacement in a Porous Medium

In this paper we propose the use of a continuous data assimilation algor...
02/09/2018

### A Viscoelastic Catastrophe

We use a differential constitutive equation to model the flow of a visco...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 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 , 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  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. , based on the work of Koltermann and Gorelick , 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 

. 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

.

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 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  .

Some macroscopic properties of porous media are mainly determined by the connectivity of the pore system . 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 , 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 , Korsawe and Starke  and Tchonkova et al. , 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. . 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.

## 2 Governing equations

The model provided by Biot’s theory of linear poroelasticity with single-phase flow  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 denote the computational domain with boundary , and . Furthermore, denotes time, belonging to a half-open time interval , with . 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]:

 −∇⋅σ′+∇p=0 \ on\ Ω×I; (1a) ∂∂t(∇⋅u)+∇⋅v=0 \ on\ Ω×I, (1b)

where and are defined by the following equations

 σ′ =λtr(ε)I+2με; (2) ε =12(∇u+∇uT); (3) v =−κη∇p. (4)

In the above relations, and

denote the effective stress and strain tensors,

the pore pressure,

the displacement vector,

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.

## 3 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])

 θ(x,t)=1−1−θ0exp(∇⋅u), (5)

with the initial porosity. Subsequently, the permeability can be determined using the Kozeny-Carman equation 

 κ(x,t)=d2s180θ(x,t)3(1−θ(x,t))2, (6)

where 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 of the total number of channels to . 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:

 κ(x,t)=⎧⎨⎩0θ<^θθ−^θθ0−^θκ0θ≥^θ, (7)

where is the initial permeability. The percolation threshold porosity , 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 , the normalised permeabilities () for the network-inspired relation and the Kozeny-Carman relation are depicted in Fig. 2 as function of the normalised porosity (). In the coming section, we will show how the permeability-porosity relation (7) is derived using a network model. Figure 2: The normalised permeability as function of the normalised porosity for the Kozeny-Carman relation and the network-inspired relation, with ^θ=0.5θ0.

### 3.1 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  and Wong , it is possible to obtain a relation between the number of open bonds (or channels) in the network and the flow rate through the network . This relation is determined by a power law as follows:

 Q∝(No−^N)bforNo>^N, (8)

in which is the number of bonds corresponding with the percolation threshold porosity and 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 imposed over the horizontal direction of the network. Second, we assume that mass conservation holds in each node of the network , hence

 ∑j∈Siqij=0, (9)

where , and is the flow rate in the channels connected to node . This flow rate is given by the Poiseuille flow

 qij=πr48ηlΔpij, (10)

where and are the radius and the length of a channel between two neighbouring nodes, respectively, and is the pressure drop.

A linear system for the pressure 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 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 of the total number of channels in the network until that 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 , the permeability of the network can be determined using Darcy’s law

 κ=−QAηLΔp, (11)

where is the cross-sectional area of the network and 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

 θ=VoVtθ0, (12)

in which 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.

#### 3.1.1 Rectangular network

We start by describing the results obtained for a rectangular network (see Fig. 1). In this case, we use a network with horizontal nodes and vertical nodes. We computed the fraction of closed channels , that is needed to obtain a normalised permeability such that , where . In Fig. 3, the histograms for , and are shown. The mean

and the standard deviation

of the distribution of depend on the normalised permeability as presented in Table 1.

#### 3.1.2 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 and the number of nodes in the vertical direction . The coordination number of the interior nodes is eight or four (see Fig. 4). In Fig. 5, the histograms for , and are depicted. In addition, the values of the mean and the standard deviation  of the distribution of are presented in Table 2.

#### 3.1.3 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 . 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 are given.

This method is also used to compute for . The mean of the smallest value of for which holds is depicted in Table 4. Since , 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 . 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.

## 4 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.

### 4.1 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 .

The computational domain is a two-dimensional rectangular surface with Cartesian coordinates . In order to solve this problem, Biot’s consolidation model, as described in Sect. 2, is applied on the computational domain  with width and height . The fluid is injected into the soil through a filter placed on boundary segment . More precisely, the boundary conditions for this problem are given as follows:

where is the unit tangent vector at the boundary, the outward unit normal vector and 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:

 u(x,0)=0forx∈Ω. (14)

### 4.2 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 and height . The fluid is injected into the soil through a filter placed on boundary segment . The boundary conditions for this problem are given as follows:

where is the unit tangent vector at the boundary, the outward unit normal vector, is a prescribed pump pressure due to the injection of the fluid and is the intensity of a uniform vertical load. Figure 10 shows the definition of the boundary segments. Initially, the following condition is imposed:

 u(x,0)=0forx∈Ω. (16)

## 5 Numerical method

To present the variational formulation of problem (1), we first introduce the appropriate function spaces. Let be the Hilbert space of square integrable scalar-valued functions, with inner product . And let denote the subspace of of functions with first derivatives in . We further introduce the function spaces

 W(Ω) ={w∈(H1(Ω))2:(w⋅n)|Γ1∪Γ3∪Γ4=0}; Q(Ω) ={q∈H1(Ω):q|Γ2=ppump\ and\ q|Γ4=0}.

In addition, we consider the bilinear forms

 a(u,w) =λ(∇⋅u,∇⋅w)+2μ2∑i,j=1(ϵij(u),ϵij(w)); b(p,w) =(p,∇⋅w); c(p,q) =2∑i=1(κη∂p∂xi,∂q∂xi).

The variational formulation for problem (1) with boundary and initial conditions (13)-(14), consists of the following, using the notation : Find such that

 a(u(t),w)−b(p(t),w) =h(w) ∀w∈W; (17a) b(q,˙u(t))+c(p(t),q) =0 ∀q∈Q0, (17b)

with the initial condition and where

 h(w) =−ppump∫Γ2w⋅n dΓ; Q0(Ω) ={q∈H1(Ω):q|Γ2∪Γ4=0}.

Subsequently, problem (17) is solved by applying the finite element method. Let be a function space of piecewise polynomials on of degree . Hence, we define finite element approximations for and as with basis and with basis , respectively . Afterwards, we approximate the functions and with functions and , defined as

 uh(t)=nu∑i=1ui(t)ϕi,ph(t)=np∑j=1pj(t)ψj, (18)

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 , then the discrete Galerkin scheme of (17) is formulated as follows: For , find such that

 a(umh,wh)−b(pmh,wh) =h(wh) ∀wh∈Wkh; (19) b(qh,umh)+τc(pmh,qh) =b(qh,um−1h) ∀qh∈Qk′0h, (20)

while for : . 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 . 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.

### 5.1 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 , both in the continuous as in the discrete system. Since we are interested in the case , we can assume without loss of generality that is bounded and that it fulfils the following relation (see Fig. 8):

 κ(x,t)=κmaxf(θ), (21)

where is constant and . Consequently, define the bilinear form

 cθ(p,q)=2∑i=1(f(θ)η∂p∂xi,∂q∂xi), (22)

and define , which gives . This transformation results in the following variational formulation: For , find such that

 (Pτκ):={a(˜u,w)−b(pm,w)=g(u,w)∀w∈W;b(q,˜u)+τκmaxcθ(pm,q)=0∀q∈Q0,

where . The well-posedness of this problem has been analysed in [34, 46].

Using the inequality for , we can show that

 \abscθ(pm,q) ≤12η[2∑i=1(∂pm∂xi,∂pm∂xi)+2∑i=1(∂q∂xi,∂q∂xi)] ≤12η(∥pm∥2H1(Ω)+∥q∥2H1(Ω))<∞, (23)

where the second inequality is a result of the definition of the -norm

 ∥p∥2H1(Ω)=(p,p)+2∑i=1(∂p∂xi,∂p∂xi). (24)

Furthermore, define the spaces and . Hence, assuming that , the variational formulation becomes: For , find such that

 (P0):={a(˜u0,w)−b(pm0,w)=g(u,w)∀w∈W;b(q,˜u0)=0∀q∈R0.

Brenner and Scott  showed, using Korn’s inequality, that the bilinear form is coercive in . In addition, the bilinear form satisfies the inf-sup condition on as proven in . The continuity of both bilinear forms follows from the Cauchy-Schwarz inequality. Hence, from Brezzi’s splitting theorem  we can conclude that a unique solution exists for problem . This problem is referred to as a saddle point problem. We assume that our domain allows classical solutions to exist for problem , and as a result the solution of this problem is in the spaces . Then we prove that solutions to converge to solutions of , which illustrates that the solutions of inherit the saddle point structure of .

###### Theorem 1.

The solution of converges to the solution of for .

###### Proof.

Subtracting from yields

 a(˜u−˜u0,w)−b(pm−pm0,w) =0 ∀w∈W; (25a) b(q,˜u−˜u0)+τκmaxcθ(pm,q) =0 ∀q∈Q0. (25b)

Take in (25a), hence we get

 a(˜u−˜u0,˜u−˜u0)−b(pm−pm0,˜u−˜u0)=0. (26)

Using (25b) and the bilinearity of the form gives

 a(˜u−˜u0,˜u−˜u0)+τκmaxcθ(pm,pm)=τκmaxcθ(pm,pm0).

From the coercivity of and the definition of , we can state that

 0 ≤a(˜u−˜u0,˜u−˜u0)+τκmaxcθ(pm,pm) (27a) =τκmaxcθ(pm,pm0)<∞, (27b)

where the last inequality follows from (23). Let , then follows that . Consequently, the coercivity of implies that as . Therefore, we have for all . Subsequently, this last result together with Eq. (25a) lead to . Hence, we can conclude that as . ∎

Numerically, the discrete Galerkin scheme (19)-(20) can be expressed as a linear equations system

 (PDτκ):={Aum−BTpm=h,Bum+τκmaxCpm=Bum−1,

with and . The matrix is the symmetric positive definite elasticity matrix, the divergence matrix and is the diffusive matrix. The vector is the right-hand side vector with components . Initially, . For , we have

 (PD0):={Aum0−BTpm0=h,Bum0=Bum−10.
###### Theorem 2.

The solution of converges to the solution of as .

###### Proof.

The convergence proof in the discrete problem is similar to the continuous problem. Hence, it is sufficient to show that is coercive in . Since , the coercivity property is automatically verified in the discrete case. ∎

Since Biot’s poroelasticity problem converges to the related saddle point problem as , 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 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 . 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 .

## 6 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 -direction and that it can be reduced to a one-dimensional problem. Aguilar et al.  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 and height . The domain is discretised using a regular triangular grid, with . 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 and Poisson’s ratio by

 λ=νE(1+ν)(1−2ν),μ=E2(1+ν). (28)

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 at a distance  from the injection filter. The initial permeability 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 .

### 6.1 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 . 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 , 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 , 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 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 can be explained by the permeability profiles shown in Figs. (b)b-(b)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. (c)c-(c)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 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 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 , 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 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. Figure 14: The time average of the volumetric flow rate ¯¯¯¯Qout as a function of the percolation threshold pc, using Δx=Δy=0.02. Figure 15: The time average of the volumetric flow rate ¯¯¯¯Qout as a function of the percolation threshold pc, using Δx=Δy=0.04.

### 6.2 Numerical results for the squeeze problem

The impact of the imposed vertical load on boundary segments  and  is shown in Figs. 1618, 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 . 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 , 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 , corresponding with a rectangular network, are depicted.

In Figs. (a)a-(a)a, 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 . 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 network-inspired relation respectively. The abrupt transition in the boundary condition between boundary segments and (and between and ), 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. (b)b, is larger than the decrease in the porosities obtained using the network-inspired relation, Figs. (b)b and (b)b.

In Fig. 19, the time average of the volumetric flow rate 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. Figure 19: The time average of the volumetric flow rate ¯¯¯¯Qout as a function of the percolation threshold pc.

## 7 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 (), triangular structured network () and triangular unstructured network (), are distinguished. However, since the topology of macro-scale porous media is not known, computations are also performed with percolation thresholds in the interval 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 . 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.

##### Acknowledgements

The work of M. Rahrah was supported by the Netherlands Organisation for Scientific Research NWO (project number 13263). The work of L.A. Lopez-Peña was supported by the Mexican Institute of Petroleum (IMP) through the Programa de Captación de Talento, Reclutamiento, Evaluación y Selección de Recursos Humanos (PCTRES) grant.

## References

•  G. Aguilar, F. Gaspar, F. Lisbona, and C. Rodrigo (2008) Numerical stabilization of Biot’s consolidation model by a perturbation on the flow equation. Int. J. Numer. Meth. Engng. 75, pp. 1282 – 1300. Cited by: §1, §2, §5.1, §5, §6.
•  I. Balberg (1987) Recent developments in continuum percolation. Philosophical Magazine B 56 (6), pp. 991 – 1003. Cited by: §1, §3.1.
•  M. Bause, F.A. Radu, and U. Köcher (2017) Space-time finite element approximation of the Biot poroelasticity system with iterative coupling. Comput. Methods Appl. Mech. Engrg. 320, pp. 745 – 768. Cited by: §1.
•  L. Berger, R. Bordas, K. Burrowes, V. Grau, S. Tavener, and D. Kay (2016) A poroelastic model coupled to a fluid network with applications in lung modelling. Int. J. Numer. Meth. Biomed. Engng. 32 (1). Cited by: §1.
•  B. Berkowitz and R. P. Ewing (1998) Percolation theory and network modeling applications in soil physics. Surveys in Geophysics 19, pp. 23 – 72. Cited by: §1, §1, §1.
•  Y. Bernabe, W. F. Brace, and B. Evans (1982) Permeability, porosity and pore geometry of hot-pressed calcite. Mechanics of Materials 1 (3), pp. 173 – 183. Cited by: §1.
•  M. A. Biot (1941) General Theory of Three-Dimensional Consolidation. J. Appl. Phys. 12, pp. 155 – 164. Cited by: §1, §2.
•  D. Boffi, F. Brezzi, and M. Fortin (2013) Mixed Finite Element Methods and Applications. Vol. 44, Springer. Cited by: §5.1.
•  J. W. Both, M. Borregales, J. M. Nordbotten, K. Kumar, and F. A. Radu (2017) Robust fixed stress splitting for Biot’s equations in heterogeneous media. Appl. Math. Lett. 68, pp. 101 – 108. Cited by: §1.
•  D. Braess (2007) Finite elements: Theory, fast solvers, and applications in solid mechanics. Cambridge University Press. Cited by: §5.1.
•  S. C. Brenner and L. R. Scott (2007) The Mathematical Theory of Finite Element Methods. Springer Science & Business Media. Cited by: §5.1.
•  S. R. Broadbent and J. M. Hammersley (1957) Percolation processes: I. Crystals and mazes. In Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 53, pp. 629 – 641. Cited by: §1.
•  P. C. Carman (1937) Fluid Flow through Granular Beds. Trans. Inst. Chem. Eng. 15, pp. 150 – 166. Cited by: §1.
•  A. H. -D. Cheng (2016) Poroelasticity. theory and applications of transport in porous media. Vol. 27, Springer. Cited by: §1.
•  A. Ern and J. Guermond (2013) Theory and Practice of Finite Elements. Vol. 159, Springer Science & Business Media. Cited by: §5.1.
•  F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich (2003) A finite difference analysis of Biot’s consolidation model. Appl. Numer. Math. 44 (4), pp. 487 – 506. Cited by: §1.
•  F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich (2006) Staggered grid discretizations for the quasi-static Biot’s consolidation problem. Appl. Numer. Math. 56 (6), pp. 888 – 898. Cited by: §1.
•  J. B. Haga, H. Osnes, and H. P. Langtangen (2012) On the causes of pressure oscillations in low-permeable and low-compressible porous media. Int. J. Numer. Anal. Meth. Geomech. 36 (12), pp. 1507 – 1522. Cited by: §1.
•  Q. Hong, J. Kraus, M. Lymbery, and F. Philo (2019) Conservative discretizations and parameter-robust preconditioners for biot and multiple-network flux-based poroelasticity models. Numer. Linear Algebra Appl. 26, pp. e2242. Cited by: §1.
•  Q. Hong, J. Kraus, M. Lymbery, and M. F. Wheeler (2019) Parameter-robust convergence analysis of fixed-stress split iterative method for multiple-permeability poroelasticity systems. arXiv preprint arXiv:1812.11809v2. Cited by: §1.
•  X. Hu, C. Rodrigo, F. J. Gaspar, and L. T. Zikatanov (2017) A nonconforming finite element method for the Biot’s consolidation model in poroelasticity. J. Comput. Appl. Math. 310, pp. 143 – 154. Cited by: §1.
•  C. E. Koltermann and S. M. Gorelick (1995) Fractional packing model for hydraulic conductivity derived from sediment mixtures. Water Resour. Res. 31 (12), pp. 3283 – 3297. Cited by: §1.
•  J. Korsawe and G. Starke (2005) A Least-Squares Mixed Finite Element Method for Biot’s Consolidation Problem in Porous Media. SIAM J. Numer. Anal. 43 (1), pp. 318 – 339. Cited by: §1.
•  J. Kozeny (1927) Uber Kapillare Leitung der Wasser in Boden. Royal Academy of Science, Vienna, Proc. Class I 136, pp. 271–306. Cited by: §1.
•  P. Luo, C. Rodrigo, F. J. Gaspar, and C. W. Oosterlee (2015) Multigrid method for nonlinear poroelasticity equations. Comput. Visual. Sci. 17, pp. 255 – 265. Cited by: §1.
•  G. Mavko and A. Nur (1997) The effect of a percolation threshold in the Kozeny-Carman relation. Geophysics 62 (5), pp. 1480 – 1482. Cited by: §1.
•  A. Mikelić and M. F. Wheeler (2013) Convergence of iterative coupling for coupled flow and geomechanics. Computational Geosciences 17 (3), pp. 455 – 461. Cited by: §1.
•  P. J. Phillips and M. F. Wheeler (2007) A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case. Comput. Geosci. 11 (2), pp. 131 – 144. Cited by: §1.
•  L. B. Porter, R. W. Ritzi, L. J. Mastera, D. F. Dominic, and B. Ghanbarian-Alavijeh (2013) The Kozeny-Carman Equation with a Percolation Threshold. Ground Water 51 (1), pp. 92 – 99. Cited by: §1.
•  M. Rahrah and F. Vermolen (2018) Monte Carlo Assessment of the Impact of Oscillatory and Pulsating Boundary Conditions on the Flow Through Porous Media. Transp. Porous Med. 123 (1), pp. 125 – 146. Cited by: §3, §6.
•  M. Rahrah and F. Vermolen (2019) Uncertainty Quantification in Injection and Soil Characteristics for Biot’s Poroelasticity Model. In European Conference on Numerical Mathematics and Advanced Applications ENUMATH 2017, pp. 645 – 652. Cited by: §1.
•  C. Rodrigo, F.J. Gaspar, X. Hu, and L.T. Zikatanov (2016) Stability and monotonicity for some discretizations of the Biot’s consolidation model. Comput. Methods Appl. Mech. Engrg. 298, pp. 183 – 204. Cited by: §1, §5.1.
•  A. Segal (2012) Finite element methods for the incompressible Navier-Stokes equations. DIAM, Delft. Cited by: §6.
•  R. E. Showalter (2000) Diffusion in Poro-Elastic Media. J. Math. Anal. Appl. 251, pp. 310 – 340. Cited by: §5.1.
•  M. Spiegelman (1993) Flow in deformable porous media. Part 1. Simple analysis. J. Fluid Mech. 247, pp. 17 – 38. Cited by: §1.
•  C. Stover and E. W. Weisstein Percolation threshold. Note: From MathWorld - A Wolfram Web ResourceAccessed on 04 February 2019 External Links: Link Cited by: §3.1.3.
•  K. H. Støverud, M. Darcis, R. Helmig, and S. M. Hassanizadeh (2012) Modeling Concentration Distribution and Deformation During Convection-Enhanced Drug Delivery into Brain Tissue. Transp. Porous Med. 92 (1), pp. 119 – 143. Cited by: §1.
•  M. Tchonkova, J. Peters, and S. Sture (2008) A new mixed finite element method for poro-elasticity. Int. J. Numer. Anal. Meth. Geomech. 32 (6), pp. 579 – 606. Cited by: §1.
•  T. Tsai, K. Chang, and L. Huang (2006) Body force effect on consolidation of porous elastic media due to pumping. J. Chin. Inst. Eng. 29 (1), pp. 75 – 82. Cited by: §3.
•  P. A. Vermeer and A. Verruijt (1981) An accuracy condition for consolidation by finite elements. Int. J. Numer. Anal. Meth. Geomech. 5 (1), pp. 1 – 14. Cited by: §1.
•  J. Wan (2003) Stabilized finite element methods for coupled geomechanics and multiphase flow. Ph.D. Thesis, Stanford University. Cited by: §1.
•  H. F. Wang (2000) Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press. Cited by: §2.
•  S. Wang and K. Hsu (2009) Dynamics of deformation and water flow in heterogeneous porous media and its impact on soil properties. Hydrol. Process. 23, pp. 3569 – 3582. Cited by: §3.
•  M. Wheeler, G. Xue, and I. Yotov (2014) Coupling multipoint flux mixed finite element methods with continuous Galerkin methods for poroelasticity. Comput. Geosci. 18, pp. 57 – 75. Cited by: §1.
•  P. Wong (1988) The Statistical Physics of Sedimentary Rock. Physics Today 41, pp. 24 – 32. Cited by: §1, §3.1.
•  A. Ženíšek (1984) The existence and uniqueness theorem in Biot’s consolidation theory. Aplikace Matematiky 29 (3), pp. 194 – 211. Cited by: §5.1.