A Novel Approach to the Analysis of the Soil Consolidation Problem by Using Non-Classical Rheological Schemes

: The paper presents classical and non-classical rheological schemes used to formulate constitutive models of the one-dimensional consolidation problem. The authors paid special attention to the secondary consolidation effects in organic soils as well as the soil over-consolidation phenomenon. The systems of partial differential equations were formulated for every model and solved numerically to obtain settlement curves. Selected numerical results were compared with standard oedometer laboratory test data carried out by the authors on organic soil samples. Additionally, plasticity phenomenon and non-classical rheological elements were included in order to take into account soil over-consolidation behaviour in the one-dimensional settlement model. A new way of formulating constitutive equations for the soil skeleton and predicting the relationship between the effective stress and strain or void ratio was presented. Rheological structures provide a ﬂexible tool for creating complex constitutive relationships of soil.


Introduction
Soil is a complex porous three-phase material in which many phenomena take place simultaneously. Because of the complexity a lot of soil mechanics problems have to be solved numerically or by laboratory investigation [1]. For numerical calculations, mainly using the finite element method (FEM), appropriate material models have to be implemented. As there are many kinds of soil types, differing in particle size fractions, geological origin and the conditions of loading, there is no universal solution. Many commercial FEM software packages provide their own material models [2]. However, they also allow for extending the models as there is often a need to analyse some less typical soil mechanics problems or soil types more thoroughly. In case of porous media, constitutive equations are generally formulated in the effective stress and applied for coupled pore pressure-stress numerical calculations of geotechnical problems [3].
Consolidation and settlement calculations are very important aspects of geotechnical design. The settlement of soil can be divided into three parts: (i) Elastic response, (ii) primary consolidation settlement related to the drainage of water, (iii) secondary settlement, related to the creep of the soil skeleton. In most cases only the primary consolidation is considered as it causes the largest settlement for the majority of soils. However, for organic or soft clay soils, the secondary compression can be also large and has to be taken into account [4]. Consolidation of the soft soil layer is also important in the case of column supported embankments, where the load transferred to the column and geosynthetic layer strongly depends on the settlement of the soft layer [5,6].
Rheological schemes provide a flexible tool for creating complicated constitutive relationships of various materials. By applying rheological schemes to model the behaviour of the soil skeleton and coupling them with the equation of pore water flow, it is possible to describe the transient (time-dependant) response of soil taking into account secondary consolidation of viscoelastic soils [7][8][9][10][11][12]. In the study, a method of modelling the transient process of primary and secondary consolidation of soil based on Burgers rheological scheme is presented. The Burgers structure is commonly used to model the behaviour of asphalt mixtures [13,14]. Here, it will be applied for the soil skeleton in order to simulate the transient process of consolidation of an organic soil layer. It should be noted that the presented procedure can be applied for any viscoelastic rheological structures. This shows that the great potential of constitutive models developed in different research fields can be straightforwardly utilized to solve geotechnical problems. An example of such a problem is the modelling of secondary consolidation. The 1D constitutive equations derived using rheological schemes can be implemented as a part of the evolution law of the Cap yield surface of the Modified Drucker-Prager/Cap model available in commercial FEM codes.
On the other hand, the behaviour of soil in the effective stress space (e − log σ curve) can be modelled by using non-classical rheological schemes as will be demonstrated in this paper. To do this, the non-classical Kepes element [15] will be introduced in a rheological scheme of the soil skeleton in effective stress, and based on that scheme a system of algebraic-differential equations will be formulated. This system of equations has not been presented previously in the geotechnical literature. Moreover, it should be emphasised that the Kepes element has not been applied to the geotechnical problem either. A detailed description of the Kepes rheological model for both 1D and 3D cases can be found in the monograph by Zbiciak [16].

One-Dimensional Consolidation Model
The one-dimensional consolidation problem is the consolidation of an infinite saturated homogeneous clay layer, loaded by an infinite and uniform pressure. The infinite layer of soil undergoing consolidation can be drained on both sides (see Figure 1) or only on one side and rests on a non-deformable base (e.g., rock). Only excess pore water pressure, which arises after application of the external load, is considered. Despite the fact that one-dimensional formulation seems to be a simplistic approach it is often used for the prediction of the settlement of engineering structures, such as road embankments [17][18][19].
The simplest model of one-dimensional consolidation was presented by Karl Terzaghi [20]. The partial differential equation of flow can be expressed as: where k is the coefficient of permeability, γ w is the unit weight of water, σ is the effective stress function of depth and time, and ε is strain function of depth and time.
Assuming the linear constitutive law for the soil's skeleton, in which the strain function ε(z, t) is related to the effective stress σ (z, t) by the following equation where M denotes the oedometer modulus, Equation (1) can be solved with appropriate initial and boundary conditions. For the two-side drainage (see Figure 1) the boundary conditions are as follows where q(t) is a function of the applied load and H is the thickness of the soil layer. The initial condition is σ (z, 0) = 0. Equation (1) can be solved both analytically and numerically [20]. In the paper [10] numerical solutions of the Terzaghi problem were demonstrated. Moreover, some more advanced rheological schemes modelling the soil skeleton were also taken into consideration and the results of numerical solutions were compared with known analytical formulas [8].
The simple linear relationship between the effective stress and strain (Equation (2)) does not take into account secondary consolidation to be observed in organic or soft clay soils. This can be shown by comparison of the model prediction with oedometer laboratory test data carried out on a remoulded sample of soil with about 7% organic matter content (see Figure 2). Taking different values of parameters M and k, it was impossible to obtain a satisfying fit to the oedometer test results.

Classical Rheological Schemes for Consolidation
Various classical rheological structures, represented by connections of springs and dashpots, are widely used to model the rate-dependent behaviour of asphaltic materials (see [21][22][23][24][25]) and polymers [26]. Such schemes for modelling problems of soil mechanics were adopted, e.g., in papers [7,8]. Gibson and Lo [8] presented an analytical solution of one-dimensional consolidation in which the relationship between the effective stress and vertical strain was formulated according to the standard model.

Burgers Model for Soil Skeleton
In this section, we introduce another rheological scheme not widely used in soil mechanics-the Burgers model (see Figure 3). This model is commonly used to model the behaviour of asphalt mixtures [13,14]. Here, it will be applied for the soil skeleton to simulate the transient soil response during consolidation and to obtain the settlement change in time. Assuming the Burgers model for the soil skeleton (see Figure 3), the constitutive relationship can be formulated as follows [10]: where the strain function ε v (z, t) can be determined by solving the following differential equation and the superimposed dot denotes the derivative with respect to time. In Equations (4) and (5) the rheological structure's parameters (see Figure 3) M 1 and M 2 can be interpreted as the oedometric (constrained) moduli and η 1 and η 2 are the viscosity constants.
Equations (4) and (5) are an extension of the linear relationship (2). Thus, in order to simulate the process of transient consolidation, the system of coupled differential equations (1), (4) and (5) has to be solved with the boundary conditions (3) and 'zero' initial conditions for ε(z, t) and ε v (z, t). After solving the initial boundary value problem, the function of displacement (settlement) can be found by integrating the strain function.
Advanced elastic-plastic models implemented in various geotechnical finite element method software packages can be extended by taking into account soil's skeleton creep. For example the Drucker-Prager/Cap model available in Abaqus softwere [27] gives the possibility to implement the constitutive relationship between the creep volumetric strain and effective pressure (first invariant of the effective stress tensor). Thus, the onedimensional considerations shown in this section can be applied to the existing threedimensional models. In order to implement such extension, in Equations (4) and (5), 1D strain, ε, should be replaced by volumetric strain, ε vol = tr ε ε ε, and σ by effective pressure p = − 1 3 tr σ σ σ .

Numerical Solution
The system of partial differential equations was solved using the numerical method of lines implemented in Mathematica software [10,28]. The assumed parameters of the Burgers model are shown in Table 1. The comparison with the example test data is shown in Figure 4. Comparing Figures 2 and 4 shows that, with this approach, the oedometer test results can be predicted a lot better than they were in the case of the simple elastic relationship (2). In order to achieve even better fit of the model, the parameters of the model can be found by an automatic curve fitting to the laboratory test data. The procedure and results will be shown in future works.

Kepes Element for Effective Stress Behaviour
In Section 3 transient response of soil was analysed, that is the change of settlement in time related to the drainage of pore water and rheological effects in the soil skeleton. Models of transient consolidation can be used to predict settlement of a structure after a certain amount of time, for example, in order to calculate how much settlement occurs during the time of construction. In this section models for different conditions will be presented. Here, we introduce a non-classical rheological element to model the behaviour of soil in effective stress-that is only the response of the soil skeleton in drained conditions. Such models of steady-state consolidation can be used for calculation of long-term (maximum) settlement which occurs when there is no excess pore pressure.

Formulation
Kepes model describes materials which dissipate energy. It is a special case of a plastic material, in which the set of the admissible stress depends on the strain state. The process of energy dissipation for the Kepes material as well as the plastic material is rate independent. The phenomenon of rate independency makes Kepes and plastic elements different from the classic viscous element (dashpot). The assumption is made that for an undeformed body stresses are zero. The one-dimensional constitutive relationship for the Kepes model is similar to the relationship of a classical slider element which describes ideal plasticity but with the plastic limit not constant but dependent on the strain: σ 0 (ε). For the one-dimensional case the continuous set of admissible stress for the Kepes material is as follows [16] Θ kp := [−σ 0 ; σ 0 ] where σ 0 = κ|ε| and κ is a material constant. The Kepes element can be effectively used to model the behaviour of soil in effective stress as observed in oedometer or isotropic consolidation tests. The typical e − log σ curve obtained from these tests is shown in Figure  5. While we are going to formulate the relationship between vertical strain (equivalent to the volumetric strain for the assumed conditions) and stress, the void ratio e can be calculated, for a given soil's initial void ratio e 0 , using the well-known formula: The first model we will analyse is a connection in series of the Kepes element and a spring (see Figure 6a). To formulate the constitutive relationship in effective stress for such a scheme, let us first denote the logarithm of the effective stress as In the definition (8) natural logarithm can also be used. The choice depends on the way results of oedometer or isotropic consolidation tests are presented. Here, as the common logarithm is used, the slopes of the virgin compression line and the unloadingreloading line will be related to the compression index (C c ) and to the swelling index (C s ), respectively. Now, for the Kepes element we can define σ L = βσ 0 where σ 0 = κ|ε p |, so: and the derivative with respect to time is: .
For elements connected in a series, as shown in Figure 6a, the stress in the spring is where M is the oedometric (constrained) modulus. Thus, the strain in the Kepes element is hence .
Substituting Equation (10) to Equation (13) we get After some algebra, Equation (14) can be rewritten as It can be shown that assuming only compression and no tension in the model (which is the case for soils), Equation (15) can be reduced to the following form [16]: where the projection onto the set of negative numbers was used: Numerical implementation needs the case of ε p = 0 to be treated in a special manner to avoid indeterminate formulas (see Equation (17)). Only for this case can an appropriate spring element be substituted instead of the Kepes element, so where Equation (17) can be solved for a given strain excitation function ε(t). After determining ε p (t), the logarithmic stress can be calculated using Equation (11). Then, linear effective stress σ = 10 σ L (compare Equation (8)).
The model presented in Figure 6a is rate independent, as was stated in the introduction of Section 4. Thus, the results of simulation will not depend on the velocity of applied strain excitation. In order to take into consideration the effect of rate dependency, the scheme presented in Figure 6a can be modified by adding an additional element, for example, a dashpot. This way, secondary consolidation phenomena can be taken into account.
The extended model is presented in Figure 6b. After introducing a dashpot connected in parallel, Equation (11) becomes: where η is viscosity. For the viscosity parameter η = 0 model presented in Figure 6a is obtained.

Numerical Tests
The equations were tested for the following function of strain ε(t) = At sin(ωt) + Bt (22) in which A = 3 × 10 −7 1 min and ω = 0.025 1 min were assumed. In order to test the model's prediction for different strain rates, two values of the parameter B were assumed: B = 1 × 10 −5 1 min and B = 1 × 10 −6 1 min . The functions of the strain excitation are shown in Figure 7. The model was tested for parameters shown in Table 2. The idea was to find a function which can approximate a consolidation test and use it to verify the applicability of the model.   Figure 8 shows the relation between log σ and the vertical strain obtained by the model for η = 0, that is for the rate-independent scheme (see Figure 6a). The produced curve is the same as the idealised graph obtained in consolidation laboratory tests. The relation between vertical effective stress and strain is shown in Figure 9. Corresponding curves after adding rate-dependency, by setting the viscosity parameter η = 20 MPa · min, are shown in Figures 10 and 11. In Figure 10 one can see that the normal consolidation line shifts to the right-hand side for the lower strain rate. This phenomenon is confirmed by many researchers for soil exhibiting secondary compression [29].  In order to take into account pre-consolidated soils, a simple modification of Equation (16) should be done by introducing an additional term in the denominator as follows: where σ p is the preconsolidation pressure. The simulation results of the preconsolidated soil model, assuming σ p = 50 kPa, are visualized in Figure 12.

Discussion
FEM software with an open structure like, for example, Abaqus gives the possibility to implement parts of constitutive models (e.g., evolution of the yield surface) by applying user-defined subroutines. In particular, in case of the Modified Drucker-Prager/Cap model, the evolution of the Cap yield surface is defined in terms of the volumetric inelastic strain, that being the sum of plastic volumetric strain and creep volumetric strain (ε cr vol ). The ε cr vol is related to modelling the secondary consolidation phenomenon. The consolidation creep mechanism in the Modified Drucker-Prager/Cap is defined by specifying a constitutive law for ε cr vol in a form of differential equation in terms of an equivalent creep pressure and ε cr vol . Three laws for creep are implemented in Abaqus, one of them being the Singh-Mitchell law. However, any equation for calculating ε cr vol can be given by a user subroutine CREEP. Thus, the equations presented in the manuscript for transient consolidation can be implemented as an extension to the Modified Drucker-Prager/Cap model. Rheological structures serve as a useful tool for building creep laws. In the article, the Burgers structure was analysed in detail; however, even more complicated viscoelastic models can be used in order to relatively easily find equations describing creep.
The steady-state consolidation model presented in the manuscript, in which the Kepes element was used, can be also utilized in order to build and to implement a complete material model for consolidation analysis. Geotechnical FEM software (i.e., Plaxis) allows user definitions of material models programmed in Fortran. The authors plan to extend the presented approach to a three-dimensional case and implement the constitutive model using the UMAT user subroutine available in Abaqus software.
The right-hand branch of the steady-state consolidation model can be further modified to obtain a more suitable solution. This will be done in future works. Moreover, both incremental loading (IL) and constant rate of strain (CRS) oedometer tests will be carried out by the authors in order to further validate the proposed model.
The equation of flow can be modified to simulate the flow of water more realistically [4]. In a transient FEM analysis the equations of flow are solved numerically. Moreover, according to test data, the coefficient of permeability changes in time as the consolidation progresses and the density of the material increases [30]. The authors plan to assume a non-constant function of this parameter in the equation of flow in future works.

Conclusions
Transient one-dimensional consolidation of soil was modelled by using the Burgers rheological scheme and numerical integration of differential equations. Comparison with laboratory test data shows that the model provides good prediction of the soil's behaviour and an iterative curve fitting procedure can be implemented for finding the model's parameters based on laboratory tests. Constitutive equations of the soil skeleton, formulated here for the one-dimensional transient consolidation case, can be also applied for the isotropic consolidation conditions. Additionally, a new way of formulating constitutive equations for the soil skeleton and predicting the relationship between the effective stress and strain in drained conditions (e − log σ curve) was presented. The set of equations for the non-classical Kepes rheological element connected in series with a spring was presented. The model was then extended by a simple addition of a rate-dependent dashpot element which allowed to take into account secondary consolidation effects. Testing the model for two different strain rates showed the well-known phenomena of the normal consolidation line shift. The presented solution can be treated as an extension or alternative to the Isotache model being developed by many researchers [29,31].
The main findings of this paper can be summarized in the following points: • A method of modelling the transient process of primary and secondary consolidation of soil based on rheological schemes was presented. In the paper Burgers structure was used as an example; however, the presented procedure can be applied for any viscoelastic rhelogical structures. The 1D constitutive equations derived using rheological schemes can be implemented as part of the evolution law of the Cap yield surface in commercial finite element method codes. • A new way of formulating constitutive equations for the soil skeleton in steady-state (drained) conditions using the non-classical Kepes element was presented. • An explicit form of differential equations was formulated for the non-classical Kepes rheological element. • Presented rheological models are able to describe the material response for cyclic loading conditions: Both for loading and unloading.
• Formulated equations can be solved using known algorithms for solving ordinary differential equations (e.g., the Runge-Kutta method). Funding: This research received no external funding.

Institutional Review Board Statement:
This study did not require ethical approval.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data supporting the findings of this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest.