Elastoplastic Coupled Model of Saturated Soil Consolidation under Effective Stress

: Soil’s consolidation is a geotechnical problem resulting from a stress-transfer process that initiates when the load is applied to the water contained in the soil, producing a reduction in pore water pressure and rearranging the solid particles, and thus causing a decrease in soil volume. Therefore, consolidation is a coupled ﬂow–mechanical problem. Coupled models have been developed to simulate this phenomenon while considering different theories, providing consistent results. This paper presents an elastoplastic coupled model of consolidation under Terzaghi’s effective stress formulated using the equations of transient ﬂow, balance moment, motion, and the critical state model that considered elastoplastic strains. The coupled model algorithm provided fast and easy results due to its ﬂexibility, as it allowed combinations in loading and boundary conditions. Additionally, it considered the external/internal water ﬂow as an inﬂow or outﬂow, which modiﬁed the pore water pressure and produced changes in the horizontal and vertical displacements. The numerical results obtained showed an appropriate behavior of the consolidation phenomenon, as well as the evolution of the vertical U y and horizontal U x displacements, water pressure p w , volumetric ε v and deviatoric ε q strain, mean σ p and deviatoric σ q stress, volumetric variation ∆ ε v , and elastic/plastic behavior of the ﬁnite elements while considering the yield surface of the critical state.


Introduction
Any infrastructure development such as buildings and roads on saturated compressible cohesive soils (solid particles and water) should be avoided [1] because implications are generated for infrastructure performance and stability [2] due to their low bearing capacity [1], consolidation phenomenon [3], and settlement [4]. Consequently, the study of consolidation is one of the main problems of geotechnical engineering [5,6]. Consolidation is a decrease in soil volume that occurs over a period as a result of the water expulsion from the pores caused by the applied load [1,[6][7][8][9][10][11][12][13] Volume changes result in differential settlement, leading to pavement cracking, pipe failure and foundation damage [1,[14][15][16]. These damages can be more expensive than those associated with natural disasters [17]. The first theory of saturated soil consolidation was proposed by Terzaghi in 1930 [6][7][8][9]12,13,18,19] to predict the change in excess water pore pressure and settlement of the subsoil [10]. Based on this theory, the soil is saturated, Darcy's law is applicable, the relationship between effective stress and strain is linear, and permeability (saturated hydraulic conductivity) and compressibility remain constant [9].
Moreover, most terrestrial materials exhibit behaviors inconsistent with saturated soils [4,10]; therefore, researchers extended the theory of consolidation of saturated soils to unsaturated soils [19]. These soils consist of three phases: solid, water, and air [7,20,21], Water 2022, 14, 2958 3 of 27 fully coupled hydromechanical model under effective stresses that accurately predicted the transition between saturated and unsaturated states as well as between elastic and elastoplastic states, but with limitations in the prediction of volumetric strains. Due to this, there are no analytical solutions for hydromechanical coupling in unsaturated soils; therefore, it is necessary to employ numerical techniques such as the finite element (FE) method [6,7,13,38] in order to simplify the nonlinear behavior of the unsaturated consolidation theory, assuming that water permeability and air transmittance are constant [12]. As a result, several coupled models have been developed to simulate saturated soils consolidation [58][59][60][61]. These models were based on: (a) the soil being saturated; (b) water phase being incompressible; (c) solid phase being incompressible but with a compressible arrangement; and (d) Darcy's law governing the behavior of the flow through the soil. In addition, some authors contributed improvements by incorporating certain theories: the theory of small strains was used by Bentler [58]; the principle of virtual work was used by Manzolillo et al. [61]; and the mass conservation law was considered in the models in [59] and in Manzolillo et al. [61]. The last two employed Terzaghi's effective stress and Galerkin's solution method of weighted residuals, respectively. Additionally, Krishnamoorthy [60] considered the nonlinear behavior of soil and used the hyperbolic relationship of Duncan and Chang [62]. All these models reported consistent results with those obtained in the domain and in the laboratory while showing advantages and disadvantages, suggesting that they can be improved.
The literature shows that consolidation of compressible cohesive soils should be analyzed by means of a fully coupled numerical FE model under the concept of effective stresses able to generate the transition between the saturated and unsaturated states, as well as between the elastic and elastoplastic states, allowing researchers to obtain an adequate consolidation behavior under different boundaries conditions, applied stresses, and water flow into the soil. The above is done to reduce the time and cost of this type of analysis.
The objective of this paper was to formulate and develop the elastoplastic coupled model (flow-mechanical-critical state) for the saturated consolidation of soils under effective stresses while considering the "Moment Balance Laws", "transient flow" (flow) and "Motion" (mechanical) equations in which Terzaghi's "effective stress" was introduced. The inner product between vector functions was used while considering the "Principle of Virtual Work" associated with the "Variational Principle of Minimum Potential Energy" in combination with the finite element method (FEM) and Galerkin's method, which includes a time step that allows the evolution of displacements and water pore pressure as an approximate solution. Additionally, the model was formulated to extend it to unsaturated consolidation, as it permitted the introduction of Bishop's effective stress, which included the Terzaghi effective stress (σ-u a ) and the matrix stress (χψ) representing the product of the chi parameter and the soil suction this last parameter was obtained from the SWRC. The variation of permeability with respect to soil suction (hydraulic conductivity) could also be obtained from the SWRC.

Governing Equations
The applied stress to a saturated soil will cause it to undergo a reduction in soil volume resulting from the expulsion of water from the pores in a time-dependent manner, generating the transient process called consolidation. However, much of the planet is covered by three-phase (unsaturated) soils: solid, liquid, and air. These soils experience variations in their water content produced by climatic conditions that subject them to wetting and drying cycles. Therefore, a numerical procedure is required to obtain the solution to a coupled consolidation problem (saturated-unsaturated) using the finite element method. Thus, the proposed coupled elastoplastic model was formulated and prepared considering the consolidation of saturated soil, but it could be modified for the consolidation of unsaturated soil after coupling the variation in hydraulic conductivity with respect to the degree of saturation obtained from the soil water retention curve (SWRC). Hence, soil mechanical deformation and the flowing water pressure constitute a coupled field problem. Considering the above, the following assumptions were considered: (a) complete saturation and incompressibility of water; (b) incompressible solid phase (soil skeleton); and (c) Darcy's law governed the flow behavior through the soil. Additionally, the theory of small deformations was used. According to continuum mechanics theory, the governing equations are provided below:

Equilibrium Equations
According to the formulation of the elastoplastic model shown in Appendix A of this paper and the balance laws of continuum mechanics and considering the Piola-Kirchhoff stress tensor S in the Lagrangian configuration, this tensor provides the measured force per unit area in the reference configuration [63]: In this study, div is the divergence operator that measures the difference between the outgoing and incoming flux of a vector field on the surface surrounding a control volume. The body forces b 0 = 0 and acceleration field .. χ = 0 were considered because the consolidation process could be regarded as a quasistatic process.

Continuity Equations
In Equation (2), v is the field velocity and . ε is the strain derivative of the field. Thus, Equations (1) and (2) are the basic phenomenological equations of the coupled problem (see Appendix A). This coupling procedure depends on the application of the principle of virtual work, which indicates that the work developed by the internal forces within a system equals the work performed by the external forces acting on the system [64]. This principle is associated with the variational principle of the minimum potential energy [65].

Darcy's Law
where v i is the flow velocity field, k ij is the soil permeability coefficient, γ w is the volumetric mass of water, and h i is the piezometric level, which can be evaluated according to the reference level z; i.e., h i = p w /γ w + z.

Boundary Conditions
The boundary conditions are related to the pore water pressure at all boundaries and to the applied stress load at the surface. Two additional conditions were added: surface water flow produced by meteorological conditions and water flow into the soil due to pipe failure. All boundaries Ω evolve over time and include Dirichlet Ω D and Newman Ω N boundaries. Hence, the boundary conditions can be written in terms of prescribed values: (a) Pore pressure conditions: (b) Applied stress to the surface: (c) Superficial water flow: (d) Water flow into the soil: Water 2022, 14, 2958 5 of 27 (e) Displacements on the soil boundary:

Phenomenological Case of the Solid in the Soil
The inner product was set in Equation (1) while integrating the volume Ω of the body under study. Note that divS ∈ V (vector space V) and such inner product associates a work. For simplicity of writing, φ T is indicated at this point as ϕ, showing that it is a vector function and the product ϕdivS is a scalar function; rewriting in alternate notation gives: Since σ ij,j is equivalent to S, σ ij,j = ∂σ ij /x j is the divergence of the tensor S. Moreover, using the derivative in j of the expression ϕ i σ ij,j according to the derivative product rule, Equation (9) becomes (10): Applying the "Divergence Theorem" [63] for an arbitrary scalar field ϕ with Γ is given as: Substituting (11) into (10) gives: The stress tensor S considers Terzaghi's effective stress equation involving stresses in both the soil skeleton and water contained in soil pores; hence, in index notation it is given as in Equation (13): where S is the soil total stress tensor, S is the soil effective stress, p w is the pore water pressure (hydrostatic stresses), and I is the identity tensor = ∑ i e i ⊗ e i .
Substituting (14) into (12) results in: On the other hand, ϕ i,j associates the derivative of a vector function corresponding in turn to a tensor, whereas any tensor A = E + W can be represented as the sum of a symmetric tensor plus an antisymmetric tensor [63], which in index notation are E = ϕ (i,j) and W = ϕ [i,j] ; therefore, ϕ i,j = ϕ (i,j) + ϕ [i,j] when substituting into the first integral of (15) and recalling that the inner product W·S = 0, since W ∈ Asim.
The term ϕ (i,j) relates the strain tensor in the index notation ε ϕ ij = ϕ (i,j) . The field equations of elastostatics use the constitutive equation S = C[E], which relates the stress tensor to the infinitesimal strain tensor, which in index notation is σ ij = Cε ij Water 2022, 14, 2958 6 of 27 (derived from Hooke's law). Thus, considering the strain tensor and Hooke's law in (16) obtains: By substituting (17) into (15), it is possible to obtain: According to the strain tensor, ε ϕ ij = ϕ (i,j) and it applies for p w δ ij , as δ ij = 0 for i = j. In addition, recalling that σ ij n j = Sn in accordance with the "Cauchy hypothesis" of stress existence [63], Sn = s(n) = t ≡ vector of surface tractions, so Equation (18) takes the form: The finite element approximation (Galerkin) requires: (20) where [N] is the matrix of shape functions and [B] is the matrix of derivatives of shape functions.
Equations (20) and (21) allow their variations to be set in the following form (where δ ij is the "Kronecker Delta"): Applying the "Virtual Work Theorem" [63] and substituting the variations in the fields ϕ and ε in Equation (19) gives: Equation (23) establishes the expression of the virtual work in which a small variation in the function ϕ was applied to generate such work. When substituting the variations according to Equations (22) and the unitary isotropic tensor m [64], which correlates the degrees of freedom of the adjoint functions B ϕ T and N p , it can be written as: because {δφ} T does not involve integration variables. Since it is a common term on both sides of the equality, Equation (24) becomes: The identification of the matrices with the subscripts ϕ, u, and p allowed us to recognize where their approximation originated.

Phenomenological Case of the Flow in the Soil
By integrating the phenomenological Equation (2) in the volume of the body under study and proceeding analogously to the case of soil solids, applying an inner product with an arbitrary function θ to generate a work gives: Denoting for the moment, as was done in (9) ϑ = θ T , and reminding that divv in index notation is divv = v i,i , then (26) transforms into: Considering and applying the derivative product rule in (27) gives: According to the "Divergence Theorem" [63] in the first term of (28), where v i n i = q i , and considering Equation (3), which establishes Darcy's law, results in: Here, ε ii = u i,i ; therefore, u i,i . Furthermore, the piezometric level h = (p w /γ w ) + z can be shown in terms of the water pore pressure p w and the reference head z: The finite element approximation (Galerkin) requires: When deriving Equation (32), the result is: .
indicating that the associated variations in (34) apply the "Principle of Virtual Work" in Equation (30): When substituting (34) into (30), the unit isotropic tensor m gives: In Equation (35), the term V ϑ ,j k ij z ,j dV was not considered because it depends on the reference level at which the piezometric level is evaluated, which for the problem reviewed was not included. On the other hand, substituting the matrix [B ui ] for the product of the unit isotropic tensor in the coupling matrix and extracting the term {δϑ} T from the integrals obtains: (36) Here, [B θ ] has a correspondence analogous to B p . In addition, k ij is associated with the permeability matrix [k].
The final discrete coupled Equations (25) and (36) can thus be written as follows: U} is the nodal displacement velocity vector. The resulting matrices in Equations (37) and (38) suggest the following integral forms (Equation (39)): Matrices [B i ] contain the derivatives of the shape functions; [k] is the permeability matrix, and a vector is needed to couple the values corresponding to the degrees of freedom of the displacement {U} and pore water pressure {p w }.
The coupled model provides a solution for one-dimensional consolidation of saturated soils, which were developed and expressed as Equations (37) and (38), and can be formulated based on the concept of effective stresses. Therefore, coupling the critical state model (Cam-clay) to the above model also yields the mean effective σ p and deviatoric σ q stresses in the elastoplastic behavior of saturated soils [66].
The strain in the critical state model can be evaluated while considering the yield surface f as follows: where M is the slope of the critical state line or stress path in the CD triaxial test method, ϕ is the soil internal friction angle, and σ 0 is the preconsolidation stress, which can be considered a material hardening parameter. In the event of plastic strains, the stress can evolve the hardening law in Equation (42): where ∆ε P v is the variation in the plastic volumetric strain and λ and κ are the slopes of the virgin and unloading sections, respectively, of the compressibility curve of the saturated consolidation test method.
According to Castro Barco [67], the yield surface ( f = 0) represents an ellipse in the σ p -σ q (σ q is the deviatoric stress). The yield surface is the limit of the elastic stress states. These stresses within this limit only give rise to elastic strain increments; therefore, in this case the stiffness matrix is the stiffness matrix of the elastic component [K] = [K e ]. Conversely, while stresses that tend to exceed the limit give rise to increases in plastic strains, in this case, both elastic and plastic strains occur, where the stiffness matrix is a composite of elastic and plastic components [K] = [K e ] + K p . Therefore, the terms of Equations (37) and (38) are maintained and only the stiffness matrix [K] is modified. Then, Equation (37) can be rewritten as Equation (43): The coupled model established in Equations (37) and (38) was formulated in terms of Cartesian stresses σ x , σ y , and σ z and strains ε x , ε y , and ε z , while the critical state provides constitutive equations formulated while considering both volumetric ε v and deviatoric ε q strains. Due to the above, the constitutive Equation (43) must be transformed as a function of both the stresses σ x , σ y , and τ xy and strains ε x , ε y , and γ xy . Hence, since this case involves plane strain, then ε z = 0 even though σ z exists.
where υ and E are Poisson's ratio and the elastic modulus, respectively, of the soil; σ 1 , σ 2 , and σ 3 are the principal stresses applied in the CD triaxial test process corresponding to the stresses σ x , σ y , and σ z , respectively; ε x , ε y , and ε z are the strains along the corresponding directions; and K and G are the volumetric (bulk) and shear moduli, respectively. The last one relates shear stress and angular strain, e 0 and e f are the initial and final void ratios, respectively, and σ apl is the applied surface stress. Similarly, the stiffness matrix of the plastic component originates from the determination of the mean effective σ p and deviatoric σ q stresses as follows: Finally, the constitutive equations of the elastoplastic coupled model of saturated soil consolidation under effective stresses can be obtained as: The resulting matrices in Equations (48) and (49) thus suggest the following integral forms (Equation (50)): To solve Equations (48) and (49), a solution procedure can be applied. The variables include the soil displacement U and the pore water pressure p w at that time. One of the procedures used in this paper entailed the application of the finite difference method in time combined with Galerkin's method to obtain theta values. The approximate values refer to the external forces F acting on the soil surface, external/internal water flow in the soil Q, pore water pressure p w , displacement of the soil mass, and rate of movement in terms of the velocity.
where A is the area of the linear triangular finite element and C V is the soil consolidation coefficient. Finally, the elastoplastic coupled model algorithm was developed in Fortran code to visualize the saturated consolidation phenomenon and its temporal evolution.

Numerical Examples and Results
To validate the elastoplastic coupled model, two numerical examples were examined: (1) a foundation slab covering the entire soil surface, in which vertical flow was allowed only at the bottom boundary, representing the simulation of typical Terzaghi's one-dimensional consolidation; and (2) an isolated footing partially covering the surface, in which vertical flow was allowed at the bottom and top boundaries. In both cases, lateral displacement was restricted.

Case I: Vertical Flow (Bottom Boundary) under Lateral Restrictions
The typical case of one-dimensional consolidation due to a rectangular foundation slab of length L = 5.0 m and infinite width (plane strain state) transferring a uniform load F = −80 kPa (F = −8.0 Ton/m 2 ) onto a clayey soil stratum with thickness H = 3.0 m was analyzed. The U x and U y displacements were restricted both at the lateral boundaries and at the bottom boundary of the soil. In addition, a sand layer occurred below the clay layer, permitting the bottom boundary of the clay layer to drain and indicating a pore pressure of zero.
Similarly, the lateral boundaries did not allow water flow. Therefore, the vertical flow (bottom boundary) was studied. The mesh spacing was set to 0.25 m to use the calibration curves obtained. Figure 1 shows the soil domain under study and Table 1 lists the soil properties needed for the mechanical analysis. In this analysis, the time step was ∆t = 13.93 days.     Figure 2 shows the pore water pressure distribution at time steps of 13.93, 682.57, and 1379.07 days. The pore water pressure distribution gradually decreased with increasing depth because the change in the applied stress distribution also decreased with increasing depth. Thus, at 13.93 days, the pore pressure value reached −79.8 kPa (−7.98 Ton/m 2 ) below the foundation and was very close to the applied stress value of −80 kPa (−8 Ton/m 2 ); however, when approaching the lower boundary, the pore pressure reached zero because flow was permitted. Therefore, when considering only one impermeable boundary, the water contained in the soil domain needed a longer drainage time and covered a greater distance while the filtration rate was very low, resulting in a required analysis time of 1379.07 days.   Figure 2 shows the pore water pressure distribution at time steps of 13.93, 682.57 1379.07 days. The pore water pressure distribution gradually decreased with increa depth because the change in the applied stress distribution also decreased with increa depth. Thus, at 13.93 days, the pore pressure value reached −79.8 kPa (−7.98 Ton/m 2 low the foundation and was very close to the applied stress value of −80 kPa (−8 Ton however, when approaching the lower boundary, the pore pressure reached zero bec flow was permitted. Therefore, when considering only one impermeable boundary water contained in the soil domain needed a longer drainage time and covered a gr distance while the filtration rate was very low, resulting in a required analysis tim 1379.07 days. Additionally, when comparing Figure 2a-c, we observed that the pore pressure decreased over time until a value close to zero was reached. This behavior was appropriate and correct when analyzing the soil consolidation problems, which evolved over time. In addition, the Mandel-Cryer phenomenon as described by Conte [7] was observed due to the sensitivity of the model, which caused the excess pressure to increase to infinity. For this reason, it was required to calibrate the model and obtain calibration curves to reduce this effect by determining the imminent point at which the pore water pressure began to dissipate. Figure 3 shows the symmetrical behavior of horizontal displacement U x . The foundation produced displacement to the right (positive) and left (negative) from the centerline of the domain. The largest displacements occurred in the top third (near the top boundary) while displacement was restricted at the lateral and bottom boundaries; hence, the displacements were very small or negligible. As shown in Figure 3a-c, the maximum value of −0.0007 m (−0.700 mm) was observed at 13.93 days, but over time, this value increased to −0.000713 m (−0.713 mm). This indicated that the largest displacements occurred at the beginning of the consolidation process; i.e., at the time of foundation emplacement. Subsequently, the displacements decreased due to soil rearrangement, which was attributed to pore pressure dissipation throughout the domain. boundary) while displacement was restricted at the lateral and bottom boundaries; hence, the displacements were very small or negligible. As shown in Figure 3a-c, the maximum value of −0.0007 m (−0.700 mm) was observed at 13.93 days, but over time, this value increased to −0.000713 m (−0.713 mm). This indicated that the largest displacements occurred at the beginning of the consolidation process; i.e., at the time of foundation emplacement. Subsequently, the displacements decreased due to soil rearrangement, which was attributed to pore pressure dissipation throughout the domain. The behavior of vertical displacement is shown in Figure 4. It was evident that the vertical displacements at the top boundary were greater and decreased with increasing depth due to the distribution of the applied stress. As shown in Figure 3a, during the first 13.93 days after the beginning of the consolidation process, in the superficial corners, the magnitude of the displacements at −0.01579 m (−1.579 cm) was larger than that of the displacements in the central zones at −0.008755 m (−0.8755 cm), producing a concave profile of superficial displacements due to the large FE size; thus, with finer meshing, this effect could be reduced. This effect could only be observed through the coupled model and could not be observed via classical mathematical functions. In contrast, the displacements near the bottom boundary were already very small or zero. In addition, at 682.57 days (equivalent to 50% of the analysis time), the maximum displacements changed to 3.368 cm, and at the end of the analysis time, the displacements reached 3.373 cm (as shown in Figure 4b). This indicated that significant displacements could occur after 50% of the consolidation time, which was consistent with Terzaghi's consolidation theory. The behavior of vertical displacement U y is shown in Figure 4. It was evident that the vertical displacements at the top boundary were greater and decreased with increasing depth due to the distribution of the applied stress. As shown in Figure 3a, during the first 13.93 days after the beginning of the consolidation process, in the superficial corners, the magnitude of the displacements at −0.01579 m (−1.579 cm) was larger than that of the displacements in the central zones at −0.008755 m (−0.8755 cm), producing a concave profile of superficial displacements due to the large FE size; thus, with finer meshing, this effect could be reduced. This effect could only be observed through the coupled model and could not be observed via classical mathematical functions. In contrast, the U y displacements near the bottom boundary were already very small or zero. In addition, at 682.57 days (equivalent to 50% of the analysis time), the maximum U y displacements changed to 3.368 cm, and at the end of the analysis time, the displacements reached 3.373 cm (as shown in Figure 4b). This indicated that significant displacements could occur after 50% of the consolidation time, which was consistent with Terzaghi's consolidation theory.   Figure 5 shows the evolution of the profile of the vertical displacements of the superficial nodes, revealing that the vertical displacements of the superficial boundary produced by the foundation were negative; i.e., settlement occurred. Figure 5 shows different soil superficial profiles for each time increment (a color line was obtained every 13.93 days); i.e., after 13.93 days, the foundation produced a settlement of −0.80 cm in the center, while in the corner it was −1.50 cm (superior magenta line). At the end of the time analysis, the settlements reached a magnitude of −2.65 cm in the center and −3.30 cm in the corner (bottom magenta line). Based on this finding, we deduced that the soil stratum volume decreased. This behavior was consistent with the type of analysis performed.  Figure 5 shows the evolution of the profile of the vertical displacements of the superficial nodes, revealing that the vertical displacements of the superficial boundary produced by the foundation were negative; i.e., settlement occurred. Figure 5 shows different soil superficial profiles for each time increment (a color line was obtained every 13.93 days); i.e., after 13.93 days, the foundation produced a settlement of −0.80 cm in the center, while in the corner it was −1.50 cm (superior magenta line). At the end of the time analysis, the settlements reached a magnitude of −2.65 cm in the center and −3.30 cm in the corner (bottom magenta line). Based on this finding, we deduced that the soil stratum volume decreased. This behavior was consistent with the type of analysis performed. duced by the foundation were negative; i.e., settlement occurred. Figure 5 shows different soil superficial profiles for each time increment (a color line was obtained every 13.93 days); i.e., after 13.93 days, the foundation produced a settlement of −0.80 cm in the center, while in the corner it was −1.50 cm (superior magenta line). At the end of the time analysis, the settlements reached a magnitude of −2.65 cm in the center and −3.30 cm in the corner (bottom magenta line). Based on this finding, we deduced that the soil stratum volume decreased. This behavior was consistent with the type of analysis performed.  Figure 6 shows the finite element mesh of the soil stratum and thus the behavior (elastic and/or plastic) they were under according to the value of the parameter . When considering this meshing, the behavior of the FEs within the soil stratum could be identified at any time step during the application of the foundation-applied stress. As shown in Figure 6a, for Δt = 13.93 days, 90% of the finite elements within the domain exhibited elastic behavior (stress and strain) or blue zones. Over time (as shown in Figure 6b), 75% of the FEs exhibited plastic behavior (red zones), and at the end of the analysis time; i.e., 1379.07 days (as shown in Figure 6c), all FEs exhibited plastic behavior. Because the low permeability impeded the pore pressure dissipation process and due to the lateral restrictions; i.e., at the bottom and top, the soil was subjected to plastic strains (irreversible) at the end of the consolidation process.  Figure 6 shows the finite element mesh of the soil stratum and thus the behavior (elastic and/or plastic) they were under according to the value of the parameter f . When considering this meshing, the behavior of the FEs within the soil stratum could be identified at any time step during the application of the foundation-applied stress. As shown in Figure 6a, for ∆t = 13.93 days, 90% of the finite elements within the domain exhibited elastic behavior (stress and strain) or blue zones. Over time (as shown in Figure 6b), 75% of the FEs exhibited plastic behavior (red zones), and at the end of the analysis time; i.e., 1379.07 days (as shown in Figure 6c), all FEs exhibited plastic behavior. Because the low permeability impeded the pore pressure dissipation process and due to the lateral restrictions; i.e., at the bottom and top, the soil was subjected to plastic strains (irreversible) at the end of the consolidation process.

Case II: Vertical Flow (Top and bottom Boundary) under Lateral Restrictions
Additionally, a rectangular isolated footing of length = 2.5 and infinite width (plane strain state) was analyzed. The footing was placed at the center below the top boundary of the domain and applied a uniform stress = −80 ( = −8.0 / ) on the clayey soil stratum with thickness = 3.0 .
The top boundary permitted water flow in the unloaded zones (to the left and right of the footing). Similarly, the bottom boundary also allowed water flow, indicating that the pore pressure at these boundaries was zero.
displacements were restricted at the lateral boundaries, while at the bottom soil boundary, and were also restricted. The top boundary permitted water flow in the unloaded zones (to the left and right of the footing). Similarly, the bottom boundary also allowed water flow, indicating that the pore pressure at these boundaries was zero. U x displacements were restricted at the lateral boundaries, while at the bottom soil boundary, U x and U y were also restricted. The meshing was set to 0.25 m in order to use the calibration curves obtained. Figure 7 shows the study domain loaded by the footing, while the soil properties needed for mechanical analysis are listed in Table 2. According to the calibration curve (see Figure A3 in Appendix B), to perform the analysis, the time step was set to ∆t = 4.385 days, while the analysis time was set to t = 877 days. Through numerical analysis, Figure 8 reveals that the pore water pressure of the soil subjected to the above conditions exhibited a very consistent behavior. After the analysis time and time step of the calibration curve were obtained (see Figures A3 and A4 in Appendix B), it was observed that the pore water pressure did not approach infinity. Conversely, it was evident that the pressure began to decrease over time. Thus, Figure 8a shows that after 4.385 days of consolidation, the pore water pressure was −57.89 kPa (−5.789 Ton/m 2 ), which represented a reduction of 27.64% with respect to the stress applied by the footing (−80 kPa). Similarly, after 434.115 days (as shown in Figure 8b), the pore water pressure decreased 99.65% (−0.276 kPa) compared to the applied stress. At the end of the analysis time; i.e., 872.615 days (as shown in Figure 8c), it was evident that the pore water pressure approached zero, representing a typical consolidation phenomenon. In addition, Figure 8 shows that the pore water pressure at the top boundary, which was not loaded (to the left and right of the footing), reached zero. Similarly, the bottom boundary exhibited the same behavior because these boundaries allowed water to drain and the pore pressure within the domain to dissipate. Figure 9 shows the evolution of the horizontal displacements with the consolidation time. Figure 9 shows symmetric behavior of the horizontal displacements . On the basis of performed numerical analyses, we obtained, at the center of the top boundary,  Through numerical analysis, Figure 8 reveals that the pore water pressure of the soil subjected to the above conditions exhibited a very consistent behavior. After the analysis time and time step of the calibration curve were obtained (see Figures A3 and A4 in Appendix B), it was observed that the pore water pressure did not approach infinity. Conversely, it was evident that the pressure began to decrease over time. Thus, Figure 8a shows that after 4.385 days of consolidation, the pore water pressure was −57.89 kPa (−5.789 Ton/m 2 ), which represented a reduction of 27.64% with respect to the stress applied by the footing (−80 kPa). Similarly, after 434.115 days (as shown in Figure 8b), the pore water pressure decreased 99.65% (−0.276 kPa) compared to the applied stress. At the end of the analysis time; i.e., 872.615 days (as shown in Figure 8c), it was evident that the pore water pressure approached zero, representing a typical consolidation phenomenon.
versely, it was evident that the pressure began to decrease over time. Thus, Figure 8a shows that after 4.385 days of consolidation, the pore water pressure was −57.89 kPa (−5.789 Ton/m 2 ), which represented a reduction of 27.64% with respect to the stress applied by the footing (−80 kPa). Similarly, after 434.115 days (as shown in Figure 8b), the pore water pressure decreased 99.65% (−0.276 kPa) compared to the applied stress. At the end of the analysis time; i.e., 872.615 days (as shown in Figure 8c), it was evident that the pore water pressure approached zero, representing a typical consolidation phenomenon. In addition, Figure 8 shows that the pore water pressure at the top boundary, which was not loaded (to the left and right of the footing), reached zero. Similarly, the bottom boundary exhibited the same behavior because these boundaries allowed water to drain and the pore pressure within the domain to dissipate. Figure 9 shows the evolution of the horizontal displacements with the consolidation time. Figure 9 shows symmetric behavior of the horizontal displacements . On the basis of performed numerical analyses, we obtained, at the center of the top boundary, produced displacements to the left (negative) and right (positive) from the centerline of the stratum. After 4.385 days (as shown in Figure 9a), the maximum value ranged from −0.004102 m (−0.4102 cm) to 0.004205 m (0.4205 cm) just at the strip below the footing. However, with increasing depth, the displacements were very small and almost negligible. Additionally, over time (as shown in Figure 9b,c), this value decreased to 0.00314 m (0.314 cm). This indicated that the largest displacements occurred at the beginning of In addition, Figure 8 shows that the pore water pressure at the top boundary, which was not loaded (to the left and right of the footing), reached zero. Similarly, the bottom boundary exhibited the same behavior because these boundaries allowed water to drain and the pore pressure within the domain to dissipate. Figure 9 shows the evolution of the horizontal displacements U x with the consolidation time. Figure 9 shows symmetric behavior of the horizontal displacements U x . On the basis of performed numerical analyses, we obtained, at the center of the top boundary, produced displacements to the left (negative) and right (positive) from the centerline of the stratum. After 4.385 days (as shown in Figure 9a), the maximum value ranged from −0.004102 m (−0.4102 cm) to 0.004205 m (0.4205 cm) just at the strip below the footing. However, with increasing depth, the displacements were very small and almost negligible. Additionally, over time (as shown in Figure 9b,c), this value decreased to 0.00314 m (0.314 cm). This indicated that the largest U x displacements occurred at the beginning of the consolidation process. Thus, when the consolidation time exceeded 50% of the analysis time, the displacements continued to increase but at a smaller magnitude due to rearrangement of the solid soil particles and pore water pressure dissipation throughout the domain. the consolidation process. Thus, when the consolidation time exceeded 50% of the analysis time, the displacements continued to increase but at a smaller magnitude due to rearrangement of the solid soil particles and pore water pressure dissipation throughout the domain.   202 cm). However, in the unloaded zones, the vertical displacements were positive, with a value of +0.006679 m (+0.6679 cm). This indicated that the isolated footing settled across the entire contact surface with the soil, which produced a combination of vertical and horizontal displacements within the domain that generated, in turn, upward (positive) vertical displacements in the unloaded zones at the top boundary. Additionally, after 434.115 (Figure 10b) and 872.615 days (Figure 10c), the footing generated settlement within the central strip of the domain; while at the unloaded surface, lateral boundaries, and bottom boundary, the displacements were zero. Thus, by analyzing the displacements shown in Figure 10, it could be inferred that the domain expe-   (Figure 10c), the footing generated settlement within the central strip of the domain; while at the unloaded surface, lateral boundaries, and bottom boundary, the displacements were zero. Thus, by analyzing the U y displacements shown in Figure 10, it could be inferred that the domain experienced settlement in the loaded zone and uplift in the unloaded zones, indicating permanent strains.
the top boundary of the domain; i.e., the loaded zone (central strip and interior of the domain), the displacements were negative (typical of loaded zones) and reached a maximum value of −0.01202 m (−1.202 cm). However, in the unloaded zones, the vertical displacements were positive, with a value of +0.006679 m (+0.6679 cm). This indicated that the isolated footing settled across the entire contact surface with the soil, which produced a combination of vertical and horizontal displacements within the domain that generated, in turn, upward (positive) vertical displacements in the unloaded zones at the top boundary. Additionally, after 434.115 (Figure 10b) and 872.615 days (Figure 10c), the footing generated settlement within the central strip of the domain; while at the unloaded surface, lateral boundaries, and bottom boundary, the displacements were zero. Thus, by analyzing the displacements shown in Figure 10, it could be inferred that the domain experienced settlement in the loaded zone and uplift in the unloaded zones, indicating permanent strains. The above is clearly shown in Figure 11. This figure shows the variation in the soil superficial profile of the vertical displacements of the superficial nodes, providing a clear view of the evolution of the vertical displacements at the top boundary produced by load application to this boundary and pore water pressure dissipation throughout the interior of the domain (soil). Similarly, we observed how the superficial profile was strained, so we could deduce that the loaded area experienced settlement (displacements below the The above is clearly shown in Figure 11. This figure shows the variation in the soil superficial profile of the vertical displacements of the superficial nodes, providing a clear view of the evolution of the vertical displacements at the top boundary produced by load application to this boundary and pore water pressure dissipation throughout the interior of the domain (soil). Similarly, we observed how the superficial profile was strained, so we could deduce that the loaded area experienced settlement (displacements below the red dotted line), revealing a decrease in volume; while the unloaded area experienced uplifts (displacements above the red line), revealing an increase in volume.   Figure 12 shows the meshing domain of the FEs and indicates whether the behavior of the FEs is elastic or plastic. This determination was based on the yield surface ; for ≤ 0, the FEs exhibited elastic behavior, but for > 0, the FEs exhibited plastic behavior. As a result, at a time step of 4.385 days, a significant portion of the finite elements within the domain experienced plastic stress and strain behavior (red zones in Figure 12a). Over time until the analysis time was reached, the finite elements that exhibited plastic behavior (below the unloaded zones) changed to exhibit elastic behavior (change from red to blue in Figure 12b,c).  Figure 12 shows the meshing domain of the FEs and indicates whether the behavior of the FEs is elastic or plastic. This determination was based on the yield surface f ; for f ≤ 0, the FEs exhibited elastic behavior, but for f > 0, the FEs exhibited plastic behavior. As a result, at a time step of 4.385 days, a significant portion of the finite elements within the domain experienced plastic stress and strain behavior (red zones in Figure 12a). Over time until the analysis time was reached, the finite elements that exhibited plastic behavior (below the unloaded zones) changed to exhibit elastic behavior (change from red to blue in Figure 12b,c). of the FEs is elastic or plastic. This determination was based on the yield surface ; for ≤ 0, the FEs exhibited elastic behavior, but for > 0, the FEs exhibited plastic behavior. As a result, at a time step of 4.385 days, a significant portion of the finite elements within the domain experienced plastic stress and strain behavior (red zones in Figure 12a). Over time until the analysis time was reached, the finite elements that exhibited plastic behavior (below the unloaded zones) changed to exhibit elastic behavior (change from red to blue in Figure 12b Based on the above, a direct relationship could be observed between the vertical and horizontal displacements because the area below the foundation (loaded) exhibited plastic behavior. Furthermore, the secondary effects produced by loading also subjected the unloaded areas to plastic behavior. The reason for this phenomenon was that Based on the above, a direct relationship could be observed between the vertical U y and horizontal U x displacements because the area below the foundation (loaded) exhibited plastic behavior. Furthermore, the secondary effects produced by loading also subjected the unloaded areas to plastic behavior. The reason for this phenomenon was that the soil flowed downward in the central loaded area, displacing the soil toward the sides and finally lifting the sides and deforming the surface soil outward in the unloaded areas. When the pore water pressure was very low-practically zero (at 434.115 days)-the finite elements no longer exhibited plastic behavior but changed to exhibit elastic behavior and remained so until the end of the time analysis process. The reason for this phenomenon was that there were no longer any overpressure-producing permanent strains.
Notably, the proposed elastoplastic coupled model also captured the behavior of (a) the volumetric strain ε v , (b) the deviatoric strain ε q , (c) the mean stress σ p , (d) the deviatoric stress σ q , and (e) the volumetric variation ∆ε v ; the latter reflected the final value compared to the total volume change in the domain as either positive (an increase in volume) or negative (a decrease in volume).
Finally, a comparative analysis between the case studies in the present work revealed that under the same soil properties, magnitude of the applied stress, and restrictions at the lateral and bottom boundaries in both cases, but with different permeable boundaries, the analysis time increased from Cases I to II (872.615 and 1393 days, respectively) and the time step also increased from 4.385 to 13.93 days.
Based on Table 3, we deduced that when flow was permitted at the top boundary, the time needed for the pore water pressure to decrease or dissipate was less than that in the analysis case that did not allow surficial flow. This result was congruent with the behavior of the consolidation phenomenon under both conditions; i.e., the time needed for the pore water pressure to dissipate was shorter because water flowed toward the two boundaries, which suggested that the travel distance of water through the soil was also shorter; hence, the pore water pressures were lower. In the opposite case, with the top boundary sealed, water traveled a greater distance to be discharged, and as a result, the pore water pressures were higher.

Discussion
When analyzing the pore water pressure (p w ) behavior, it was possible to observe the Mandel-Cryer phenomenon described by Conte [7], which explains that an excess pressure exists at the beginning of such a process due to the applied stress; this excess pressure increases and then begins to decrease until it is completely dissipated. The proposed model showed this excess pressure but it increased to infinity; for this reason, it was necessary to calibrate the model and obtain the calibration curves shown in Appendix B. Thus, the behavior of the pore water pressure was consistent with those reported in the literature reviewed. However, a direct comparison was not possible because the reviewed models showed the analysis of a specific point of the stratum while the proposed model analyzed the entire soil stratum by means of color maps.
Regarding horizontal displacements (U x ), the classical theory of consolidation describes that these displacements are zero; however, the proposed model demonstrated small displacements, but these did exist. Likewise, after analyzing the horizontal displacements, it was observed that in the two cases of study, the settlements were symmetrical when tracing a line center. In addition, it was remarkable to identify higher displacements occurring in the first quarter of the analysis time that were becoming lower and lower due to this. This behavior was correct and consistent with the consolidation phenomenon.
The vertical displacements (U y ) in the two case studies showed higher-magnitude displacements occurring in the superficial zone (below the foundation) that decreased as the depth increased. This was consistent with the distribution of the applied stress because it decreased as the depth increased, having less of an effect each time. Likewise, the displacements evolved over time, presenting those of higher magnitude at the beginning of the process, but as time progressed, they decreased considerably. Additionally, due to the Mandel-Cryer phenomenon and the combination of the vertical and horizontal displacements, it was possible to observe the superficial profile of the soil at each time step. The profile revealed that when the soil was subjected to an applied stress, it experienced settlements (negative displacements), but in nonloaded zones (around the foundation), expansions or swellings (positive displacements) were manifested. Such behavior cannot be determined using classical theories, nor did the literature models reviewed report it.
On the other hand, the elastoplastic model evaluated the value of the yield surface f to identify at each time step the behavior (elastic and/or plastic) of each finite element in the mesh. As a result of this analysis, at the beginning of the process in the two case studies, the FEs under the foundation were in elastic behavior (blue zones); then the FEs evolved to plastic behavior (red zones). Therefore, at the end of the analysis process, the loaded zones experienced permanent strains. Likewise, in the periphery of the foundation, the FEs were under plastic behavior, which evidenced that the swellings on the surface were permanent strains. Finally, the calibration curves revealed that when there were two permeable vertical boundaries instead of one, the analysis time and time step were shorter. This indicated that the consolidation process would take less time because the water could flow through two boundaries and permit the water pore pressure to dissipate faster. This behavior is typical in consolidation processes and was in accordance with the classical theory and literature models reviewed.

Conclusions
In this paper, an elastoplastic coupled model (flow-mechanical-critical state) was formulated and developed to study the consolidation of saturated soils under effective stresses and their transition from elastic to elastoplastic behavior and/or vice versa. Based on this investigation, the following conclusions could be derived: The variation in consolidation in saturated soils was determined by the rate at which the fluid pressure in the soil could be reduced. Thus, the process of strain (mechanical behavior) and the distribution of the water pressure (hydraulic behavior) in soils is a coupled flow-mechanical process. Thus, the use of the coupled model represented an integral solution to the coupling problem, yielding approximate solutions that reproduced the hydromechanical behavior of the soil in a complete way.
The proposed elastoplastic coupled model indicated that the resulting coupling equations confirmed that the approach and formulation were adequate and exhibited correct coupling and congruence in the matrix product process.
Simulations using the proposed elastoplastic coupled model algorithm provided fast and easy results due to its flexibility, since it permitted an infinite combination of conditions: (1) boundaries, (2) loads, and (3)  The numerical results revealed an appropriate behavior of the consolidation phenomenon under the conceptual framework of the effective stress of saturated soil, supporting its feasibility and reliability. The water pore pressure behavior presented an overpressure at the beginning (Mandel-Cryer effect), but as time passed it tended to zero. The vertical and horizontal displacements were greater below the foundation and decreased with depth and decreased with time. The soil superficial profile evolved with time and showed the permanent strains at the surface (settlements and/or swellings). The model showed that the FE in the zones below the foundation were initially in elastic behavior, but at the end of the process, they evolved to the plastic regime and showed permanent strains. Therefore, the model was consistent with Terzaghi's consolidation theory, and if two boundaries were permeable instead of one, the analysis time was shorter and would permit water to flow more quickly.

Future Work
The analysis process does not end here: the different combinations of boundary conditions, loading levels, flow patterns, and soil properties considered in the many simulations we performed to improve the model must be examined. However, due to space limitations, only the most recurrent and significant consolidation problems were included.
The proposed elastoplastic coupled model was formulated and prepared while considering saturated soil consolidation but could be modified for unsaturated soil consolidation after coupling the variation in the hydraulic conductivity with respect to the saturation degree obtained from the soil water retention curve (SWRC). This process should be incorporated into the permeability matrix.  During a movement, mechanical interactions between parts of a body or between a body and its environment are described by five types of forces: (1) contact forces between separate parts of a body; (2) contact forces exerted at the boundary of a body by its environment; (3) body forces exerted at points inside the body by the environment; (4) contact forces exerted at points inside the body by the environment; and (5) contact forces exerted at points inside the body by the environment [63].
One of the most important and far-reaching axioms in the continuum mechanics is Cauchy's Hypothesis about the mode of contact forces. Thus, Cauchy assumed the existence of a density of surface forces s(n, x, t) defined by each vector n and all (x, t) on the trajectory T of the motion ( Figure A1a). The field has the following properties: let L be a surface oriented in B t with a positive unit normal vector n at the point x. Then, s(n, x, t) is the force per unit area exerted on the surface L of the material on the negative side of L by the material on the positive side. In fact, if C is an oriented surface tangent to L at x and has the same positive unit normal vector, then the force per unit area at x is the same at both C and L ( Figure A1b). by the material on the positive side. In fact, if is an oriented surface tangent to at and has the same positive unit normal vector, then the force per unit area at is the same at both and ( Figure A1b). Figure A1. Existence of the surface force density ( , , ) between the surfaces and [63].
To determine the contact forces between two separate parts and D at time (Figure A2a), it is sufficient to integrate at time over the contact surface. Thus, Equation (53) gives the force exerted on at D at time . Here, is the external unit normal vector for at .  To determine the contact forces between two separate parts P and D at time t (Figure A2a), it is sufficient to integrate at time t over the contact surface. Thus, Equation (A1) gives the force exerted on P at D at time t. Here, n is the external unit normal vector for ∂P t at x. To determine the contact forces between two separate parts and D at time (Figure A2a), it is sufficient to integrate at time over the contact surface. Thus, Equation (53) gives the force exerted on at D at time . Here, is the external unit normal vector for at .
For points on the boundary of , ( , , )-with the external unit normal vector for at -gives the surface force per unit area applied to the body by the environment. This force is referred to as the surface tractions. In either case, given a part , it represents the total contact forces exerted on at time ( Figure A2b).
( ) The environment can exert forces at points inside the body; such is the case of gravity. These forces are determined by the vector field on the trajectory T; ( , ) gives the force per unit volume exerted by the surroundings on . Thus, (55) gives that part of the force caused by the environment on .
Following the above discussion, we can denote as the surface forces, as the body forces, and the force ( , ) on a part at time by: Thus, the basic axioms connecting motion and force are the momentum balance laws. These laws state that for every part at time : Figure A2. (a) Contact forces; (b) forces at points in the interior caused by the surrounding environment [63].
For points on the boundary of B t , s(n, x, t)-with n the external unit normal vector for ∂B t at x-gives the surface force per unit area applied to the body by the environment. This force is referred to as the surface tractions. In either case, given a part P, it represents the total contact forces exerted on P at time t ( Figure A2b).
The environment can exert forces at points inside the body; such is the case of gravity. These forces are determined by the vector field b on the trajectory T; b(x, t) gives the force per unit volume exerted by the surroundings on x. Thus, (A3) gives that part of the force caused by the environment on P.
Following the above discussion, we can denote s as the surface forces, b as the body forces, and the force f(P, t) on a part P at time t by: Thus, the basic axioms connecting motion and force are the momentum balance laws. These laws state that for every part P at time t: An obvious consequence of (A5), and considering that the linear momentum of a body B is the same as that of a particle of mass m(B) associated to the center of mass of B, is: ..
Thus, by introducing (A6) into (A4), the momentum balance law can be written as: Cauchy's Theorem: This theorem is one of the central results of continuum mechanics. Its main statement is that s(n) is linear in n. It specifies the existence of stresses. It says: let (s, b) be a system of forces for body B during a motion. Then, a necessary and sufficient condition to satisfy the momentum balance laws is the existence of a tensor field T (named Cauchy stress) such that [63]: (a) For each unit vector n s(n) = T(n) (A9) (b) The tensor T is symmetric; (c) The tensor T satisfies the motion equation: Piola-Kirchhoff Stress Tensor: However, the Cauchy's stress tensor T measures the contact forces per unit area in the deformed configuration. In many problems of interest, it is not convenient to work with the Cauchy's stress tensor T because the deformed configuration is not known in advance. For this reason, the stress tensor that gives the measured force per unit area in the reference configuration was introduced [63].
Let (X, T) be a dynamic process. Then, given the part P, writing it in terms of the total surface force on P at time t as an integral over ∂P gives: This stress tensor S is denoted as the Piola-Kirchhoff stress tensor in (A12) and Sn in (39) is the surface force per unit area measured in the reference configuration [63].
Similarly, if b is the body forces corresponding to (X, T), then: On the other hand, let F be the gradient tensor of the deformation function f , which indicates: The determinant of F establishes the ratio between the volume of the deformed field with respect to the original volume [63]. This gives: This expression is general and involves the body with its internal phases. Likewise, when evaluating V f y V i for a unit volume, we obtain: where ε is the volumetric strain of the body.
Deriving expression (A19) in time, we obtain: However, according to the theorem of the "Transport of Volume" [63]: where v is the velocity of the particle; therefore: Equations (A16) and (A23) were the basic phenomenological equations for this study. The procedure to determine the coupling equations consisted of the application of the "Principle of Virtual Work", which expresses "The work developed by the internal force in a system is equal to the work developed by the external forces acting on it". This principle associates the "Variational Principle of Minimum Potential Energy". We proceeded in this study in an analogous way and established an inner product through vector functions.

Calibration Curves of the Elastoplastic Coupled Model (Flow-Mechanical-Critical-State)
A calibration was conducted to identify the variables that significantly influenced the results. To accomplish this task, the following variables were considered: (a) the saturated permeability coefficient k y , (b) the type of lateral restriction, and (c) the time step ∆t. In contrast, θ = 0.9 was defined as an adequate value to prevent numerical oscillations and the elastic modulus E was defined by Equation (46). Finally, Figures A3 and A4 show the calibration curves for the different combinations of flow conditions with and without lateral restrictions, respectively. When analyzing the results obtained from the simulations and as shown in Figure A3, the dependence of the time step ∆t and analysis time t on the permeability coefficient k was evident. For example, for k = 8.64 × 10 −4 m/day (k = 1 × 10 −8 m/s), the time step should be ∆t = 13.93 days, and the minimum analysis time should reach t = 139.3 days. Thus, if the permeability were reduced 10 times; i.e., for k = 8.64 × 10 −5 m/day, the time step should have been increased by the same order; i.e., ∆t = 139.3 days, and the minimum analysis time should have been 10 times longer (t = 1393 days) to ensure that the pore pressure dissipated to zero. Three equations representing the model calibration curves for the different vertical flow conditions without lateral flow and with lateral restrictions are shown in Figure A1.  Figure A4 shows behavior analogous to that depicted in Figure A3, but in these calibration curves, flow across the lateral boundaries was added and no lateral restrictions were considered. Finally, in the simulations under these conditions, the calibration curves provided equations that allowed for the determination of the required time step ∆t and analysis time t as a function of the permeability coefficient k.