Model order reduction for deformable porous materials in thin domains via asymptotic analysis

We study fluid-saturated porous materials that undergo poro-elastic deformations in thin domains. The mechanics in such materials are described using a biphasic model based on the theory of porous media (TPM) and consisting of a system of differential equations for material’s displacement and fluid’s pressure. These equations are in general strongly coupled and nonlinear, such that exact solutions are hard to obtain and numerical solutions are computationally expensive. This paper reduces the complexity of the biphasic model in thin domains with a scale separation between domain’s width and length. Based on standard asymptotic analysis, we derive a reduced model that combines two sub-models. Firstly, a limit model consists of averaged equations that describe the fluid pore pressure and displacement in the longitudinal direction of the domain. Secondly, a corrector model re-captures the mechanics in the transverse direction. The validity of the reduced model is finally tested using a set of numerical examples. These demonstrate the computational efficiency of the reduced model, while maintaining reliable solutions in comparison with original biphasic TPM model in thin domain.

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

Full biphasic model
We consider a fluid-saturated porous material ϕ = ϕ S ∪ ϕ F in a thin rectangular domain = (0, H ) × (0, L) that undergoes small elastic deformation. Without loss of generality, and for more illustration, we consider an remain constant, where n α defines the volume fraction of the phase ϕ α with α ∈ {S, F}. The full-saturation condition is expressed with the sum of volume fractions via n S + n F = 1.
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 the balance equation of volume and Darcy's law In Eq. (3), T is the Cauchy stress tensor, T S E is the solid's extra stress and λ 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 where ε S is the linearized strain with and μ S and λ S are the Lamé constants. Consequently, Eq. (3) can be expressed by the solid displacement field u S . In Eqs. (4) and (5), w F is the seepage velocity, k SF the intrinsic (isotropic) permeability and μ F R the fluid's viscosity. The primary variables in the biphasic model are the fluid's pressure λ and the solid displacement u S = (u S 1 , u S 2 ) T , where u S 1 and u S 2 are the displacement components in the horizontal (transverse) and vertical (longitudinal) directions, respectively. We set w F = (w F 1 , w F 2 ), where w F 1 and w F 2 are the seepage velocity components in the horizontal and vertical direction, respectively. Then, the volume balance equation (4) can be written as Note that we scaled second term in the first equation in (8) by γ −2 , where γ = H L is the domain's widthlength 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 γ 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 into the momentum balance equation (3) produces and Throughout the paper, we call Eqs.
where the vector n in (11) 2 is the outer normal to the boundary. We also set periodic boundary conditions on the left and right boundaries The initial conditions for this problem are defined as 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.

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 γ > 0 explicitly. Then, by applying standard asymptotic analysis, based on the small γ , to the dimensionless full model we derive a reduced biphasic model consisting of two main parts: The first is a limit model for γ approaches zeros, while the second is a corrector that re-captures the process in the horizontal direction.

Dimensionless model
We define the geometrical parameter γ = H/L and the dimensionless variables where ψ = . Substituting these variables into the volume balance equation (8), then omitting the par-signs yields Substituting the dimensionless variables (14) into momentum balance equations (9) and (10), then multiplying the resulting equations by H and L, respectively, produces and In these two equations, the par-signs are omitted and the parametersλ S andμ S are the dimensionless Lamé constants defined asλ The Dirichlet and Neumann boundary conditions for the dimensionless model are now summarized as follows where = (0, 1) × (0, 1) is the rescaled domain.

Asymptotic analysis
We consider a solution (u (15)- (17) and define the asymptotic expansions Note that the choice that w F 1 = O(γ ) 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 Equating the terms of the highest power of γ , namely O(γ −2 ), yields Assuming that the intrinsic permeability K SF is strictly positive in the domain, Eq. (21) implies that λ 0 is independent of the horizontal direction Similarly, equating the terms of order O(γ −1 ) implies The terms of order O(γ 0 ) in (20) fulfill the equation Now, we substitute the expansions (19) into the vertical component of Darcy equation (15) 3 and consider the terms of order O(γ 0 ). Then, we have Finally, we substitute (19) into (15) 1 . This leads to Equating the terms of order O(γ −2 ) gives ∂ x u S 1,0 = 0. Then, using the boundary condition (11) 3 , we obtain Similarly, the terms of O(γ −1 ) satisfy However, equating the terms of order O(γ 0 ) gives Setting Eqs. (24) and (25) into (28) produces Integrating this equation over the horizontal direction from 0 to 1 produces 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 Applying the equilibrium result (22) for λ 0 and (35) for u S 2,0 leads to the one-dimensional equation for λ 0 and u S where K SF av = 1 0 K SF dx. Now, we integrate Eq. (29) again over the horizontal direction from 0 to x, which leads to Reordering the terms in the this equation leads to a new formula for the pressure gradient ∂ z λ 2 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 ∂ z λ 2 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 Momentum Balance Equations: We substitute the asymptotic expansions (19) into the momentum balance equation (17). Then, we have Equating the term of order O(γ −2 ) leads to This equation together with the periodic boundary condition (12) implies that Similarly, equating the terms of order O(γ −1 ) leads to Equating the terms of order O(γ 0 ) in (34), then using the result (26) leads to Note that u S 1,0 = 0 inside the dimensionless domain (0, 1)×(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 Then, the periodic boundary condition (12) produces the one-dimensional equation This equation together with (31) represents a coupled system of one-dimensional equations for the primary variables λ 0 and u S 2,0 . Finally, we substitute the asymptotic expansions (19) into the dimensionless momentum balance equation (16). Then, we have Equating the terms of order O(γ 0 ), then using the result (31) produces Note that this equation emphasizes that u S 1,0 = 0 (see Eq. (26)) in the dimensionless domain (0, 1) × (0, 1), as a consequence of (35). However, Eq. (40) includes also the mechanics of u S 1,0 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 S 1,0 at the top boundary of the domain (z = 1) Similarly, the terms of order O(γ ) satisfy This equation, however, implies that u 1,1 = 0 inside the domain (0, 1) × (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 O(γ 2 ) and use (26). Then, we obtain Setting result (32) into this equation produces 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 = λ 0 , U S 2 = u S 2,0 and U S 1 = 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 S 2 = u S 2,2 and u S 1 = u S 1,2 . The reduced model can be written as reduced model = limit model + γ 2 corrector model.
The limit model is a coupled system of the one-dimensional equations while the corrector model is a coupled system of two-dimensional equations Remark 3 It is remarkable that the right side of the corrector model vanishes whenever the permeability K SF and the material parametersλ S andμ 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 where 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 u S 0 = (u S 1,0 , u S 2,0 ) T and second u S 1 = (u S 1,1 , u S 2,1 ) T terms in the asymptotic expansions of the displacement vector u S = u S 0 + γ u S 1 + γ 2 u S 2 are periodic in the horizontal direction, but not necessarily the third term 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 Then, we obtain a more general reduced model with the limit part and a corrector part 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.

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.

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 , U S 2 and U S 1 . Then, the corrector model (46) is solved for the components u S 2 and u S 1 . The final solution for the reduced model is then given as U S 2 + γ 2 u S 2 for the vertical displacement, U S 1 + γ 2 u S 1 for the horizontal displacement and + γ 2 λ for the pressure, where λ = λ 2 satisfies (32). We define a partition 0 = t 0 < t 1 < · · · < t N = T of the time interval [0, T ] with constant step size t = t n − t n−1 , for n, N ∈ 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 Then, Eq. (45) 3 is solved for U S,n+1 Finally, the corrector model (46) for u S 1 and u S 2 is solved using 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 where the fluid pressure λ is set to zero, allowing an outflow. The load is also assumed to linearly increase 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 TM i7-8565U CPU at 1.80GHz×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 m 2 and load q = −1 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 loadbearing 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×30 squared elements. This choice guarantees a convergence of the solutions for all numerical examples. The time step size is t = 10 −5 , and simulation end time is 2 × 10 −4 seconds.
In Fig. 3, 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 (third row) for both full and reduced models. To demonstrate the effect of decreasing the geometrical parameter γ on the mechanics in the domain, solutions for the dimensionless full model are presented with γ = 1 (left column) and γ = 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  For the sake of preciseness, we compare in Fig. 4 the numerical solutions for the reduced model in the onedimensional 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 ∈ {1, 5, 10, 15} meters. Figure 4 shows that the horizontal displacement where K SF m = 2 3 10 −4 + 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 (third row) for the reduced model and the full model with different geometrical parameters (γ = 1 and γ = 1/10). Figure 6 shows that the reduced model is a good approximation of the full model in thin domains (γ = 1/10). However, for non-thin domains (γ = 1), the reduced model well approximates the vertical displacemet and pressure, but not the horizontal displacement.
For more preciseness, we present in Fig. 7 material displacement and fluid pore pressure in the onedimensional 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 ∈ {1, 5, 10, 15}. Figure 7 shows, as in the previous scenario, that horizontal displacement u S,γ 1 converges to that of the reduced model U S 1 as the domains' lengths L get larger, or in other words, the parameter γ = 1/L becomes smaller. In addition, the figure shows that the vertical displacements u S,γ 2 and pressures λ γ for the full model are almost identical to those for the reduced model, for all γ > 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 m 2 in the left and right layers, while k SF = 10 −6 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 S 1 + γ 2 u S 1 (first row), vertical displacement U S 2 + γ 2 u S 2 (second row) and pressure + γ 2 λ (third row) for the reduced model (right column) and compare them with the corresponding solutions for the full model with γ = 1 (left column) and γ = 1/10 (middle column). The Cartesian grid has 30 × 30 elements, time step size 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 γ = 1/10, but different from those for the full model with γ = 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 ∈ {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.
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.
In Fig. 14, we compare the volumetric strain for the reduced model with that for the full model with a small geometrical parameter γ = 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 γ .
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 γ .
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.

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 q = −3 Pa, 0 < x < 1 3 , 0 Pa, 1 3 ≤ x < 1.
(53) Then, the mean load value equals q m = 1 Pa. We use constant permeability k SF = 10 −5 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 γ = 1 (left column) and the dimensionless full model with γ = 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 γ = 1. However, the difference significantly decreases as the geometrical parameter decreases γ = 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 γ = {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).

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 γ > 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 γ tends to zero. It consists of one-dimensional equations that describe the transversely averaged fluid pressure and material displacement 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  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 (γ = 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