A moving finite element framework for fast infiltration in nonlinear poroelastic media

Poroelasticity theory can be used to analyse the coupled interaction between fluid flow and porous media (matrix) deformation. The classical theory of linear poroelasticity captures this coupling by combining Terzaghi’s effective stress with a linear continuity equation. Linear poroelasticity is a good model for very small deformations; however, it becomes less accurate for moderate to large deformations. On the other hand, the theory of large-deformation poroelasticity combines Terzaghi’s effective stress with a nonlinear continuity equation. In this paper, we present a finite element solver for linear and nonlinear poroelasticity problems on triangular meshes based on the displacement-pressure two-field model. We then compare the predictions of linear poroelasticity with those of large-deformation poroelasticity in the context of a two-dimensional model problem where flow through elastic, saturated porous media, under applied mechanical oscillations, is considered. In addition, the impact of introducing a deformation-dependent permeability according to the Kozeny-Carman equation is explored. We computationally show that the errors in the displacement and pressure fields that are obtained using the linear poroelasticity are primarily due to the lack of the kinematic nonlinearity. Furthermore, the error in the pressure field is amplified by incorporating a constant permeability rather than a deformation-dependent permeability.


Introduction
A saturated porous medium is composed of a porous solid material, fully saturated by a viscous fluid, flowing through connected pores. In deformable porous materials such as soils, rocks and tissues, the flow of the pore fluid and the deformation of the solid matrix are tightly coupled to each other. Poromechanics involves fluid flow in porous media that can deform when subjected to external forces and to variations in pressure of the saturating fluid. Moreover, poromechanical deformations are poroelastic when they are controlled by the reversible storage and release of elastic energy. Menel  In the last few decades, the mechanics of porous media has been of great interest due to its potential application in many geological and biological systems across a wide range of scales such as civil engineering [7,13,20,21,32,40,43,45], energy and environmental technologies [14,17,26,37,39,51], material science [29] and biophysics [18,30], where poromechanics plays an important role in modelling bones and soft tissues [1,15,41]. In physical chemistry, poromechanical processes include mass and heat transfer [53]. Additionally, poromechanics has been studied intensely in geophysics, in the context of consolidation of aquifers [33,36] and in the context of enhanced oil or gas recovery [4,38,44]. Our current main motivation is related to the development of a comprehensive model for fluid injection into a monolayer of soil particles subjected to surface mechanical oscillations. Therein, we include dynamic effects, especially a time and space dependent porosity and permeability.
Due to the high complexity and the unknown geometry of porous media, a fully resolved model is nearly impossible to obtain. Classically, in the theory of linear poroelasticity [5,48], fluid flow is described by Darcy's law and fluid mass conservation, and matrix deformation is described by Terzaghi's effective stress and linear elasticity. This theory is originally emerged in soil mechanics with the work of Terzaghi [46], but its general statement was given by Biot [5,6] using an elastic formulation for the solid matrix and Darcy's law for the fluid flow. This approach, as formulated in Biot's consolidation model [5], is valid for infinitesimal deformations of the solid. However, it becomes increasingly inappropriate for moderate to large deformations. The well-known theory of large-deformation poroelasticity [28] combines Darcy's law with Terzaghi's effective stress and nonlinear elasticity in a rigorous kinematic framework, leading to a strongly nonlinear coupling between the pore structure and the fluid flow [18,37]. Another nonlinear poroelasticity model that takes large deformations into account is considered in [8]. In this model, the mechanical deformation follows the Saint Venant-Kirchhoff constitutive law for hyperelastic solid materials and the fluid compressibility in the fluid equation is assumed to be nonlinear. In the current paper, the fluid phase is assumed to be incompressible and a linear stressstrain constitutive law is considered.
It is difficult to obtain analytical solutions for poroelasticity problems. Therefore, solving these problems relies mainly on numerical methods. In addition, numerical methods are necessary to solve large deformation poromechanical problems since these problems are inherently nonlinear. Poroelasticity problems have been attracting attention from the scientific computing community [25,31,50] (and references therein). Some recent work can be found in [19,24,25,35,42]. The simplicity of the displacement-pressure two-field formulation [31,52] is attractive and hence pursued by this paper. In addition, a numerical model has been developed to solve the poroelasticity equations following the continuous Galerkin (CG) finite element method. When using the finite element method to solve the poroelastic equations, the main challenge is to ensure convergence of the method and to prevent numerical instabilities that often manifest themselves in the form of spurious oscillations in the pressure field. It has been proved that this problem is caused by the saddle point structure in the coupled equations resulting in a violation of the famous Ladyzhenskaya-Babuska-Brezzi (LBB) condition [34], thus highlighting the need for a stable combination of mixed finite elements [22] such as the popular Taylor-Hood element.
In this work, we propose finite element methods for the resolution of the governing equations both in the theory of linear poroelasticity as in the largedeformation poroelasticity. The fluid-mass balance equation is discretised in time by a backward Euler scheme. The resulting system of nonlinearly coupled equations is solved by a standard Picard iterative procedure, which is linearly convergent. In the literature, this system is also solved by Newton's method [8], which is quadratically convergent. The drawbacks of the Newton-Raphson method are that the method is only locally convergent and that the computation of derivatives is needed. Another valuable alternative to Picard's method is the L-scheme [8,9]. The L-method is robust and linearly convergent, and does not involve the computation of any derivatives. Moreover, the convergence rate does not depend on the mesh size. Only a relatively mild constraint on the time step size is required when the hydraulic conductivity is not taken constant [27]. The L-scheme contains a constant parameter L > 0 which mimics the Jacobian from Newton iteration. However, in order to determine the parameter L, for any given problem, it is necessary to use apriorily derived convergence estimates [12]. In the current paper, Picard's iterative method is used since it is easy to understand and to implement and since it does not involve the computation of derivatives or unknown parameters. Furthermore, monolithic approaches for solving the quasi-static two-field poroelasticity equations are adopted. Another approach that is widely used in coupling the flow and the mechanics in porous media is the fixedstress split method [10], which falls within the class of segregated approaches. This method can be combined with a linearisation scheme for the nonlinear poroelasticity models (see [9] where the L-scheme is used for the linearisation).
To assess the accuracy of aforementioned constitutive laws and the performance of finite element methods presented herein, we consider a two-dimensional simulation. The obtained numerical results from both linear and nonlinear poroelasticity theories are compared with each other. The research problem we address in the present paper is for which applied mechanical oscillations it is sufficient to solve the linear poroelasticity model, which is computationally cheaper and simpler to solve.
The rest of this paper is organised as follows. In Section 2 the considered constitutive equations are summarised, including the employed permeability model. Section 3 presents a two-dimensional numerical example that is used to demonstrate the difference between the linear and nonlinear poroelasticity models. In Section 4, the nonlinear equations are discretised and the finite elements are combined with the first-order implicit Euler temporal discretisation to establish a solver for linear and nonlinear poroelasticities on triangular meshes, which couples the solid displacement and fluid pressure in a monolithic system. Furthermore, the nonlinear equation systems are solved using a Picard iterative method. Section 5 discusses the numerical results and some concluding remarks are reported in Section 6.

Governing equations
In the following, we briefly recall the equations governing the problem of a porous material subjected to oscillating mechanical deformations characterised by displacements u of the solid skeleton.
We consider a two-phase mixture composed of an elastic solid matrix whose voids are continuous and completely saturated by an incompressible Newtonian fluid. In this study, it is further assumed that the porous material is in an initial state of hydraulic and mechanical equilibrium, gravitational body force remains constant and the matrix grains are incompressible. Let t ⊂ R 3 denote a bounded domain occupied by a homogeneous and isotropic elastic body with boundary t and x = (x, y, z) ∈ t . Denote by 0 the reference domain corresponding to the poroelastic medium in the initial state and t the deformed domain. Furthermore, t denotes time, belonging to a half-open time interval I = (0, T ], with T > 0. To determine the local displacement of the skeleton of a porous medium as well as the fluid flow through the pores, the poroelastic equations with single-phase flow can be expressed as [28]: where σ and v f are defined by the following equations Stress-strain constitutive law: σ = λtr(ε)I + 2με; (2) Darcy's law: where σ is Terzaghi's effective stress tensor for the porous medium, p is the fluid pressure, ρ is the fluid density, g is the gravitational acceleration vector, u is the solid displacement vector, v f is Darcy's velocity, λ and μ are the Lamé constants, ε is the effective strain tensor, κ is the permeability of the porous medium and η is the fluid viscosity. Note that in Eq. 1b we have used the material derivative, which reads as: where v s = ∂u ∂t 0 is the solid velocity. This system needs to be complemented by appropriate boundary and initial conditions that will be specified in Section 3.
In the infinitesimal deformation range, corresponding to the assumptions that u 1 and ∂u/∂x 1, the model provided by Biot's theory of linear poroelasticity [5] is used: In the finite deformation range, the deformations are not very small and can not be neglected. Hence, the poroelastic equations with single-phase flow are expressed as: Furthermore, for the solid skeleton, we consider linear and nonlinear constitutive laws for the relationship between strain and displacement. Assuming that the solid deforms elastically, these relationships are quasi-static and reversible. Hencky elasticity is a nonlinear hyperelastic model that is based on a logarithmic strain measure and provides good agreement for the elastic behaviour of a wide variety of materials under moderate to large deformations [2,23]. The Hencky strain tensor can be written as [28]: where F = (I − ∇u) −1 is the deformation gradient tensor, with I denoting the identity tensor. The natural logarithm ln(·) is computed in each element of the tensor FF T . On the other hand, assuming that the porous material is linearly elastic, the linear strain tensor can be defined by: Considering Hencky elasticity law (7) combined with nonlinear poroelasticity (6) is the most appropriate model for moderate to large deformations. However, this model is computationally expensive. According to Auton and MacMinn [3], using nonlinear poroelasticity (6) combined with linear elasticity (8) offers a good compromise between accuracy, robustness and computational efficiency, demonstrating the same qualitative behaviour as the fully nonlinear model. Hence, in this paper, we adopt linear elasticity (8) for all models: ε = ε L . In addition, we assume that the solid and fluid phases are individually incompressible, such that deformation occurs only through rearrangement of the solid skeleton with corresponding changes in the local porosity. This is then likely to alter the permeability of the material. The deformation-dependent permeability can be determined using the Kozeny-Carman equation [49]: where d s is the mean grain size of the soil and the porosity θ is computed from the displacement vector using the porosity-dilatation relation (see [35,47]): with θ 0 the initial uniform porosity. In this study, the term porosity refers to the entire connected porosity. Since in linear poroelasticity it is assumed that the deformations are infinitesimal, this model is in the literature often combined with constant permeability. In this study, we consider three models: linear poroelasticity, where Eq. 5 is used, combined with constant permeability κ(x, t) = κ(x, 0) which we abbreviate as (L C ), linear poroelasticity combined with the Kozeny-Carman equation (L KC ) and nonlinear poroelasticity, where Eq. 6 is used, combined with the Kozeny-Carman equation (N KC ).

Numerical experiment
This section presents a numerical example that verifies the proposed finite element formulation and highlights the differences between the infinitesimal deformation and the finite deformation regimes in a poromechanical problem. In this example, we consider the effect of applied mechanical oscillations with small and large amplitudes.
In this problem, a poroelastic medium is instantaneously subjected to uniform boundary pressure on the left boundary. As soon as the boundary pressure is applied, excess pore pressure develops inside the domain, and so, the pore fluid starts to drain through the right boundary. The boundary pressure is maintained constant throughout. Figure 1 illustrates the setup of the problem. We consider a rectangular domain with initial width L and initial height H . At the top of the domain, an oscillating mechanical deformation is applied. No-flow conditions are imposed on the top and bottom boundaries. The material is assumed to be fully saturated and free of gravitational forces throughout. The boundary conditions for this problem are as follows: where t is the unit tangent vector at the boundary, n the outward unit normal vector and p pump is a prescribed pump pressure. Figure 1 shows the definition of the boundary segments. Initially, the following condition is fulfilled: For the boundary displacement u vib , a standing wave is considered, represented by: with γ the amplitude of the oscillation and t the time increment.

Numerical procedure
In this section, we outline the numerical procedures used to discretise the poroelastic models presented in Section 2 and to solve the resulting coupled fluid/solid finitedimensional problem. The weak form of the governing equations will be derived and discretised using a continuous Galerkin finite element approach with displacements and fluid pressures as primary variables. The suitability of the proposed methodology to model flow through elastic, saturated porous media under finite deformations will be demonstrated using the illustrative numerical example described in the previous section.

Weak formulation
We present a finite element framework for Eqs. 5 and 6, using the continuously deforming domain t with initially p pump and q| 4 = 0} and W = {w ∈ (H 1 ( t )) 2 : w| 1 = (0, u vib ) T and (w · n)| 3 = 0}. Furthermore, we consider the bilinear forms [35]: Using the notationu = ∂u/∂t, the variational formulation for Eq. 5 with boundary and initial conditions (11)- (12) consists of the following: Find (u(x, t), p(x, t)) ∈ (W ×Q) such that: with the initial condition u(0) = 0, and where In this work, we use the implicit Euler method for the discretisation in time, which is unconditionally stable and commonly used in computational poromechanics literature while for m = 0: u 0 = 0. At each time step, we solve the equations as a fully coupled system. Now we derive the variational formulation for Eq. 6b, first we multiply this equation by a basis function q ∈ Q 0 , and integrate the result over t , to get: First, we introduce J = [∇u − (∇ · u)I]v s + v f , then using Clairaut's theorem on equality of mixed partials for C 2 -functions over time and position, we obtain The definition of the material derivative (4) gives According to Dziuk and Elliott [16], it holds that Dq Dt = 0 for the Lagrangian basis functions that we will use in this study. Hence, we get using the definition of the material derivative (4): The divergence theorem gives: Furthermore, from Reynolds' theorem, it follows: This results into: Applying the divergence theorem again yields: From boundary conditions (11) it follows, using the definition of J, that: After applying the implicit Euler method for the temporal discretisation of v s : it follows that: We thus obtain, for each time step, a nonlinear system to be solved using an iterative scheme. Nonlinear algebraic system Eq. 27 can be solved by Picard's iterative procedure, where the subscripts (·) k−1 and (·) k denote the values of the previous and the current iterations respectively. In addition, we choose the initial guess u m 0 = u m−1 and the stopping criterion u k −u k−1 u k Picard's iterative scheme is also applied for solving the models that use the Kozeny-Carman equation. Thus, after having obtained the numerical approximation for the displacement in the previous iteration u k−1 , we update the porosity using Eq. 10. Subsequently, the Kozeny-Carman relation (9) is used to calculate the permeability in the current iteration k. Equation 14 and Eqs. 14a and 28 are solved by applying the finite element method, with triangular Taylor-Hood elements [11]. Regarding spatial discretisation, the displacement field is approximated using finite elements with quadratic basis functions, whereas continuous piecewise linear approximation is used for the pressure field. Time discretisation of the above dynamical equations is performed by using the implicit Euler method. Let P l h ⊂ H 1 ( t ) be a function space of piecewise polynomials on t of degree l. Hence, we define finite element approximations for W and Q as W l

Finite element formulation
. . , n u } and Q l h = Q ∩ P l h with basis {ψ j ∈ Q l h : j = 1, . . . , n p }, respectively. Subsequently, we introduce a spatial approximation for the functions u(x, t) and p(x, t), writing this in the form: in which the Dirichlet boundary conditions are imposed. When a continuously deforming grid is used, each trial function is time-dependent due to the motion of the grid. Hence, the finite-element trial solution is of the form: For the mesh points, it holds x(t) = x 0 + u(t), where x 0 are the coordinates of the reference domain 0 .

Numerical results
The Galerkin finite element method with triangular Taylor-Hood elements is employed for the solution of the discretised quasi-two-dimensional problems (5) and (6). The numerical investigations are carried out using the matrixbased software package MATLAB (version R2017a). The computational domain is a rectangular surface with initial width L = 1.0 m and initial height H = 1.0 m. The domain is discretised using a regular triangular grid, with x = y = 1/200. Mesh refinement did not yield any significant changes of the numerical solution. In addition, the hydraulic and mechanical properties used in the simulation can be found in Table 1. The solid material properties are characteristic of an unconsolidated, sandy formation.
Furthermore, the Lamé constants λ and μ in Eq. 2 are related to elasticity modulus E and Poisson's ratio ν by: λ = νE (1+ν)(1−2ν) and μ = E 2(1+ν) . The suitability of the proposed methodology to model flow through elastic porous media under infinitesimal and finite deformations is investigated in this study by means of the L 2 -norm of the computed displacements u and pressure field p . Subsequently, to compare the results from the different models, we compute the percentage change as follows: where A and B are two different models from the three models considered in this study: L C , L KC and N KC . In the generations of the simulation results, the time increment is chosen to be t = 0.1. In order to obtain some insight into the impact of the applied mechanical oscillations on the solid displacements and the fluid pressure, we present an overview of the simulation results in Figs. 2 and 3. In these simulations, water is injected into the soil at a constant pump pressure of 5 bar. We start with the simulation results for the nonlinear model (N KC ) without any oscillations applied, i.e. γ = 0. The simulated pressure and displacement profiles are provided in Fig. 2a and b. Mechanically, the deformations in the porous medium are negligible, other than a small  shift of the grains to the right, as a result of the force exerted on the grains by the injected water. As shown in Fig. 2b, the simulated pressure is almost linear. This means that the injected water flows in a horizontal direction through the domain from the left to the right boundary. In Fig. 2c and d, the numerical solutions at t = 0.9 are shown for the nonlinear model using applied oscillations with γ = 0.1. Figure 2c and d show the impact of the applied oscillations, imposed on the top of the domain, on the water flow. In contrast to the pressure shown in Fig. 2b, the numerical solution for the pressure in the problem with oscillations is no longer linear, but shows an oscillatory behaviour, as depicted in Fig. 2d. In this figure, we can see that the fluid pressure increases when the grains are pressed together by the applied oscillation. The norm of the simulated displacement and pressure profiles that have been obtained using an applied oscillation with γ = 0 (no oscillation) are depicted in Fig. 3a and b, while the simulated results that have been obtained using an applied oscillation with γ = 0.1 are provided in Fig. 3c and d.
In Fig. 3a and b, the behaviours of the displacement and pressure fields as function of the time, without applied oscillations, are shown. The difference between the linear models L C and L KC is negligibly small in the displacement field, while this difference is more visible in the pressure field for small times. This is a consequence of the different permeability relationships that are used in these models and that have more effect on the pressure field than on the solid deformations. Over time, the difference between the three models becomes smaller in the pressure field. However, the value of the norm of the displacement as a result of the nonlinear model is larger than the norm of the displacements in the linear models. Hence, for larger times, the nonlinearity has more effect on the displacement than on the pressure field. This is also expected from Eqs. 5 and 6. The impact of the applied mechanical oscillation is shown in Fig. 3c and d, where we notice an oscillatory behaviour in the displacement and pressure profiles. In these figures, we notice more similarity in the results of the three models for the displacement than for the fluid pressure. Furthermore, it is clear that the applied oscillation, which has an amplitude equal to 10% of the height of the domain, has a larger impact on the results than the adopted mathematical models.
The percentage change (31) in the norms of the simulated displacement u % and pressure p % is depicted in Figs. 4 and 5, for different values of the amplitude of the applied oscillations.
When comparing the linear models in Figs. 4 and 5, we notice that the impact of the Kozeny-Carman relation on the displacement is small, whereas this impact on the fluid pressure field becomes larger with increasing amplitude of the applied oscillation. The reason for this behaviour is that the permeability relationship is directly related to the pressure through Eq. 3. Moreover, the influence of the deformations on the porosity (see Eq. 10), and thus on the Kozeny-Carman permeability, is greater for larger deformations. In addition, the comparison between the linear and the nonlinear models both combined with Kozeny-Carman equation leads to the conclusion that the nonlinearity has a larger impact on the displacement field more than the permeability relationship. In contrast, the pressure field is more influenced by the permeability relation that is used than the nonlinearity in the models. In the extrema, the percentage change between the nonlinear model and the linear model with constant permeability is as large as the percentage change between the linear models and the percentage change between the nonlinear model and the linear model combined with the Kozeny-Carman relation.
Picard's iterative scheme is used to solve the nonlinear poroelasticity model (N KC ). For the previous values of the amplitude of the applied oscillations, the number of iterations per time step as function of the time t is depicted in Fig. 6. It can be seen from Fig. 6 that the number of Picard iterations stabilises with time and depends mildly on the magnitude of γ . The aim of this study was to quantify the amplitudes for which the linear poroelasticity model is accurate enough and demonstrates the same qualitative behaviour as the nonlinear poroelasticity model. As expected, for small applied mechanical oscillations, the difference between the linear and the nonlinear models is small both in the displacement as in the pressure fields. In Table 2, the upper limits of the absolute values of the percentage change for the different models are depicted.
For instance, given an applied oscillation with an amplitude γ = 0.08 (8% of the height of the domain), an accuracy of 10% can be obtained in both the displacement and the pressure fields using the linear poroelasticity model combined with Kozeny-Carman equation. However, this accuracy can only be obtained using the linear model combined with constant permeability for applied oscillations with an amplitude of 0.04 or smaller. For applied mechanical oscillations with amplitudes larger than γ = 0.12, the nonlinear poroelasticity model becomes unstable resulting in negative porosities. Hence, we will not consider these oscillations in this paper.

Discussion and conclusions
In this work, we have developed a mathematical model for fluid injection into a monolayer of soil particles subjected to surface mechanical oscillations, based on the two-field model (solid displacement and fluid pressure). This two-field mixed formulation is employed to calculate the solid displacement and the fluid pressure directly, in a monolithic system. In addition, we have included dynamic effects, such as a time-and space-dependent porosity and permeability. Firstly, we have summarised the governing equations both in the theory of linear poroelasticity as in the large-deformation poroelasticity. Subsequently, we have presented a finite element solver for the linear and nonlinear poroelasticity models, combined with constant and deformation-dependent permeability. This solver is developed on a triangular mesh and relies on quadratic basis functions for the discretisation of the displacement in elasticity and the continuous piecewise linear approximation for discretisation of pressure in Darcy's flow. These spatial discretisations are combined with the backward Euler temporal discretisation of the fluidmass balance equation. For the nonlinear poroelasticity equation, a weak formulation based on the motion of the solid was first presented, then linearisation of the resulting nonlinear coupled equation systems has been made using a standard Picard iterative procedure, which is subsequently implemented in a finite element code that is based on Taylor-Hood elements.
The suitability of the proposed methodology to model flow through elastic, saturated porous media under finite deformations is demonstrated using an illustrative numerical example. In this example, injection of a fluid into a two-dimensional fully saturated porous medium is considered, assuming that the solid material is subjected to surface mechanical oscillations with different amplitude sizes and that the fluid and solid constituents are individually incompressible.
Linear poroelasticity is a good model for very small deformations, besides that it is a simple model to solve and is computationally cheap. On the other hand, the wellknown large-deformation theory is more suitable to solve poroelasticity problems with moderate to large deformations. However, adopting this nonlinear mathematical model increases the computational complexity and cost, especially because the basis functions in the finite element code have to be updated in every iteration within the time integration. For this reason, the two-dimensional numerical example is used in this study to investigate the accuracy of the linear poroelasticity model for applied mechanical oscillations with different sizes. In addition, the impact of introducing a deformation-dependent permeability according to the Kozeny-Carman equation is explored.
The numerical example is solved using mechanical oscillations with amplitudes in the range Since the nonlinear poroelasticity model was unstable for applied mechanical oscillations with amplitudes larger than 0.12, resulting in negative porosities, we did not do any simulations for larger oscillations. In order to remove the nonphysical behaviour, one could analyse the improvement of the initial condition for the Picard iteration scheme of the displacement by the use of the linear model. As an alternative, one could investigate the performance of different time-integration schemes and time stepping. Another aspect is that, in the current work, we first applied time integration, followed by the finite element discretisation and finally the Picard method was implemented. One could reverse the order between the application of the finite element discretisation and Picard's method and investigate whether this reversal gives any improvement of the results. On the other hand, we could also analyse whether the model would actually predict these negative porosities and look whether the model (in terms of the Kozeny-Carman relation Eq. 9 and in the wake of Eq. 10) is entirely appropriate for this regime.
Firstly, by considering the numerical example without applied oscillations, we have shown that the nonlinearity has more effect on the displacement than on the pressure field. This is a consequence of the nonlinearity in the displacement that can be seen in Eq. 6. On the contrary, the impact of including a deformation-dependent permeability was larger on the pressure field than on the solid deformations. The reason for this behaviour is that the permeability relationship is directly related to the pressure through Darcy's law (3). Secondly, the impact of the applied mechanical oscillation was investigated by applying standing waves on the top surface of the solid matrix. From the numerical results, we noticed that the oscillatory behaviour was visible in the displacement and pressure profiles. Moreover, the differences between the three models (L C ), (L KC ) and (N KC ) are small for small applied oscillations, while these differences become larger by an increasing amplitude of the applied oscillation. Hence, the errors in the simulated displacement and pressure as a result of solving the linear poroelasticity model in the finite-deformation range increase when we choose applied oscillations with large amplitudes. The difference between the linear models can be explained by the impact of the large deformations on the porosity, which in turn has a larger impact on the Kozeny-Carman permeability. While this influence is not taken into account in the linear model combined with constant permeability.
In summary, for the studied problem and the set of parameters chosen, the use of the linear poroelasticity model for solving physical problems with finite deformations results in errors in the displacement and pressure fields that are mainly the consequence of the lack of the kinematic nonlinearity. To reduce these errors, especially in the pressure field, the linear poroelasticity model can preferably be combined with a deformation-dependent permeability, such as the Kozeny-Carman relationship.

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