1 Introduction

Fluid-saturated, deformable, porous materials are ubiquitous in many environmental, industrial and biological applications. Examples of such materials are soil, polymers, foams and living organs like muscle, liver, cartilage or even brain tissues. These types of materials can be characterized as porous biphasic materials with additional substances. The phases are understood to be the immiscible components, such as the porous solid matrix and the saturating fluid. However, the solutes are not considered further in this work. We reduce ourselves to the two-phase model, with the aim of demonstrating the performance of the developed model reduction method. Mathematical models are the classical tool to describe and predict the mechanics in these biphasic materials. They are derived from different multiphase homogenization approaches. The most common approaches are the Biot theory [6], the mixture theory [22] and the theory of porous media (TPM) [9]. Following the latter approach, this paper considers a biphasic model of coupled quasi-static balance equations of mass and momentum. We refer to [10] for a comprehensive introduction to the TPM.

For the aim of a simple presentation, we consider a quasi-static biphasic model with negligible constituents’ acceleration. We also assume isothermal conditions, negligible viscous extra stresses of the fluid and negligible gravity forces. Despite this choice of a simplified biphasic model, it is still a strongly coupled system of differential equations for fluid pressure and solid displacement. Hence, solving this model is of high computational complexity, mainly for nondeterministic models [29] or in data-driven evaluations [20], where the model has to be repeatedly solved. Therefore, it is our goal in this paper to reduce the complexity of the full biphasic model, in particular for any kind of limit regimes.

In the last decades, several model order reduction methods have been established. These are typically applied to the discretized models, with the advantage of preserving the details of the original continuous models. We refer here to the most popular projection-based model order reduction approaches, namely the singular value decomposition (SVD)-based methods. Examples of the SVD-based reduced methods are the method of balanced truncation [21] and the method of proper orthogonal decomposition (POD) [24] for biomechanical applications. The first method aims to approximate a system by identifying the important system modes and neglecting those of low impact, and the latter method approximates a given data set with a reduced basis in a low-dimensional subspace by applying a Galerkin projection. Further modifications of the POD method have been suggested to treat nonlinear terms in models. Examples of such extensions are: the Gappy-POD method [11], the discrete empirical interpolation method (DEIM) [8] and for biomechanical applications [13]. In addition to this, is the approach of hyper-reduction (HR) [27, 28]. We also refer to further developing modified approaches with a weighted inner product to perform a “goal-oriented” POD [7, 17]. For a general overview of these reductions methods, we refer to [12].

We focus in this paper on knowledge-driven model order reduction methods for saturated materials in thin domains, where the domain’s width and length are scale separated. In such domains, fluid flow can be assumed essentially hydrostatic in the transverse direction and the dynamics in the domain are almost in the longitudinal direction only. In petroleum application and CO2 sequestration, this assumption is known as the vertical equilibrium assumption [14, 19, 30], while in the field of hydrogeology, it is called the Dupuit’s assumption [16]. It has been utilized to reduce the computational complexity of two-phase flow models in nondeformable porous materials in thin domains. We refer here to the asymptotic approach proposed in [30] for Darcy regimes and extended in [4] for Brinkman regimes. This approach aims to reduce the complexity of the models by means of asymptotic analysis based on the geometrical width–length ratio of the domain. It has been tested for accuracy and computational efficiency, see, e.g., [2, 4]. In Brinkman regimes, the well-posedness of the resulting reduced model is investigated in [3]. Moreover, it is proved in [1] that the reduced model is the analytical limit of the full two-phase flow model as the geometrical ratio approaches zero. For other approaches that also utilize the hydrostatic assumption, we refer to [5, 14, 15, 23]. We also refer to [25, 26] for coupled two-scale PDE–ODE approaches and to [20] for model order reduction approaches in biomechanical applications.

The main contribution of this paper is a mathematical derivation of a reduced biphasic model consisting of two sub-parts. The first part is a limit model that describes the mechanics in the longitudinal direction of the domain. It results from the asymptotic analysis by letting the geometrical parameter of domain’s width–length ratio approach zero. It consists of three coupled one-dimensional equations for fluid pressure and material displacement in the longitudinal direction, in addition to displacement in the transverse direction at the boundary where load is applied. The second part is a corrector model (of small order) that re-captures the mechanics in the transverse direction. It consists of two full-dimensional equations for material displacements in both transverse and longitudinal directions. A further contribution is a numerical validation of the reduced model based on a comparison study with the full model for accuracy and computational efficiency. For the comparison, four scenarios with different input parameters, like permeability and load position, are designed.

The paper has the following structure. Section 2 introduces the original biphasic model (which we call here the full model) together with a set of suitable initial and boundary conditions. In Sect. 3, we derive the reduced model by, first, rescaling the full model into a dimensionless one that explicitly depends on the geometrical parameter of domain’s width–length ratio. Then, we apply formal asymptotic analysis to the dimensionless full model as the geometrical parameter approaches zero. Finally, we propose in Sect. 4 a numerical scheme, based on the finite difference method, to the reduced model. The reduced model is then numerically tested for accuracy and computational efficiency in comparison with the dimensionless full model.

2 Full biphasic model

We consider a fluid-saturated porous material \(\varphi \,=\,\varphi ^S\cup \varphi ^F\) in a thin rectangular domain \(\Omega =(0,H)\times (0,L)\) that undergoes small elastic deformation. Without loss of generality, and for more illustration, we consider an initial boundary value problem of a classical consolidation scenario, where the longitudinal direction is the z-direction and the transverse direction is the x-direction, see Fig. 1. The fluid \(\varphi ^{F}\) and solid \(\varphi ^{S}\) constituents are material incompressible. Thus, the true densities \(\rho ^{\alpha R}\) of the partial densities \(\rho ^{\alpha }\) with

$$\begin{aligned} \rho ^{\alpha }= n^{\alpha }\,\rho ^{\alpha R} \end{aligned}$$
(1)

remain constant, where \(n^{\alpha }\) defines the volume fraction of the phase \(\varphi ^{\alpha }\) with \(\alpha \in \{S,F\}\). The full-saturation condition is expressed with the sum of volume fractions via

$$\begin{aligned} n^S+n^F=1. \end{aligned}$$
(2)

As mentioned before, we assume isothermal conditions, negligible extra stresses of the low-viscous fluid, in addition to negligible gravity forces and constituents’ acceleration terms. Consequently, the biphasic model is quasi-static and governed by the balance equation of momentum

$$\begin{aligned} {\mathbf{0}} \,=\, {\mathrm{div}}\,{\mathbf{T}}\,=\, {\mathrm{div}}\,({\mathbf{T}}^S_E\,-\,{\lambda }\,{\mathbf{I}})\,, \end{aligned}$$
(3)

the balance equation of volume

$$\begin{aligned} 0\,=\,{\mathrm{div}}\,({\mathbf{u}}^S)^\prime _S\,+\,{\mathrm{div}}\,(n^{F}{\mathbf{w}}^F) \end{aligned}$$
(4)

and Darcy’s law

$$\begin{aligned} n^F{\mathbf{w}}^F\,=\,-\,\displaystyle \frac{k^{SF}}{\mu ^{FR}}\,{\mathrm{grad}}\,\lambda . \end{aligned}$$
(5)
Fig. 1
figure 1

A sketch of a thin rectangular domain \(\Omega =(0,H)\times (0,L)\) with load applied at the drained top boundary. The left, right and bottom boundaries are undrained and fixed in their perpendicular direction

Full size image

In Eq. (3), \({\mathbf{T}}\) is the Cauchy stress tensor, \({\mathbf{T}}^S_E\) is the solid’s extra stress and \(\lambda \) denotes the fluid’s pore pressure. A further simplification of the problem follows by using a materially and geometrically linear model. Hence, the Cauchy extra stress simplifies to

$$\begin{aligned} {\mathbf{T}}^{ S}_{E}\,=\,2\,\mu ^S{{\varvec{\varepsilon }}^S}\,+\,\lambda ^S({{\varvec{\varepsilon }}^S}\,\cdot \,{\mathbf{I}})\,{\mathbf{I}} \end{aligned}$$
(6)

where \({{\varvec{\varepsilon }}^S}\) is the linearized strain with

$$\begin{aligned} {{\varvec{\varepsilon }}^S}\,=\,\frac{1}{2}\,(\text {grad }\,{\mathbf{u}}^S\,+\,\text {grad }^T{\mathbf{u}}^S), \end{aligned}$$
(7)

and \(\mu ^S\) and \(\lambda ^S\) are the Lamé constants. Consequently, Eq. (3) can be expressed by the solid displacement field \({\mathbf{u}}^S\). In Eqs. (4) and (5), \({\mathbf{w}}^F\) is the seepage velocity, \(k^{SF}\) the intrinsic (isotropic) permeability and \(\mu ^{FR}\) the fluid’s viscosity.

The primary variables in the biphasic model are the fluid’s pressure \(\lambda \) and the solid displacement \({\mathbf{u}}^S=(u_1^S,u_2^S)^T\), where \(u_1^S\) and \(u_2^S\) are the displacement components in the horizontal (transverse) and vertical (longitudinal) directions, respectively. We set \({\mathbf{w}}^F = (w_1^F, w_2^F)\), where \(w_1^F\) and \(w_2^F\) are the seepage velocity components in the horizontal and vertical direction, respectively. Then, the volume balance equation (4) can be written as

$$\begin{aligned} \partial _x (n^F w_1^F) + \gamma ^{-2}\partial _x (u_1^S)' + \partial _z (n^F w_2^F) + \partial _z (u_2^S)'&= 0,\nonumber \\ n^F w_1^F&= -\frac{k^{SF}}{\mu ^{FR}} \partial _x \lambda ,\nonumber \\ n^F w_2^F&= -\frac{k^{SF}}{\mu ^{FR}} \partial _z \lambda . \end{aligned}$$
(8)

Note that we scaled second term in the first equation in (8) by \(\gamma ^{-2}\), where \(\gamma =\tfrac{H}{L}\) is the domain’s width–length ratio. This approach allows deriving balance equation for the horizontal processes in the system. It implies that the gradient of the solid velocity in this direction is large enough and will not lead to a trivial horizontal displacement as \(\gamma \) approaches zero. Similar approaches have been repeatedly used in the literature, and we refer to [18], where Darcy’s law is derived from the Navier–Stokes equation via homogenization. Substituting the explicit form of the strain

$$\begin{aligned} {{\varvec{\varepsilon }}^S}=\left( \begin{array}{c c} \partial _x u_1^S &{} 0.5(\partial _x u_2^S+\partial _z u_1^S)\\ 0.5(\partial _z u_1^S + \partial _x u_2^S) &{} \partial _z u_2^S \end{array}\right) \end{aligned}$$

into the momentum balance equation (3) produces

$$\begin{aligned} -\partial _x \lambda + \lambda ^S (\partial _{xx} u_1^S + \partial _{xz}u_2^S)+ \mu ^S(2 \partial _{xx}u_1^S + \partial _{zz}u_1^S+\partial _{xz}u_2^S)=0, \end{aligned}$$
(9)

and

$$\begin{aligned} -\partial _z \lambda + \lambda ^S (\partial _{zx} u_1^S + \partial _{zz}u_2^S)+ \mu ^S( \partial _{xx}u_2^S + \partial _{xz}u_1^S+2\partial _{zz}u_2^S)=0. \end{aligned}$$
(10)

Throughout the paper, we call Eqs. (8), (9) and (10) the full model. It is noticeable for the consolidation problem in Fig. 1 that the fluid is essentially hydrostatic in the horizontal (transverse) direction and flows almost only in the vertical (longitudinal) direction. The right, bottom and left sides of the domain are undrained, and fixed in the direction perpendicular to the boundary. A load function \(q=q(x,t)\) is applied to the top side. The Dirichlet and Neumann boundary conditions are summarized as

$$\begin{aligned} 2\mu ^S {{\varvec{\varepsilon }}^S} + \lambda ^S({{\varvec{\varepsilon }}^S} \cdot I)I - \lambda I&= q \quad \text { on } \partial \Omega _{top}\times (0,T),\nonumber \\ k^{SF}\nabla \lambda \cdot \mathbf {n}&=0 \quad \text { on } \partial \Omega _{bot}\cup \partial \Omega _{left}\cup \partial \Omega _{right}\times (0,T),\nonumber \\ u_1^S&=0 \quad \text { on } \partial \Omega _{left}\cup \partial \Omega _{right}\times (0,T),\nonumber \\ u_2^S&=0 \quad \text { on } \partial \Omega _{bot}\times (0,T), \end{aligned}$$
(11)

where the vector \({\mathbf{n}}\) in (11)\(_2\) is the outer normal to the boundary. We also set periodic boundary conditions on the left and right boundaries

$$\begin{aligned} {\mathbf{u}}^S(0,z,t)&={\mathbf{u}}^S(H,z,t) \quad \text { for } z\in (0,L),\,t\in (0,T),\nonumber \\ \nabla {\mathbf{u}}^S(0,z,t)&=\nabla {\mathbf{u}}^S(H,z,t) \quad \text { for } z\in (0,L),\,t\in (0,T). \end{aligned}$$
(12)

The initial conditions for this problem are defined as

$$\begin{aligned} {\mathbf{u}}^S&= {\mathbf{u}}^S_I \quad \text { in }\Omega ,\nonumber \\ \lambda&= \lambda _I \quad \text { in }\Omega . \end{aligned}$$
(13)

Remark 1

The choice of periodic boundary conditions in Eq. (12) fits with the classical consolidation problem under consideration. However, and as will be seen in the next sections, the derivation of the reduced model can be extended to more general boundary conditions.

3 Reduced biphasic model

In this section, we rescale the full model (8), (9) and (10) using dimensionless variables. This leads to a dimensionless full model that expresses the importance of the domain’s geometrical ratio \(\gamma >0\) explicitly. Then, by applying standard asymptotic analysis, based on the small \(\gamma \), to the dimensionless full model we derive a reduced biphasic model consisting of two main parts: The first is a limit model for \(\gamma \) approaches zeros, while the second is a corrector that re-captures the process in the horizontal direction.

3.1 Dimensionless model

We define the geometrical parameter \(\gamma =H/L\) and the dimensionless variables

$$\begin{aligned} \begin{array}{ccll} {\overline{x}}=\dfrac{x}{H},\quad \quad &{}{\overline{z}}=\dfrac{z}{L},\quad \quad &{} {\overline{t}}=\dfrac{t}{L/\psi }, \quad \quad &{}{\overline{K}}^{SF}=\dfrac{k^{SF}}{k_{m}^{SF}},\quad \\ {\overline{w}}_1^F =\dfrac{w_1^F}{\psi },\quad \quad &{} {\overline{w}}_2^F =\dfrac{w_2^F}{\psi },\quad \quad &{} {\overline{\lambda }}=\dfrac{\lambda }{q_m}, &{} \end{array} \end{aligned}$$
(14)

where \(\psi =\tfrac{q_m K_m^{SF}}{\mu ^{FR}L}\), \(k_m^{SF}=\frac{1}{|\Omega |}\int _{\Omega }k^{SF}\,dx\,dz\) is the mean value of the intrinsic permeability in the domain, and \(q_m=\frac{1}{H}\int _0^H |q|\,dx\) is the mean load applied to the top of the domain (see Eq. (11)\(_1\)). Substituting these variables into the volume balance equation (8), then omitting the par-signs yields

$$\begin{aligned} \gamma ^{-2}\partial _x (u_1^S)'+ \gamma ^{-1} \partial _x(n^F w_1^F)+ \partial _z (u_2^S)'+ \partial _z( n^F w_2^F )&=0,\nonumber \\ n^F w_1^F&= -\gamma ^{-1} K^{SF} \partial _x \lambda ,\nonumber \\ n^F w_2^F&= -K^{SF}\,\partial _z \lambda . \end{aligned}$$
(15)

Substituting the dimensionless variables (14) into momentum balance equations (9) and (10), then multiplying the resulting equations by H and L, respectively, produces

$$\begin{aligned} -\partial _x \lambda + {\hat{\lambda }}^S (\partial _{xx} u_1^S + \partial _{xz} u_2^S)+ {\hat{\mu }}^S(2 \partial _{xx}u_1^S + \gamma ^{2}\partial _{zz}u_1^S+\partial _{xz}u_2^S)=0. \end{aligned}$$
(16)

and

$$\begin{aligned} -\partial _z \lambda + {\hat{\lambda }}^S (\partial _{zx} u_1^S + \partial _{zz}u_2^S)+ {\hat{\mu }}^S( \gamma ^{-2} \partial _{xx}u_2^S + \partial _{xz}u_1^S+2\partial _{zz}u_2^S)=0. \end{aligned}$$
(17)

In these two equations, the par-signs are omitted and the parameters \({{\hat{\lambda }}}^S\) and \({{\hat{\mu }}}^S\) are the dimensionless Lamé constants defined as

$$\begin{aligned} {\hat{\lambda }}^S = \frac{\lambda ^S}{q_m},\quad \text { and }\quad {\hat{\mu }}^S = \frac{\mu ^S}{q_m}. \end{aligned}$$

The Dirichlet and Neumann boundary conditions for the dimensionless model are now summarized as follows

$$\begin{aligned} 2{{\hat{\mu }}}^S {{\varvec{\varepsilon }}^S} + {{\hat{\lambda }}}^S({{\varvec{\varepsilon }}^S} \cdot I)I - \lambda I&= \frac{q}{q_m} \quad \text { on } \partial \Omega _{top}\times (0,T),\nonumber \\ K^{SF}\nabla \lambda \cdot \mathbf {n}&=0 \quad \text { on } \partial \Omega _{bot}\cup \partial \Omega _{left}\cup \partial \Omega _{right}\times (0,T),\nonumber \\ u_1^S&=0 \quad \text { on } \partial \Omega _{left}\cup \partial \Omega _{right}\times (0,T),\nonumber \\ u_2^S&=0 \quad \text { on } \partial \Omega _{bot}\times (0,T), \end{aligned}$$
(18)

where \(\Omega =(0,1)\times (0,1)\) is the rescaled domain.

3.2 Asymptotic analysis

We consider a solution (\(u_1^{S,\gamma },u_2^{S,\gamma },\,w_1^{F,\gamma },\,w_2^{F,\gamma },\,\lambda ^{\gamma }\)) of the dimensionless full model (15)–(17) and define the asymptotic expansions

$$\begin{aligned} \begin{array}{cll} Z^{\gamma }&{}=Z_0+\gamma Z_1+{\mathcal {O}}(\gamma ^2)&{}\quad \quad Z\in \{u_1^S,\,u_2^S,\,w_2^F,\,\lambda \},\\ Z^{\gamma }&{}=\gamma Z_1+{\mathcal {O}}(\gamma ^2)&{}\quad \quad Z\in \{w_1^F\}. \end{array} \end{aligned}$$
(19)

Note that the choice that \(w_1^F={\mathcal {O}}(\gamma )\) corresponds to the assumption that the fluid is hydrostatic in the x-direction.

Mass Balance Equations: Substituting the expansions (19) into the horizontal component of Darcy equation (15)\(_2\) yields

$$\begin{aligned} n^F \big ( w_{1,1}^F + \gamma ^1 w_{1,2}^F \big )= - \gamma ^{-2} K^{SF} \partial _x(\lambda _0+\gamma \lambda _1+\gamma ^2 \lambda _2). \end{aligned}$$
(20)

Equating the terms of the highest power of \(\gamma \), namely \({\mathcal {O}}(\gamma ^{-2})\), yields

$$\begin{aligned} K^{SF}\,\partial _x \lambda _0= 0 . \end{aligned}$$
(21)

Assuming that the intrinsic permeability \(K^{SF}\) is strictly positive in the domain, Eq. (21) implies that \(\lambda _0\) is independent of the horizontal direction

$$\begin{aligned} \lambda _0 = \lambda _0(z,t). \end{aligned}$$
(22)

Similarly, equating the terms of order \({\mathcal {O}}(\gamma ^{-1})\) implies

$$\begin{aligned} \lambda _1=\lambda _1(z,t). \end{aligned}$$
(23)

The terms of order \({\mathcal {O}}(\gamma ^0)\) in (20) fulfill the equation

$$\begin{aligned} n^F w^F_{1,1}=-K^{SF} \partial _x \lambda _2. \end{aligned}$$
(24)

Now, we substitute the expansions (19) into the vertical component of Darcy equation (15)\(_3\) and consider the terms of order \({\mathcal {O}}(\gamma ^0)\). Then, we have

$$\begin{aligned} n^F w^F_{2,0}=-K^{SF} \partial _z \lambda _0. \end{aligned}$$
(25)

Finally, we substitute (19) into (15)\(_1\). This leads to

$$\begin{aligned}&\gamma ^{-2}\partial _x (u_{1,0}^S+\gamma u_{1,1}^S+\gamma ^2 u_{1,2}^S)'+ \partial _x(n^F ( w_{1,1}^F + \gamma w_{1,2}^F))\\&\quad + \, \partial _z (u_{2,0}^S+\gamma u_{2,1}^S+\gamma ^2 u_{2,2}^S)'+ \partial _z( n^F (w_{2,0}^F +\gamma w_{2,1}^F + \gamma ^2w_{2,2}^F))=0. \end{aligned}$$

Equating the terms of order \({\mathcal {O}}(\gamma ^{-2})\) gives \(\partial _x u_{1,0}^S=0\). Then, using the boundary condition (11)\(_3\), we obtain

$$\begin{aligned} u_{1,0}^S = 0 . \end{aligned}$$
(26)

Similarly, the terms of \({\mathcal {O}}(\gamma ^{-1})\) satisfy

$$\begin{aligned} u_{1,1}^S = 0 . \end{aligned}$$
(27)

However, equating the terms of order \({\mathcal {O}}(\gamma ^{0})\) gives

$$\begin{aligned} \partial _x ( u_{1,2}^S)' + \partial _x(n^F w_{1,1}^F )+ \partial _z (u_{2,0}^S)'+ \partial _z( n^F w_{2,0}^F )=0. \end{aligned}$$
(28)

Setting Eqs. (24) and (25) into (28) produces

$$\begin{aligned} \partial _x (u_{1,2}^S)'- \partial _x (K^{SF}\,\partial _x \lambda _2) +\partial _z (u_{2,0}^S)'- \partial _z (K^{SF}\, \partial _z \lambda _0)=0. \end{aligned}$$
(29)

Integrating this equation over the horizontal direction from 0 to 1 produces

$$\begin{aligned} \partial _z \int _0^1 (u_{2,0}^S)'\,dx - \partial _z \int _0^1 (K^{SF}\,\partial _z \lambda _0)\,dx = \Psi _{right} -\Psi _{left} - (u_{1,2}^S)'|_{x=0}^{x=1} , \end{aligned}$$

where \(\Psi _{right}= (K^{SF}\,\partial _x \lambda _2)|_{x=1}\) and \(\Psi _{left}= (K^{SF}\,\partial _x \lambda _2)|_{x=0}\) are the flux terms through the left and right boundaries. Then, using the boundary conditions (11)\(_2\) and (11)\(_3\), the previous equation reduces to

$$\begin{aligned} \partial _z \int _0^1 (u_{2,0}^S)'\,dx - \partial _z \int _0^1 (K^{SF}\,\partial _z \lambda _0)\,dx =0 . \end{aligned}$$
(30)

Applying the equilibrium result (22) for \(\lambda _0\) and (35) for \(u_{2,0}^S\) leads to the one-dimensional equation for \(\lambda _0\) and \(u_{2,0}^S\)

$$\begin{aligned} \partial _z (u_{2,0}^S)' - \partial _z (K_{av}^{SF}\,\partial _z \lambda _0) =0 , \end{aligned}$$
(31)

where \(K^{SF}_{av}= \int _0^1 K^{SF}\,dx\). Now, we integrate Eq. (29) again over the horizontal direction from 0 to x, which leads to

$$\begin{aligned} (u_{1,2}^S)'- K^{SF} \partial _x \lambda _2+ x \,\partial _z (u_{2,0}^S)' - \partial _z\big (\partial _z \lambda _0\int _0^x K^{SF} \,dx\big ) =-\Psi _{left} + (u_{1,2}^S)'|_{x=0}. \end{aligned}$$

Reordering the terms in the this equation leads to a new formula for the pressure gradient \(\partial _z \lambda _2\)

$$\begin{aligned} K^{SF} \partial _x \lambda _2 = x \,\partial _z (u_{2,0}^S)' + (u_{1,2}^S)' - \partial _z\big (\partial _z \lambda _0\int _0^x K^{SF} \,dx\big ) +\Psi _{left} - (u_{1,2}^S)'|_{x=0}. \end{aligned}$$

Then, using Eq. (22), the later coming result in Eq. (35), in addition to the boundary conditions (11)\(_2\) and (11)\(_3\), we obtain the new formula \(\partial _z \lambda _2\)

$$\begin{aligned} K^{SF} \partial _x \lambda _2 = x \,\partial _z (u_{2,0}^S)' + (u_{1,2}^S)' - \partial _z\big (\partial _z \lambda _0\int _0^x K^{SF} \,dx\big ) . \end{aligned}$$
(32)

Remark 2

Note that whenever the intrinsic permeability \(K^{SF}\) is independent of the horizontal direction x and using Eq. (31), Eq. (32) reduces to

$$\begin{aligned} K^{SF} \partial _x \lambda _2 = (u_{1,2}^S)'. \end{aligned}$$
(33)

Momentum Balance Equations: We substitute the asymptotic expansions (19) into the momentum balance equation (17). Then, we have

$$\begin{aligned}&- \partial _z (\lambda _0\,+\, \gamma \lambda _1+\gamma ^2 \lambda _2) + {\hat{\lambda }}^S \big (\partial _{zx} (u_{1,0}^S +\gamma u_{1,1}^S + \gamma ^2 u_{1,2}^S)\nonumber \\&\quad + \partial _{zz}(u_{2,0}^S+\gamma u_{2,1}^S + \gamma ^2 u_{2,2}^S)\big )+{\hat{\mu }}^S\big (\gamma ^{-2}\partial _{xx}(u_{2,0}^S+\gamma u_{2,1}^S+\gamma ^2 u_{2,2}^S) \nonumber \\&\quad + \partial _{xz}(u_{1,0}^S+\gamma u_{1,1}^S + \gamma ^2 u_{1,2}^S) + 2 \partial _{zz} (u_{2,0}^S + \gamma u_{2,1}^S+\gamma ^2 u_{2,2}^S)\big )=0. \end{aligned}$$
(34)

Equating the term of order \({\mathcal {O}}(\gamma ^{-2})\) leads to

$$\begin{aligned} \partial _{xx}u_{2,0}^S=0. \end{aligned}$$

This equation together with the periodic boundary condition (12) implies that

$$\begin{aligned} u_{2,0}^S=u_{2,0}^S(z,t). \end{aligned}$$
(35)

Similarly, equating the terms of order \({\mathcal {O}}(\gamma ^{-1})\) leads to

$$\begin{aligned} u_{2,1}^S=u_{2,1}^S(z,t). \end{aligned}$$
(36)

Equating the terms of order \({\mathcal {O}}(\gamma ^0)\) in (34), then using the result (26) leads to

$$\begin{aligned} -\partial _z \lambda _0 + ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{zz}u^S_{2,0} +({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{xz}u^S_{1,0} +{{\hat{\mu }}}^S \partial _{xx} u_{2,2}^S=0. \end{aligned}$$
(37)

Note that \(u^S_{1,0}=0\) inside the dimensionless domain \((0,1)\times (0,1)\) as shown in (26), but not at the top boundary as will be shown later. Integrating this equation over the x-direction gives

$$\begin{aligned} -\partial _z \lambda _0 + \Bigg ( \int _0^1({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,dx\Bigg ) \partial _{zz}u_{2,0}^S&= ~({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\big (\partial _{z}u^S_{1,0}|_{x=1} -\partial _x u^S_{1,0}|_{x=0}\big )\nonumber \\&\quad +{{\hat{\mu }}}^S \big (\partial _{x} u_{2,2}^S|_{x=1}-\partial _x u_{2,2}^S|_{x=0}\big ). \end{aligned}$$
(38)

Then, the periodic boundary condition (12) produces the one-dimensional equation

$$\begin{aligned} -\partial _z \lambda _0 + \Bigg ( \int _0^1({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,dx\Bigg ) \partial _{zz}u_{2,0}^S = 0. \end{aligned}$$
(39)

This equation together with (31) represents a coupled system of one-dimensional equations for the primary variables \(\lambda _0\) and \(u_{2,0}^S\).

Finally, we substitute the asymptotic expansions (19) into the dimensionless momentum balance equation (16). Then, we have

$$\begin{aligned}&- \partial _x (\lambda _0+\gamma \lambda _1 + \gamma ^2 \lambda _2) + {\hat{\lambda }}^S \big (\partial _{xx} (u_{1,0}^S + \gamma u_{1,1}^S+\gamma ^2 u_{1,2}^S) \\&\quad + ~\partial _{xz}(u_{2,0}^S+\gamma u_{2,1}^S + \gamma ^2 u_{2,2}^S)\big )+ {\hat{\mu }}^S\big (2 \partial _{xx}(u_{1,0}^S+\gamma u_{1,1}^S+\gamma ^2 u_{1,2}^S) \\&\quad + \gamma ^{2}\partial _{zz}(u_{1,0}^S+\gamma u_{1,1}^S )+ \partial _{xz}(u_{2,0}^S+\gamma u_{2,1}^S+\gamma ^2 u_{2,2}^S)\big )=0. \end{aligned}$$

Equating the terms of order \({\mathcal {O}}(\gamma ^{0})\), then using the result (31) produces

$$\begin{aligned} \big ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\big )\partial _{xx} u_{1,0}^S + \big ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\big )\partial _{xz} u_{2,0}^S =0. \end{aligned}$$
(40)

Note that this equation emphasizes that \(u_{1,0}^S=0\) (see Eq. (26)) in the dimensionless domain \((0,1)\times (0,1)\), as a consequence of (35). However, Eq. (40) includes also the mechanics of \(u_{1,0}^S\) at the upper boundary of the domain where the load \(q<0\) is applied. Hence, this equation can be written as a one-dimensional equation for \(u_{1,0}^S\) at the top boundary of the domain (\(z=1\))

$$\begin{aligned} \big ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\big )\partial _{xx} u_{1,0}^S|_{z=1} =- \big ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\big )\partial _{xz} u_{2,0}^S|_{z=1}. \end{aligned}$$
(41)

Similarly, the terms of order \({\mathcal {O}}(\gamma )\) satisfy

$$\begin{aligned} \big ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\big )\partial _{xx} u_{1,1}^S + \big ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\big )\partial _{xz} u_{2,0}^S =0. \end{aligned}$$
(42)

This equation, however, implies that \(u_{1,1} = 0\) inside the domain \((0,1)\times (0,1)\) and at the upper boundary too, because the load q is fulfilled by the zero component of the asymptotic expansions only. Finally, we equate the terms of order \({\mathcal {O}}(\gamma ^{2})\) and use (26). Then, we obtain

$$\begin{aligned} - \partial _x \lambda _2 + ({\hat{\lambda }}^S+2{{\hat{\mu }}}^S) \partial _{xx} u_{1,2}^S + ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S) \partial _{xz} u_{2,2}^S =0. \end{aligned}$$

Setting result (32) into this equation produces

$$\begin{aligned}&-( u_{1,2}^S)'+ ({\hat{\lambda }}^S+2{{\hat{\mu }}}^S)K^{SF} \partial _{xx} u_{1,2}^S + ({\hat{p}}^S + {{\hat{\mu }}}^S) K^{SF}\partial _{xz} u_{2,2}^S \nonumber \\&\qquad \qquad = ~x\,\partial _z (u_{2,0}^S)' - \partial _z\left( \partial _z p_0 \int _0^x K^{SF}\,dx\right) . \end{aligned}$$
(43)

In the following, we summarize the reduced biphasic model as a combination of two sub-models. The first is a limit model that describes the mechanics mainly in the vertical (longitudinal) direction. It consists of three one-dimensional equations for the components \(\Lambda =\lambda _0,\,U_2^S=u^S_{2,0}\) and \(U_1^S=u^S_{1,0}\). The second model is a corrector of the limit model and describes the mechanics in the horizontal (transverse) direction. It is a coupled system of two two-dimensional equations for the components \(u_2^S=u^S_{2,2}\) and \(u_1^S=u^S_{1,2}\). The reduced model can be written as

$$\begin{aligned} \text {reduced model} = \text {limit model } + \gamma ^2 \,\text {corrector model}. \end{aligned}$$
(44)

The limit model is a coupled system of the one-dimensional equations

$$\begin{aligned} \partial _z (U_2^S)' - \partial _z \big (K_{av}^{SF}\,\partial _z \Lambda \big )&=0 ,\nonumber \\ -\partial _z \Lambda + \Bigg ( \int _0^1({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,dx\Bigg ) \partial _{zz}U_2^S&=0,\nonumber \\ \bigg ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\bigg )\partial _{xx} U_1^S|_{z=1}&= - \bigg ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\bigg )\partial _{xz} U_2^S|_{z=1}, \end{aligned}$$
(45)

while the corrector model is a coupled system of two-dimensional equations

$$\begin{aligned}&-( u_1^S)'+ ({\hat{\lambda }}^S+2{{\hat{\mu }}}^S)K^{SF} \partial _{xx} u_1^S + ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S) K^{SF}\partial _{xz} u_2^S \nonumber \\&\qquad \quad \quad = x\,\partial _z (U_2^S)' - \partial _z\left( \partial _x \Lambda \int _0^x K^{SF}\,dx\right) ,\nonumber \\&\qquad \quad \qquad {{\hat{\mu }}}^S \partial _{xx} u_2^S = \partial _z \Lambda - ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S) \partial _{zz}U_2^S - ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{xz} U_1^S|_{z=1}. \end{aligned}$$
(46)

Remark 3

It is remarkable that the right side of the corrector model vanishes whenever the permeability \(K^{SF}\) and the material parameters \({{\hat{\lambda }}}^S\) and \({{\hat{\mu }}}^S\) are independent of the transverse direction x. In such cases, and assuming the initial and boundary condition in (18), the reduced model simplifies to the limit model only.

Remark 4

The derivation of the reduced model can be extended to more general boundary conditions with less periodicity constrains. As an example, we assume the following boundary conditions

$$\begin{aligned} 2\mu ^S {{\varvec{\varepsilon }}^S} + \lambda ^S({{\varvec{\varepsilon }}^S} \cdot I)I - \lambda I&= q \quad \text { on } \partial \Omega _{top}\times (0,T),\nonumber \\ k^{SF}\nabla \lambda \cdot \mathbf {n}&=\Psi \quad \text { on } \partial \Omega _{bot}\cup \partial \Omega _{left}\cup \partial \Omega _{right}\times (0,T), \end{aligned}$$
(47)

where \({\mathbf{n}}\) is the outer normal to the boundary. In addition, we set periodic boundary conditions on the first two components in the asymptotic expansions of the displacement vector, which fits to our consideration of almost unidirectional dynamics in thin domains. So, we assume that the first \({\mathbf{u}}^S_0=(u^S_{1,0},u^S_{2,0})^T\) and second \(\mathbf{{u}}^S_1=(u^S_{1,1},u^S_{2,1})^T\) terms in the asymptotic expansions of the displacement vector \({\mathbf{u}}^S = {\mathbf{u}}^S_0 + \gamma {\mathbf{u}}^S_1 + \gamma ^2{\mathbf{u}}^S_2\) are periodic in the horizontal direction, but not necessarily the third term \({\mathbf{u}}^S_2=(u^S_{1,2},u^S_{2,2})^T\) as it accounts to the processes in the horizontal direction. These can be summarized as follows

$$\begin{aligned} \begin{array}{c c} {\mathbf{u}}_0^S(0,z,t) ={\mathbf{u}}_0^S(H,z,t) &{}\quad \text { for } z\in (0,L),\,t\in (0,T),\\ {\mathbf{u}}_1^S(0,z,t) ={\mathbf{u}}_1^S(H,z,t) &{}\quad \text { for } z\in (0,L),\,t\in (0,T),\\ \nabla {\mathbf{u}}_0^S(0,z,t) =\nabla {\mathbf{u}}_0^S(H,z,t) &{}\quad \text { for } z\in (0,L),\,t\in (0,T),\\ \nabla {\mathbf{u}}_1^S(0,z,t) =\nabla {\mathbf{u}}_1^S(H,z,t) &{}\quad \text { for } z\in (0,L),\,t\in (0,T). \end{array} \end{aligned}$$
(48)

Then, we obtain a more general reduced model with the limit part

$$\begin{aligned} \partial _z (U_2^S)' - \partial _z \bigg (K_{av}^{SF}\,\partial _z \Lambda \bigg )&=\Psi |_{x=1} -\Psi |_{x=0} - (u_1^S)'|_{x=1} + (u_1^S)'|_{x=0} ,\nonumber \\ -\partial _z \Lambda + \Bigg ( \int _0^1({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,dx\Bigg ) \partial _{zz}U_2^S&= 0,\nonumber \\ \big ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\big )\partial _{xx} U_1^S|_{z=1}&= - \bigg ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\bigg )\partial _{xz} U_2^S|_{z=1} \end{aligned}$$
(49)

and a corrector part

$$\begin{aligned}&-( u_1^S)' + ({\hat{\lambda }}^S + 2{{\hat{\mu }}}^S)K^{SF} \partial _{xx} u_1^S + ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S) K^{SF}\partial _{xz} u_2^S \nonumber \\&\quad \qquad \qquad = x\,\partial _z (U_2^S)' - \partial _z\left( \partial _x \Lambda \int _0^x K^{SF}\,dx\right) +\Psi |_{x=0} - (u_1^S)'|_{x=0},\nonumber \\&\quad {{\hat{\mu }}}^S \partial _{xx} u_2^S =\partial _z \Lambda - ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S) \partial _{zz}U_2^S- ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{xz} U_1^S|_{z=1}. \end{aligned}$$
(50)

Remark 5

The derivation of the reduced model can be systematically extended to three-dimensional thin domains. For this, similar periodic boundary conditions have to be set for the third thin direction. The resulting limit model will then have an extra one-dimensional equation for the transverse displacement at the load surface, while the corrector model will consist of three three-dimensional equations for the displacement components in the whole domain. However, for one-dimensional domains the reduced model, which simplifies to the limit model, is then the one-dimensional full biphasic model itself.

4 Numerical validation of the reduced model

In this section, we investigate the validity of the reduced model (45), (46) by comparing its numerical solutions with those of the dimensionless full model (15), (16) and (17), which is considered here as a reference of accuracy. For this, we use a finite difference scheme to discretize the spatial derivatives in both models, while the time derivatives are approximated using the Euler method.

4.1 Numerical treatment of the reduced model

We solve the reduced model numerically by discretizing the spatial derivatives using centered finite differences, and the time derivatives of the displacements using the backward Euler method. First, the limit model (45) is solved for the components \(\Lambda \), \(U_2^S\) and \(U_1^S\). Then, the corrector model (46) is solved for the components \(u_2^S\) and \(u_1^S\). The final solution for the reduced model is then given as \(U_2^S+\gamma ^2 u_2^S\) for the vertical displacement, \(U_1^S+\gamma ^2 u_1^S\) for the horizontal displacement and \(\Lambda +\gamma ^2 \lambda \) for the pressure, where \(\lambda =\lambda _2\) satisfies (32).

We define a partition \(0=t^0<t^1<\dots <t^N=T\) of the time interval [0, T] with constant step size \(\Delta t = t^n-t^{n-1}\), for \(n,N\in {\mathbb {N}}\). For the sake of a lighter notation, we keep the continuous operators of the spatial derivatives and integrals. Then, a time-discrete version of the limit model (45) is given as

$$\begin{aligned}&\partial _z \frac{U_2^{S,n+1} - U_2^{S,n}}{\Delta t} - \partial _z \Big (\partial _z \Lambda ^{n+1} K^{SF}_{av}\Big ) = 0, \\&-\partial _z \Lambda ^{n+1} + \left( \int _0^1({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,dx\right) \partial _{zz}U_2^{S,n+1} = 0. \end{aligned}$$

Then, Eq. (45)\(_3\) is solved for \(U_1^{S,n+1}\) using

$$\begin{aligned} \big ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S\big )\partial _{xx} U_1^{S,n+1}|_{z=1} = - \big ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S\big )\partial _{xz} U_2^{S,n+1}|_{z=1}. \end{aligned}$$

Finally, the corrector model (46) for \(u_1^S\) and \(u_2^S\) is solved using

$$\begin{aligned}&- \frac{u_1^{S,n+1}-u_1^{S,n}}{\Delta t}+ ({\hat{\lambda }}^S+2{{\hat{\mu }}}^S)\partial _{xx} u_1^{S,n+1} + ({{\hat{\lambda }}}^S + {{\hat{\mu }}}^S) K^{SF}\partial _{xz} u_2^{S,n+1} \\&\quad \quad \quad \quad = x\,\partial _z \frac{ U_2^{S,n+1}-U_2^{S,n}}{\Delta t} - \partial _z\left( \partial _z \Lambda ^{n+1} \int _0^x K^{FS}\,dx \right) ,\\&\quad {{\hat{\mu }}}^S \partial _{xx} u_2^{S,n+1}=\partial _z \Lambda ^{n+1} - ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{zz}U_2^{S,n+1}- ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\partial _{xz} U_1^{S,n+1}|_{z=1} . \end{aligned}$$

In the following, we present numerical solutions for the reduced model and compare them with those of the dimensionless full model (reference model). Note that the dimensionless domain is the squared unit \((0,1)\times (0,1)\). At the domain’s upper boundary, we choose the load q, see Eq. (11)\(_1\), to be an enforced displacement satisfying

$$\begin{aligned} q = ({{\hat{\lambda }}}^S + 2{{\hat{\mu }}}^S)\,\partial _z u_2=-1~ \text {Pa}, \end{aligned}$$
(51)

where the fluid pressure \(\lambda \) is set to zero, allowing an outflow. The load is also assumed to linearly increase over 10 time increments. For accuracy testing, the solutions of the dimensionless full model are considered in domains with decreasing geometrical parameter \(\gamma \in \{1,1/5,1/10, 1/15\}\).

For the numerical examples, we consider a domain consisting of a porous material with material parameters \(\lambda ^S = 100 \) Pa and \( \mu ^S = 200\) Pa. The material is saturated with a fluid with viscosity \(\mu = 10^{-3}\). The domain has a fixed width \(H=1\) m, and we variate the parameter \(\gamma \) by changing the length \(L=1,\,5,\,10\) or 15 meters. Choosing a material with intrinsic permeability \(K^{SF}=10^{-5} \mathrm{m}^2\) in a domain with length \(L=10\) m leads the factor \(\frac{\psi }{L}=10^{-4}\), see Eq. (14), that maps the real time t to the dimensionless time \({\bar{t}}\). For example, if the end time \(t=2\) seconds with time step size \(\Delta t=0.1\), then end time of the dimensionless full model (and consequently in the reduced model) is \(t = 2\times 10^{-4}\) with time step size \(\Delta t = 10^{-5}\).

Fig. 2
figure 2

Scenario 1: Classical consolidation problem with constant permeability and load applied to whole upper boundary

Full size image

For the comparison of the two models, four scenarios are designed. These have an academic character with different boundary conditions that can occur in applications of bio- and soil mechanics. The first scenario is a classical consolidation problem with constant permeability and load applied to the whole upper boundary. In the second scenario, the load q is applied to the whole upper boundary, but the domain consists of three horizontal layers with different permeabilities, while, in the third scenario, the load is also fully applied to the upper boundary, but the domain consists of three columns with different permeabilities. Finally, the load is applied partially in the fourth scenario, while the permeability is set to be constant. Note that all numerical examples are performed on the same hardware: Intel®Core\(^{\hbox {TM}}\) i7-8565U CPU at 1.80GHz\(\times \)8 with memory 16 GB.

Remark 6

In this section, the accuracy of the reduced model compared to the full solution is verified for small strains in the reference configuration. In fact, further extensions of the reduced model and its numerical scheme to more general cases, such as finite deformations or nonlinear constitutive relation, are part of our future investigations.

Scenario 1: Full load and constant permeability

We consider a classical consolidation problem, such that the porous material has constant permeability \(k^{SF}=10^{-5}\,\text {m}^2\) and load \(q=-1\,\text {Pa}\) is applied downwards to the whole upper boundary of the domain, see Fig. 2. We take this parameter set because it shows a high sensitivity between the solid and fluid load-bearing behavior and thus provides a more challenging task for the reduced model. In this scenario, the reduced model simplifies to the limit model as the right side of the corrector model vanishes, see Remark 3. Note here that the corresponding dimensionless permeability in the dimensionless full model and the reduced model satisfies \(K^{SF} = 1\). For the numerical comparison, we choose a Cartesian grid with \(30\times 30\) squared elements. This choice guarantees a convergence of the solutions for all numerical examples. The time step size is \(\Delta t = 10^{-5}\), and simulation end time is \(2\times 10^{-4}\) seconds.

Fig. 3
figure 3

Comparison of numerical solution for the reduced model (right) and the full model with \(\gamma =1\) (left) and \(\gamma =1/10\) (middle). First row is horizontal displacement, second row is vertical displacement, and the third is pressure in a grid with \(30\times 30\) elements, \(\Delta t = 10^{-5}\) and end time \(2\times 10^{-4}\) seconds

Full size image
Fig. 4
figure 4

Convergence of horizontal displacement (left), vertical displacement (middle) and pressure (right) for the full model at the 1D column \(x=1/6\) in domains with decreasing \(\gamma \in \{1,1/5,1/10,1/15\}\) to the corresponding solution for the reduced model

Full size image

In Fig. 3, we present the numerical solutions for the horizontal displacement \(U_1^S\) (first row), the vertical displacement \(U_2^S\) (second row) and the pressure \(\Lambda \) (third row) for both full and reduced models. To demonstrate the effect of decreasing the geometrical parameter \(\gamma \) on the mechanics in the domain, solutions for the dimensionless full model are presented with \(\gamma =1\) (left column) and \(\gamma = 1/10\) (middle column). Solutions of these two cases are compared with the corresponding solutions for the reduced model (right column). Figure 3 shows that numerical solutions of reduced model are very similar to those of the full model in domains with small geometrical parameter (\(\gamma =1/10\)). For domains with \(\gamma =1\), the vertical displacement and pressure are well approximated by the reduced model, but not the horizontal displacement.

Fig. 5
figure 5

Scenario 2: A sketch of consolidation process in a unit square of three layers with different permeabilities

Full size image

For the sake of preciseness, we compare in Fig. 4 the numerical solutions for the reduced model in the one-dimensional column of the domain \(\{(x,z)|x=1/6\}\) to those for the full model in four domains with fixed width \(H=1\) m, but increasing lengths \(L\in \{1,5,\,10,\,15\}\) meters. Figure 4 shows that the horizontal displacement \(u_1^{S,\gamma }\) converges to that of the reduced model \(U_1^S\) as the domains length L gets larger. However, vertical displacement \(u_2^{S,\gamma }\) and pressure \(\lambda ^{\gamma }\) for the full model are identical to those for the reduced model for all \(\gamma >0\).

Scenario 2: Full load and permeability changes vertically

Here, the accuracy and computational efficiency of the reduced model are investigated in a domain consisting of three horizontal layers with different permeabilities (\(k^{SF} = 10^{-4}\) in the upper and lower layers, while \(k^{SF}= 10^{-6}\) in the middle), see also Fig. 5. Then, the dimensionless permeability is given as

$$\begin{aligned} K^{SF} = \left\{ \begin{array}{cl} 10^{-4}/K^{SF}_m, &{} \quad \text { in } 0<z<\frac{1}{3}\cup \frac{2}{3}<z<1 ,\\ 10^{-6}/K^{SF}_m, &{} \quad \text { in } \frac{1}{3}\le z \le \frac{2}{3} , \end{array}\right. \end{aligned}$$
(52)

where \(K^{SF}_m= \frac{2}{3}\,10^{-4} + \frac{1}{3}10^{-6}\). For this setup, the corrector model has also a zero solution and the reduced model reduces to the limit model (Remark 3).

In Fig. 6, we present the numerical solutions for the horizontal displacement \(U^S_1\) (first row), the vertical displacement \(U^S_2\) (second row) and the pressure \(\Lambda \) (third row) for the reduced model and the full model with different geometrical parameters (\(\gamma =1\) and \(\gamma = 1/10\)). Figure 6 shows that the reduced model is a good approximation of the full model in thin domains (\(\gamma =1/10\)). However, for non-thin domains (\(\gamma =1\)), the reduced model well approximates the vertical displacemet and pressure, but not the horizontal displacement.

Fig. 6
figure 6

Comparison of the reduced model and the full model in layers with different permeabilities. First row is horizontal displacement, second is the vertical displacement, and third is pressure. Here, grid is with \(30\times 30\) elements, \(\Delta t = 10^{-5}\) and end time \(2\times 10^{-4}\) seconds

Full size image
Fig. 7
figure 7

Convergence of horizontal displacement (left) vertical displacement (right) and pressure for the full model in domains with \(H=1\) and increasing length \(L\in \{1,5,10,15\}\) to the those of the reduced model in the 1D column \(x=1/6\) after 20 time loops

Full size image
Fig. 8
figure 8

Pressure distribution at different times levels for the full model (left) and the reduced model (right)

Full size image
Fig. 9
figure 9

CPU time required to solve the full and the reduced model using grid sizes \(\{30\times 30,\,60\times 60,\,120\times 120,\,420\times 420\} \) over 50 time loops

Full size image

For more preciseness, we present in Fig. 7 material displacement and fluid pore pressure in the one-dimensional column of the domain \(\{(x,z)|x=1/6\}\) for both models. For the full model, solutions are computed in domains with fixed width \(H=1\), but increasing lengths \(L\in \{1,5,\,10,\,15\}\). Figure 7 shows, as in the previous scenario, that horizontal displacement \(u_1^{S,\gamma }\) converges to that of the reduced model \(U_1^S\) as the domains’ lengths L get larger, or in other words, the parameter \(\gamma =1/L\) becomes smaller. In addition, the figure shows that the vertical displacements \(u_2^{S,\gamma }\) and pressures \(\lambda ^{\gamma }\) for the full model are almost identical to those for the reduced model, for all \(\gamma >0\).

For more investigation of the accuracy of the reduced model, we study in Fig. 8 the pressure distribution behavior over increasing time levels and compare it with that for the full model. On the left side of Fig. 8 is the pressure distribution for the full model over different times, while on the right is that for the reduced model. It is clear from the figure that the reduced model is able to capture the rate of convergence to steady-state solutions as the full model.

Finally, we investigate in Fig. 9 the computational efficiency of the reduced model. For this, we compare the CPU time required for solving the reduced model with that of the full model in grids with increasing number of elements. The CPU times are summarized in Table 1. The values in the table show a high reduction in computational time by the reduced model. For more clarity, these values are also plotted in Fig. 9, which shows that the CPU time required for solving the full model grows exponentially as the number of grid cells duplicates. On the contrary, the CPU time for the reduced model increases by a factor of 1/3 as the gird size duplicates. This is a consequence of the reduction in the dimension of the linear system resulting from discretizing the reduced model (limit model in this scenario).

Scenario 3: Full load and permeability changes horizontally

In this scenario, we investigate the effect of changing the material’s intrinsic permeability in the horizontal (transverse) direction on the accuracy of the reduced model. Hence, we consider a domain consisting of three columns with different permeabilities (\(k^{SF} = 10^{-4}~\hbox {m}^2\) in the left and right layers, while \(k^{SF}= 10^{-6}~\hbox {{m}}^2\) in the middle), see also Fig. 10. Note that in this setup, the corrector model has no trivial solutions and is part of the reduced model.

In Fig. 11, we present numerical solutions for the horizontal displacement \(U_1^S+\gamma ^2 u_1^S\) (first row), vertical displacement \(U_2^S+\gamma ^2 u_2^S\) (second row) and pressure \(\Lambda + \gamma ^2\lambda \) (third row) for the reduced model (right column) and compare them with the corresponding solutions for the full model with \(\gamma =1\) (left column) and \(\gamma = 1/10\) (middle column). The Cartesian grid has \(30\times 30\) elements, time step size \(\Delta t = 10^{-5}\), and simulation ends after 20 time loops. The figure shows that solutions for the reduced model match very well with those for the full model in domains with small ratio \(\gamma =1/10\), but different from those for the full model with \(\gamma =1\).

For more clarity, we present in Fig. 12 solutions for the full model in the one-dimensional column \(\{(x,z)|x=1/3\}\) of domains with fixed width \(H=1\) but increasing length \(L\in \{ 1,\,5,\,10,\,15 \}\). These solutions are compared with corresponding solution of the reduced model in the same column. The figure shows that solutions for the full model converge to those of the reduced model as the domain’s ratio approaches zero, even when permeability changes in the transverse direction.

Table 1 Comparison of CPU time required for solving the reduced model (limit model) and the dimensionless full model with \(\gamma = 1/10\)
Full size table
Fig. 10
figure 10

Scenario 3: Illustration of consolidation process in a porous domain of three columns with different permeabilities

Full size image
Fig. 11
figure 11

Numerical solutions for the full model (left: \(\gamma =1\) and right: \(\gamma =1/10\)) in a domain of three columns with different permeabilities against the reduced model (right). First row is horizontal displacement, second row is vertical displacement, and third is pressure. Solutions correspond to \(30\times 30\) grid, \(\Delta t = 10^{-5}\) and end time \(2\times 10^{-4}\) seconds

Full size image
Fig. 12
figure 12

Convergence of horizontal displacement (left), vertical displacement (middle) and pressure (right) for the full model to those of the reduced model as domain’s length increases \(L\in \{1,\,5,\,10,\,15\}\). Solutions in the 1D column \(x=1/6\) after 20 time loops are plotted

Full size image
Fig. 13
figure 13

Comparison of vertical displacement for the full model (left) and for the reduced model (right) over various times in the 1D column \(x=1/6\)

Full size image
Fig. 14
figure 14

The volumetric strain for the reduced model (first row) in Scenarios 1, 2 and 3 (left to right) and the full model (second row) with \(\gamma =1/10\). Third row compares both models in the 1D column with \(x=1/6\) after 20 time loops

Full size image

In Fig. 13, we study the convergence rate of the reduced model to the steady state. For this, we present the vertical displacement solutions for the reduced model (right) and the full model (left) overs various times. It is clear from the figure that the reduced model has the same convergence rate to steady-state solutions as the full model.

Fig. 15
figure 15

The von Mises stress for the reduced model (first row) in Scenarios 1, 2 and 3 (left to right) and the full model (second row) with \(\gamma =1/10\). Third row compares both models in the 1D column with \(x=1/6\)

Full size image
Fig. 16
figure 16

CPU time required to solve the full model and the reduced model using grid sizes \(\{30\times 30,\,60\times 60,\,120\times 120,\,420\times 420\} \) over 50 time loops

Full size image

In Fig. 14, we compare the volumetric strain for the reduced model with that for the full model with a small geometrical parameter \(\gamma =1/10\). In the first row (left to right), we present the volumetric strain in the Scenarios 1, 2 and 3, respectively. Similarly, the corresponding volumetric strain for the full model in the three scenarios is presented in the second row. Finally in the third row, the volumetric strain for both models is plotted in the 1D column at \(x=1/6\). The figure shows that the strain is very well approximated by the reduced model in the Scenarios 1 and 2, where no horizontal variations in the permeability occur. However, in Scenario 3, where permeability changes horizontally, the strain of the reduced model is well approximated only for domain with small enough geometrical ratio \(\gamma \).

In Fig. 15, we repeat the numerical comparison from Fig. 14 for the von Mises stress instead of the volumetric strain. It is clear from the figure, that von Mises stress is well approximated by the reduced model in Scenarios 1 and 2. On the contrary to Scenario 3, the stress calculated by the reduced model is not sufficiently close to the that of the full model, even in domains with small parameter \(\gamma \).

Finally, as the reduced model here consists of the limit model and the corrector model, we re-compare in Fig. 16 the computational efficiency of the reduced model with the full model. Here, the CPU time is plotted for both models as the number of spatial grid duplicates. It is noticeable that the CPU time of the reduced model is still significantly less than that of the full model.

Fig. 17
figure 17

Scenario 4: A sketch of a consolidation process with partial load applied to the upper boundary

Full size image

Scenario 4: Partial load and constant permeability We study the effect of applying the load q partially to upper boundary of the domain on the accuracy of the reduced model, see Fig. 17. For this, we define the load function q in Eq. (51) as

$$\begin{aligned} q = \left\{ \begin{array}{rl} -3 ~\mathrm{Pa,}&{} \quad 0< x<\frac{1}{3},\\ 0 ~\mathrm{Pa,}&{} \quad \frac{1}{3}\le x<1. \end{array}\right. \end{aligned}$$
(53)

Then, the mean load value equals \(q_m=1\) Pa. We use constant permeability \(k^{SF}=10^{-5} ~\mathrm{m}^2\) as in Scenario 1. In Fig. 18, we present the numerical solutions for the reduced model (right column), the dimensionless full model with \(\gamma =1\) (left column) and the dimensionless full model with \(\gamma =1/10\) (middle column). In the first row is the material horizontal displacement, in the second is the vertical displacement, and in the third is the fluid pressure. Figure 18 shows a large difference between solutions for the reduced model and the full model with \(\gamma =1\). However, the difference significantly decreases as the geometrical parameter decreases \(\gamma =1/10\).

In Fig. 19, we plot solutions for the reduced model in a one-dimensional column of the domain and compare them to the corresponding solution for the dimensionless full model with decreasing parameters \(\gamma =\{1,\,1/5,\,1/10,\,1/15\}\). In the first row of the figure are the solutions of the models in the column \(x=1/6\), while in the second are the solutions in the middle column \(x=1/2\). It is noticeable that the horizontal displacement solutions for the full model converge to that of the reduced model at both positions. However, for the vertical displacement and pressure, the convergence rate in the middle of the domain (\(x=1/2\)) is much faster than at the sides (\(x=1/6\)).

Fig. 18
figure 18

Numerical solutions for the full model in domains with \(\gamma =1\) (left column) and \(\gamma =1/10\) (middle column) vs. that of the reduced model (right column). Solutions correspond to a \(30\times 30\) grid, after 20 time loops

Full size image
Fig. 19
figure 19

Convergence of numerical solution of the full model with increasing length \(L\in \{1,5,\,10,\,15\}\) to the solution of the reduced model at the 1D columns \(x=1/6\) (first row) and at \(x=1/2\) (second row)

Full size image

5 Conclusion

We have derived a reduced biphasic model for saturated porous materials that undergo poro-elastic deformation in thin domains. The derivation is based on formal asymptotic analysis with respect to the small geometrical parameter \(\gamma >0\) of domain’s width–length ratio. The reduced model is a combination of a limit and a corrector model. The limit model is the asymptotic limit of the full model as the parameter \(\gamma \) tends to zero. It consists of one-dimensional equations that describe the transversely averaged fluid pressure and material displacement in the longitudinal direction. The corrector model is of small order \({\mathcal {O}}(\gamma ^2)\) and describes mechanics caused by variations in the permeability or the applied load in the transverse direction.

We have proposed a numerical scheme, based on the finite difference method, to investigate the validity of the reduced model. For this, we used the dimensionless full model as a reference of accuracy and computational efficiency. For the numerical validation, we have designed several scenarios of a classical consolidation problem. These aim to study the effect of changing permeability or position of the applied load on the accuracy of the model.

The numerical examples demonstrated the high computational efficiency of the reduced model over the full model. In addition to this, the reduced model has reliable solutions that match very well with those for the full model whenever the geometrical parameter is small enough (\(\gamma ={\mathcal {O}}(10^{-1})\)).

Finally, we emphasize that the presented analytical model reduction method is excellently suited to support challenging engineering tasks, for example, in soil mechanics and biomechanics. This includes, e.g., growth and remodeling of living tissue, uncertainty quantification, optimization tasks, as well as fatigue and life cycle investigations.