Consolidation by Vertical Drains Considering the Rheological Characteristics of Soil under Depth and Time-Dependent Loading

: Vertical drains have been widely used in soil improvement projects to accelerate consolidation and improve the strength of soft soil. Previous consolidation models are mainly based on instantaneous loading, which cannot reﬂect the actual time-variable loading condition. This paper presents a set of analytical solutions for consolidation with vertical drains under depth and time-dependent loading, involving one-step loading, multi-step loading, and cyclic loading. A four-element model, which is a combination of the Merchant model and Maxwell model, is introduced to consider the rheological characteristics of the soil. By simulating the results of a consolidation test, it is found that the four-element model is more accurate than the Merchant model in predicting the changes in pore pressure and settlement during the clay consolidation process. Based on the solutions obtained, several factors affecting consolidation behavior are investigated. It has been shown that the rheological behavior becomes more and more obvious at the later stage of consolidation with the decrease both in the modulus of the spring in the Kelvin body and the viscosity coefﬁcient of the independent dashpot. With the increase in the viscosity coefﬁcient of the dashpot in the Kelvin body, the rate of consolidation becomes faster at an initial stage but slower at a later stage. For cyclic loading, the consolidation degree in each cycle reaches a maximum at the end of unloading and reaches a minimum at the beginning of the loading.


Introduction
Vertical drains have been widely used to accelerate the consolidation process of soft clayey soils [1][2][3][4][5][6][7].Accurate prediction of the consolidation process is of great significance to prevent geological and environmental disasters from occurring on soft soil foundations (e.g., the deformation and instability of foundation construction).Various consolidation models have been established to guide the construction of vertical drain projects in soft soil foundations, which contributes to the sustainable development of underground space.Since Barron [2] proposed the first solution to investigate radial consolidation for vertical drains, many studies have been performed to extend consolidation theory for vertical drains to reflect more realistic engineering problems in the past decades [2][3][4][5][6][7].For example, based on Barron's consolidation theory, Hansbo et al. [3] developed an analytical solution for consolidation with vertical drains taking into consideration the smear effect and well resistance.Kianfar et al. [6] analyzed the consolidation behavior for vertical drains by incorporating the effects of vacuum preloading and non-Darcian flow.Deng et al. [7] proposed an analytical solution to investigate radial consolidation by vertical drains with the variation of discharge capacity.However, the above studies are based on instantaneous loading.
In most engineering projects, the surcharge loading is time-variable, the pattern of loading includes one-step loading, multi-step loading, cyclic loading, and so on.Therefore, the investigation of consolidation theories with time-dependent loadings is of theoretical and practical significance [8][9][10][11][12].Leo [9] derived a series of closed-form solutions for the consolidation of vertical drains subjected to ramp loading.Conte and Troncone [11] presented an analytical solution for radial consolidation with vertical drains and general time-dependent loading.Lu et al. [12] pointed out that the loading rate and the loading pattern have great effects on the consolidation behavior of soils with vertical drains.In the above investigations, soft clay was treated as a linear elastic material, although the rheological properties of soft soil should be considered.
The rheological properties of soil have important engineering significance.In recent years, many studies of one-dimensional consolidation considering the rheological properties of soil were greatly developed [13][14][15], but there were few studies on consolidation by vertical drains considering the rheological properties of soil.For consolidation theory by vertical drains, Liu et al. [16] employed the Merchant model to study the consolidation theory by vertical drains.Wang and Xie [17] analyzed the consolidation process by vertical drains in the consideration of the rheological properties of soil under a semi-permeable boundary.Huang et al. [18] obtained the analytical solutions for consolidation with vertical drains by incorporating a fractional derivative model.Nevertheless, few of the consolidation theories by vertical drains consider time-dependent loading conditions [16][17][18][19], such as multi-step loading.
Furthermore, the above approaches all assumed that the additional stress caused by external loads was assumed to be uniform with depth.However, many factors may lead to the non-uniform distribution of additional stress along the depth under external loads [14,[20][21][22][23][24][25][26][27][28].Lu et al. [21] obtained the solution for consolidation by vertical drains with a linear change of additional stress with depth and analyzed the influence of stress distribution patterns along the depth on consolidation behavior.Satwik et al. [25] performed a numerical analysis of the consolidation of soft clays subjected to cyclic loading.Chen et al. [27] developed a consolidation model for multilayered unsaturated ground under multistage depth-dependent stress.When the effective stress in soft clay remains unchanged, the deformation of soil does not stop but continues to develop slowly with time (i.e., rheological characteristics).In addition, the loading condition in the actual project often varies with time.If the consolidation model does not consider rheological characteristics and time-dependent loading, it will improperly estimate the consolidation process of soft clay foundations and may result in an uneven settlement or failure of upper structures.At present, there is little research on consolidation theory by vertical drains considering the rheological characteristics and time-dependent stress.
In this paper, the rheological characteristics of soil are considered by using a fourelement model, and the analytical solutions are developed for consolidation by vertical drains under depth and time-dependent loadings.Using the general solutions, detailed solutions are provided for several common constitutive models, such as the Merchant model, Maxwell model, and linear elastic model.Finally, the influence of parameters on consolidation behaviors is discussed carefully based on the four-element model.

Basic Assumptions and Mathematical Modelling
The schematic diagram for the consolidation of saturated clayey soil with a vertical drain under depth and time-dependent loading is shown in Figure 1.As shown in Figure 1, H is the thickness of the soil layer; r and z are the radial and vertical coordinates, respectively; t is the time; and the increase in additional stress σ(z, t) can be assumed to vary arbitrarily with time and linearly with depth, and it is assumed to have the following form: where 0 ≤ q(t) ≤ 1 is a function of time; σ t and σ b are the maximum additional stress at z = 0 and z = H, respectively, and the variation pattern of additional stress along depth can be depicted by the top-to-bottom additional stress ratio (σ t /σ b ). Figure 2 shows the variation of q(t) with t in different loading conditions, which represents the cases of instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.
where 0 ≤ () ≤ 1 is a function of time;  and  are the maximum additional stress at  = 0 and  = , respectively, and the variation pattern of additional stress along depth can be depicted by the top-to-bottom additional stress ratio (  ⁄ ). Figure 2 shows the variation of () with  in different loading conditions, which represents the cases of instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.As shown in Figure 3a, a four-element model is adopted to simulate the rheological characteristics of the clayey soil layer in this paper.An independent spring of modulus  , an independent dashpot of the viscosity coefficient  , a Kelvin body with a spring of modulus  , and another dashpot of the viscosity coefficient  are included in this fourelement model.It can be seen that the model is a combination of the Merchant model (shown in Figure 3b) and the Maxwell model (shown in Figure 3c).Figure 3d shows a linear elastic model that has been widely used in geotechnical engineering.where 0 ≤ () ≤ 1 is a function of time;  and  are the maximum additional stress at  = 0 and  = , respectively, and the variation pattern of additional stress along depth can be depicted by the top-to-bottom additional stress ratio (  ⁄ ). Figure 2 shows the variation of () with  in different loading conditions, which represents the cases of instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.As shown in Figure 3a, a four-element model is adopted to simulate the rheological characteristics of the clayey soil layer in this paper.An independent spring of modulus  , an independent dashpot of the viscosity coefficient  , a Kelvin body with a spring of modulus  , and another dashpot of the viscosity coefficient  are included in this fourelement model.It can be seen that the model is a combination of the Merchant model (shown in Figure 3b) and the Maxwell model (shown in Figure 3c).Figure 3d shows a linear elastic model that has been widely used in geotechnical engineering.As shown in Figure 3a, a four-element model is adopted to simulate the rheological characteristics of the clayey soil layer in this paper.An independent spring of modulus E 0 , an independent dashpot of the viscosity coefficient η 0 , a Kelvin body with a spring of modulus E 1 , and another dashpot of the viscosity coefficient η 1 are included in this four-element model.It can be seen that the model is a combination of the Merchant model (shown in Figure 3b) and the Maxwell model (shown in Figure 3c).Figure 3d shows a linear elastic model that has been widely used in geotechnical engineering.Based on the four-element model, the constitutive relationship of clayey soil can be written as follows: where ε v and σ v are the vertical strain and the vertical effective stress, respectively.Based on Barron's theory of equal strain consolidation, the basic assumptions in this study are listed as follows: (a) Compressive strain occurs only in the vertical direction, and the radial sections of the cell remain radial during the consolidation; that is, the equal vertical strain condition is valid.(b) Both vertical and radial flows obey Darcy's law.(c) The radial flow in the vertical drain is neglected and the radial flow from the soil into the vertical drain at any depth is equal to the corresponding increase in flow up the vertical drain.(d) The constitutive relationship of the soil follows Equation ( 2).(e) The radial coefficient of permeability of the smear zone is smaller than that of undisturbed soil, but the other physical properties of the soil in the smear zone are the same as those in the undisturbed zone.
According to the above assumptions, the equation governing the equal strain consolidation of vertical drains can be expressed as follows: where r w , r s , and r e are the radius of the vertical drain, the radius of the smear zone, and the radius of the undisturbed zone, respectively; k v , k h , and k s are the coefficient of vertical permeability of the soil and the coefficient of radial permeability in the undisturbed zone and in the smear zone, respectively; γ w is the unit weight of water; u s and u n are the excess pore water pressure within the smear zone and undisturbed zone at any point and at any time, respectively; and u is the average excess pore water pressure throughout the entire soil layer at a given depth, which can be obtained from the following expression: According to Terzaghi's principle of effective stress, the relationship of u and σ v can be given by In the vertical flow direction, the top surface of the soil layer can be considered to be pervious, while the bottom of the layer can be considered to be impervious.In the radial flow direction, there is no flow occurring across the external surface of the soil cylinder.Thus, the vertical boundary conditions and the radial external boundary condition for consolidation with vertical drains can be expressed as follows: where u w is the excess pore water pressure in the vertical drain.
According to the continuity condition of radial flow and assumption (c), the excess pore water pressures on both sides of the interface between the vertical drain and the smear zone are equal, and the radial flow from the soil into the vertical drain at any depth is equal to the corresponding increase in flow up the vertical drain.Therefore, the continuity conditions in the interface can be expressed as follows: where k w is the vertical permeability coefficient in the vertical drain.
Combined with the radial boundary conditions, after some mathematical processing, the partial differential equations containing only u and ε v are obtained as follows: where From Equations ( 2) and ( 5), the differential form of ε v with t can be expressed as where and the expression of G(z, t) is given by By substituting Equation (12) into Equation (11), the governing equation for consolidation with vertical drains can be expressed in the following form: where The initial condition for consolidation with vertical drains can be written as follows: As shown in Figure 3a, when it is the case of instantaneous loading, q(t) = q(0) = 1, Equation (14) will become According to the vertical boundary conditions, the solution to Equation ( 16) can be expressed as where M = (2m − 1)π/2, m = 1, 2, 3  6) and (7).If each function T m (t) satisfies Equation ( 16), the solution can be obtained: where (−1) m+1 .To obtain the analytical solution, let L[T m (t)] = F m (s), where L[T m (t)] is the Laplace transformation of T m (t), and applying a Laplace transformation to Equation (18) with the initial condition of Equation ( 15), the following expression can be obtained: Inverting the Laplace transformation of Equation ( 19), the expression of T m (t) can be obtained as follows: where Thus, the expression of the average excess pore water pressure u can be written as follows:

The Overall Average Degree of Consolidation under Instantaneous Loading
After the expression for average excess pore water pressure u is known, the overall average degree of consolidation U can be obtained as follows: As shown in Figure 2b, the function of time q(t) for one-step loading can be written as follows: The solution to Equation ( 14) under the one-step loading condition can be solved by two methods.

(i) Method 1
As the boundary conditions have not changed, the solution to Equation ( 14) can also be expressed as Equation (17).Thus, the following differential equation should be satisfied for each function T m (t): where Q(t) is given by Same as before, L[T m (t)] = F m (s), where L[T m (t)] is the Laplace transformation of T m (t), and applying a Laplace transformation to Equation (32) with the initial condition of Equation ( 15), the following expression can be obtained: where L[T m (t)] = F m (s), and the expression of F m1 (s) is as follows: Applying an inverse Laplace transformation to Equation (35), the following expression can be obtained: where expressions for λ m1 , λ m2 , and a m are given in Equations ( 21), (22) and (26), respectively.Then, applying an inverse Laplace transformation to Equation (34), the following expression can be obtained: Substituting Equations (33), (34), and (36) into Equation (37), expressions for average excess pore water pressure u will be obtained as follows: 1.
(ii) Method 2 According to Xie et al. [13], the expression for average excess pore water pressure can be obtained based on the following equation: where q(t) is the function of time shown in Equation (31), and u(t) and u t are average excess pore water pressures under instantaneous loading and under one-step loading, respectively.Substituting Equations ( 29) and (31) into Equation (42), the solution to Equation ( 14) is same as Equations (39) and (41).Hence, these two methods can both be used for obtaining solutions to Equation ( 14) under the one-step loading condition.

The Overall Degree of Consolidation under One-Step Loading
When the average excess pore water pressure u is known, the overall average degree of consolidation U can be expressed as follows: Substituting Equations ( 1), (31), (39), and (41) into Equation (43), the following expression can be obtained: 1.

Solutions for Multi-Step Loading
In coastal regions, considering that the strength of soil is usually very small, the rapid rate of loading may lead to the instability of the ground.Therefore, multi-step loading, as shown in Figure 3c, is often employed.For this case, the function of time q(t) can be written as follows: where ,q 0 = 0, t 0 = 0; i denotes the ith step of loading, i has the values of 1, 2, . . .i, and q ∞ = 1.

Average Excess Pore Water Pressure under Multi-Step Loading
Substituting Equations ( 29) and ( 46) into Equation ( 42), the solution to Equation ( 14) can be expressed in the following form: 1. For where For where expressions for λ m1 , λ m2 , A m1 , A m2 and C m are given in Equations ( 21)- (25), respectively; and ( * ) is given in Equation (48).

The Overall Degree of Consolidation under Multi-Step Loading
Once the expression for average excess pore water pressure u is found, the overall average degree of consolidation U can be obtained from Equation (43).

Solutions for Cyclic Loading
The function of time q(t) shown in Figure 2d represents the cyclic loading condition, which can be expressed in the following form:

Average Excess Pore Water Pressure under Cyclic Loading
Substituting Equations ( 29) and (52) into Equation (42), the solution to Equation ( 14) can be expressed in the following form: 1. For where where where where where expressions for λ m1 , λ m2 , A m1 , A m2 , and C m are given in Equations ( 21)-( 25), respectively, and ( * ) is given in Equation (55).

The Overall Average Degree of Consolidation under Cyclic Loading
Once the expression for average excess pore water pressure u is found, the overall average degree of consolidation U can be obtained from Equation (43).

Simplification of the Model
It is obvious that the constitutive model shown in Figure 3b-d can be simplified by the four-element model shown in Figure 3a.

Solutions for Different Constitutive Models and Loading Types
According to the above methods, all the solutions under different constitutive models and loading types can be obtained.Since the solutions for the four-element model have been given above under different loading types, only the three other constitutive models as shown in Figure 3b-d are given here.Tables 1-4 show solutions for different constitutive models under instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.

Solutions for Different Constitutive Models and Loading Types
According to the above methods, all the solutions under different constitutive models and loading types can be obtained.Since the solutions for the four-element model have been given above under different loading types, only the three other constitutive models as shown in Figure 3b-d are given here.Tables 1-4 show solutions for different constitutive models under instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and (67), respectively.

Linear elastic model
where expressions for  and  are given in Equations ( 71) and (77), respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively.

Solutions for Different Constitutive Models and Loading Types
According to the above methods, all the solutions under different constitutive models and loading types can be obtained.Since the solutions for the four-element model have been given above under different loading types, only the three other constitutive models as shown in Figure 3b-d are given here.Tables 1-4 show solutions for different constitutive models under instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and (67), respectively.

Linear elastic model
where expressions for  and  are given in Equations ( 71) and (77), respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively.
Maxwell model

Solutions for Different Constitutive Models and Loading Types
According to the above methods, all the solutions under different constitutive models and loading types can be obtained.Since the solutions for the four-element model have been given above under different loading types, only the three other constitutive models as shown in Figure 3b-d are given here.Tables 1-4 show solutions for different constitutive models under instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and (67), respectively.

Linear elastic model
where expressions for  and  are given in Equations ( 71) and (77), respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively.
Maxwell model where expressions for λ m1 and A m1 are given in Equations ( 71) and (77), respectively.

Solutions for Different Constitutive Models and Loading Types
According to the above methods, all the solutions under different constitutive models and loading types can be obtained.Since the solutions for the four-element model have been given above under different loading types, only the three other constitutive models as shown in Figure 3b-d are given here.Tables 1-4 show solutions for different constitutive models under instantaneous loading, one-step loading, multi-step loading, and cyclic loading, respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and (67), respectively.

Linear elastic model
where expressions for  and  are given in Equations ( 71) and (77), respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively; and Maxwell model where expressions for  ,  and  are given in Equations ( 71), ( 73) and (75), respectively; and Linear elastic model where expressions for  and  are given in Equations ( 71) and (77), respectively; and Table 4. Solutions for different constitutive models under cyclic loading.

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), (66) and (67), respectively; and Maxwell model Linear elastic model where expressions for  and  are given in Equations ( 71) and (77), respectively; and where expressions for λ m1 , A m1 and C m are given in Equations ( 71), ( 73) and (75), respectively; and where expressions for λ m1 and A m1 are given in Equations ( 71) and (77), respectively; and Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and ( 67 Maxwell model where expressions for  ,  and  are given in Equations ( 71 where expressions for λ m1 , λ m2 , A m1 and A m2 are given in Equations ( 21), ( 22), ( 66) and (67), respectively; and Merchant model where expressions for  ,  ,  and  are given in Equations ( 21), ( 22), ( 66) and (67), respectively; and Maxwell model where expressions for  ,  and  are given in Equations ( 71 Linear elastic where expressions for λ m1 , A m1 and C m are given in Equations ( 71), (73), and (75), respectively; and

Model Average Excess Pore Water Pressure and Overall Average Degree of Consolidation
Linear elastic model Maxwell model where expressions for  ,  and  are given in Equations ( 71 Linear elastic model where expressions for  and  are given in Equations ( 71) and (77), respectively; and where expressions for λ m1 and A m1 are given in Equations ( 71) and (77), respectively; and   4 shows the comparison between the test results and analytical solutions.It can be clearly seen that the four-element model proposed in this study can accurately predict pore pressure and settlement changes, while the Merchant model will overestimate pore pressure drop and settlement development.For example, the pore pressure drop and settlement calculated by the four-element model at 1000 s are 67 kPa and 0.52 mm, respectively, while the pore pressure drop and settlement calculated by the Merchant model increase to 75 kPa and 0.81 mm, respectively.In addition, the pore pressure remained unchanged after 10,000 s; that is, the effective stress in the soil sample remained unchanged.However, the settlement continued to increase between 10,000 s and 100,000 s, which reflects the rheological characteristics of soft clay.MPa.s, kv = 7.22 × 10 −8 m/s. Figure 4 shows the comparison between the test results and analytical solutions.It can be clearly seen that the four-element model proposed in this study can accurately predict pore pressure and settlement changes, while the Merchant model will overestimate pore pressure drop and settlement development.For example, the pore pressure drop and settlement calculated by the four-element model at 1000 s are 67 kPa and 0.52 mm, respectively, while the pore pressure drop and settlement calculated by the Merchant model increase to 75 kPa and 0.81 mm, respectively.In addition, the pore pressure remained unchanged after 10,000 s; that is, the effective stress in the soil sample remained unchanged.However, the settlement continued to increase between 10,000 s and 100,000 s, which reflects the rheological characteristics of soft clay.

Consolidation Behavior Analysis
To investigate the consolidation behavior of vertical drains considering the rheological characteristics of soil under time and depth-dependent loading, taking the one-step and cyclic loading as examples, several factors affecting the consolidation behavior are discussed.In this section, the reference parameters used for analysis are: r w = 0.07 m; r s = 0.28 m; r e = 0.7 m; H = 10 m;

Consolidation Behavior Analysis under One-Step Loading Condition
For the solutions under the one-step loading condition, the expression of the overall average degree of consolidation U can be seen in Equations ( 44) and (45).In the parameter analysis, t 1 , σ t /σ b , E 1 , η 1 , and η 0 are selected to analyze the consolidation behavior.
Figure 5 shows the influence of loading time t 1 on consolidation behavior, where t 1 varies from 0.1 to 10 days.From Figure 5, it can be seen that the consolidation rate decreases with the increase in the loading time t 1 at the initial stage.When the loading time increased from 0.1 days to 10 days, the consolidation degree decreased from 80% to 55% on the 10th day.However, with the increase in time t, the overall average degree of consolidation U under different loading time t 1 tends to be the same.

Consolidation Behavior Analysis under One-Step Loading Condition
For the solutions under the one-step loading condition, the expression of the overall average degree of consolidation  can be seen in Equations ( 44) and (45).In the parameter analysis,  ,   ⁄ ,  ,  , and  are selected to analyze the consolidation behavior.
Figure 5 shows the influence of loading time  on consolidation behavior, where  varies from 0.1 to 10 days.From Figure 5, it can be seen that the consolidation rate decreases with the increase in the loading time  at the initial stage.When the loading time increased from 0.1 days to 10 days, the consolidation degree decreased from 80% to 55% on the 10th day.However, with the increase in time t, the overall average degree of consolidation  under different loading time  tends to be the same.In Figure 6, the influence of the distribution of additional stress with depth (  ⁄ ) on consolidation behavior is investigated, where   ⁄ = ∞ and   ⁄ = 0 represents that the additional stress has an inverted triangle distribution along the depth and a triangle distribution, respectively.It shows that the overall average degree of consolidation  increases with the increase in   ⁄ and the consolidation rate is the greatest with an inverted triangle distribution of additional stress.As   ⁄ increases from 0 to 2.5, the consolidation degree increases from 72% to 81% on the 10th day.The main reason for this phenomenon is that the vertical pervious boundary is located on the In Figure 6, the influence of the distribution of additional stress with depth (σ t /σ b ) on consolidation behavior is investigated, where σ t /σ b = ∞ and σ t /σ b = 0 represents that the additional stress has an inverted triangle distribution along the depth and a triangle distribution, respectively.It shows that the overall average degree of consolidation U increases with the increase in σ t /σ b and the consolidation rate is the greatest with an inverted triangle distribution of additional stress.As σ t /σ b increases from 0 to 2.5, the consolidation degree increases from 72% to 81% on the 10th day.The main reason for this phenomenon is that the vertical pervious boundary is located on the ground surface, so the larger the additional stress of the top surface, the faster the consolidation rate.In Figure 7,  varies from 0.1 to 100 MPa; thus, it shows the influence of  (i.e. the modules of the spring in the Kelvin body) on consolidation behavior.It can be seen that the consolidation degree curves are almost identical under different  at the initia stage of consolidation.Nevertheless, the rate of consolidation decreases and th rheological behavior becomes more and more obvious with the decrease in  and th increase in time .As  increases from 0.1 to 100 MPa, the steady consolidation degre (i.e., maximum consolidation degree) increases from 90 to 97%. Figure 8 shows the influence of  (i.e., the viscosity coefficient of the dashpot in th Kelvin body) on consolidation behavior, where  varies from 1.0 × 10 to 1.0 × 10 MPa • s.As can be seen from the figure, the consolidation curves intersect with each other.In the initial stage, the consolidation is accelerated with the increase in  ; however the consolidation rate gradually slows down in the later stage.In Figure 7, E 1 varies from 0.1 to 100 MPa; thus, it shows the influence of E 1 (i.e., the modules of the spring in the Kelvin body) on consolidation behavior.It can be seen that the consolidation degree curves are almost identical under different E 1 at the initial stage of consolidation.Nevertheless, the rate of consolidation decreases and the rheological behavior becomes more and more obvious with the decrease in E 1 and the increase in time t.As E 1 increases from 0.1 to 100 MPa, the steady consolidation degree (i.e., maximum consolidation degree) increases from 90 to 97%.In Figure 7,  varies from 0.1 to 100 MPa; thus, it shows the influence of  (i.e the modules of the spring in the Kelvin body) on consolidation behavior.It can be see that the consolidation degree curves are almost identical under different  at the initi stage of consolidation.Nevertheless, the rate of consolidation decreases and th rheological behavior becomes more and more obvious with the decrease in  and th increase in time .As  increases from 0.1 to 100 MPa, the steady consolidation degre (i.e., maximum consolidation degree) increases from 90 to 97%. Figure 8 shows the influence of  (i.e., the viscosity coefficient of the dashpot in th Kelvin body) on consolidation behavior, where  varies from 1.0 × 10 to 1.0 10 MPa • s.As can be seen from the figure, the consolidation curves intersect with eac other.In the initial stage, the consolidation is accelerated with the increase in  ; howeve the consolidation rate gradually slows down in the later stage.Figure 8 shows the influence of η 1 (i.e., the viscosity coefficient of the dashpot in the Kelvin body) on consolidation behavior, where η 1 varies from 1.0 × 10 5 to 1.0 × 10 9 MPa•s.As can be seen from the figure, the consolidation curves intersect with each other.In the initial stage, the consolidation is accelerated with the increase in η 1 ; however, the consolidation rate gradually slows down in the later stage.
Figure 8 shows the influence of  (i.e., the viscosity coefficient of the dashpot in Kelvin body) on consolidation behavior, where  varies from 1.0 × 10 to 1 10 MPa • s.As can be seen from the figure, the consolidation curves intersect with e other.In the initial stage, the consolidation is accelerated with the increase in  ; howe the consolidation rate gradually slows down in the later stage.The influence of η 0 (i.e., the viscosity coefficient of the independent dashpot) on the consolidation behavior can be seen in Figure 9, where η 0 varies from 1.0 × 10 5 to 1.0 × 10 9 MPa•s.It can be seen that the consolidation rate increases with the decrease in η 0 at the initial stage, but the final overall average degree of consolidation decreases with the decrease in η 0 .As η 0 increases from 1 × 10 5 MPa•s to 1 × 10 9 MPa•s, the consolidation degree decreases from 45% to 10% on the 1st day, while the consolidation degree increases from 78 to 100% on the 100th day.The influence of  (i.e., the viscosity coefficient of the independent dashpot) on consolidation behavior can be seen in Figure 9, where  varies from 1.0 × 10 1.0 × 10 MPa • s.It can be seen that the consolidation rate increases with the decrea  at the initial stage, but the final overall average degree of consolidation decreases the decrease in  .As  increases from 1×10 5 MPas to 1×10 9 MPas, the consolida degree decreases from 45% to 10% on the 1st day, while the consolidation degree incre from 78 to 100% on the 100th day.

Consolidation Behavior Analysis under Cyclic Loading Condition
For the solutions under the cyclic loading condition, the expression of the ov average degree of consolidation  can be seen in Equations ( 62)-(65).In this chapter time parameters related to cyclic loads are  = 0.25 and  = 2.0, and the parameter   ⁄ ,  ,  , and  are also used to study consolidation behaviors.
Figure 10 shows the influence of loading time  on consolidation behavior.It sh that the maximum degree of consolidation increases with the increase in  , while minimum degree of consolidation decreases with the increase in  during the c loading.The consolidation degree curves tend to a steady state after a certain numb loading cycles.

Consolidation Behavior Analysis under Cyclic Loading Condition
For the solutions under the cyclic loading condition, the expression of the overall average degree of consolidation U can be seen in Equations ( 62)-(65).In this chapter, the time parameters related to cyclic loads are α = 0.25 and β = 2.0, and the parameters t 1 , σ t /σ b , E 1 , η 1 , and η 0 are also used to study consolidation behaviors.
Figure 10 shows the influence of loading time t 1 on consolidation behavior.It shows that the maximum degree of consolidation increases with the increase in t 1 , while the minimum degree of consolidation decreases with the increase in t 1 during the cyclic loading.The consolidation degree curves tend to a steady state after a certain number of loading cycles.
time parameters related to cyclic loads are  = 0.25 and  = 2.0, and the parameters  ,   ⁄ ,  ,  , and  are also used to study consolidation behaviors.
Figure 10 shows the influence of loading time  on consolidation behavior.It shows that the maximum degree of consolidation increases with the increase in  , while the minimum degree of consolidation decreases with the increase in  during the cyclic loading.The consolidation degree curves tend to a steady state after a certain number of loading cycles.In Figure 11, the influence of   ⁄ on the consolidation behavior is studied.It can be seen that the amount of relative variation of consolidation within each cycle becomes larger for larger   ⁄ , and the amount of relative variation of consolidation degree within each cycle is the smallest with a triangle distribution of additional stress.With In Figure 11, the influence of σ t /σ b on the consolidation behavior is studied.It can be seen that the amount of relative variation of consolidation within each cycle becomes larger for larger σ t /σ b , and the amount of relative variation of consolidation degree within each cycle is the smallest with a triangle distribution of additional stress.With σ t /σ b increasing from 0 to 1, the variation of the consolidation degree under cyclic loading increases from 8 to 15%.  ⁄ increasing from 0 to 1, the variation of the consolidation degree under cyclic loading increases from 8 to 15%.The influence of  on consolidation can be seen in Figure 12.From Figure 12, it can be seen that the consolidation rate is accelerated with the increase in  at the initial stage.However, when the number of cycles increases to a certain value, the consolidation degree curves under different  almost coincide completely.Figure 13 shows the influence of  on consolidation behavior.As can be seen in the figure, the effect of  on the consolidation behavior is similar to that of  , but the effect of  on the consolidation rate is relatively greater.The consolidation rate increases with The influence of E 1 on consolidation can be seen in Figure 12.From Figure 12, it can be seen that the consolidation rate is accelerated with the increase in E 1 at the initial stage.However, when the number of cycles increases to a certain value, the consolidation degree curves under different E 1 almost coincide completely.  ⁄ increasing from 0 to 1, the variation of the consolidation degree under cyclic loading increases from 8 to 15%.The influence of  on consolidation can be seen in Figure 12.From Figure 12, it can be seen that the consolidation rate is accelerated with the increase in  at the initial stage.However, when the number of cycles increases to a certain value, the consolidation degree curves under different  almost coincide completely.Figure 13 shows the influence of  on consolidation behavior.As can be seen in the figure, the effect of  on the consolidation behavior is similar to that of  , but the effect of  on the consolidation rate is relatively greater.The consolidation rate increases with the increase in  at the initial stage, but the influence of  becomes insignificant at the Figure 13 shows the influence of η 1 on consolidation behavior.As can be seen in the figure, the effect of η 1 on the consolidation behavior is similar to that of E 1 , but the effect of η 1 on the consolidation rate is relatively greater.The consolidation rate increases with the increase in η 1 at the initial stage, but the influence of η 1 becomes insignificant at the later stage.Figure 13 shows the influence of  on consolidation behavior.As can be seen in the figure, the effect of  on the consolidation behavior is similar to that of  , but the effect of  on the consolidation rate is relatively greater.The consolidation rate increases with the increase in  at the initial stage, but the influence of  becomes insignificant at the later stage.The influence of η 0 can be seen in Figure 14.It can be seen that the consolidation rate decreases with the increase in η 0 at the initial stage and the amount of relative variation of the consolidation degree within each cycle is larger with smaller η 0 during the entire consolidation process.As η 0 increases from 1 × 10 5 MPa•s to 1 × 10 9 MPa•s, the variation of consolidation degree decreases from 55 to 13% under cyclic loading.The influence of  can be seen in Figure 14.It can be seen that the consolidation rate decreases with the increase in  at the initial stage and the amount of relative variation of the consolidation degree within each cycle is larger with smaller  during the entire consolidation process.As  increases from 1 × 10 5 MPas to 1 × 10 9 MPas, the variation of consolidation degree decreases from 55 to 13% under cyclic loading.From these figures, it should be noted that the consolidation degree in each cycle reaches a maximum at the end of unloading and the consolidation degree is the minimum at the beginning of the loading.

Discussion
For the one-dimensional consolidation problem of soft soil foundation under timedependent loading, most of the previous studies obtained numerical solutions [25][26][27][28], while the relevant analytical solutions were few and complicated [24].In this study, a general solution for the one-dimensional consolidation of a soft soil foundation was derived.Compared with existing analytical solutions, the variable substitution method was used to obtain the governing equation, which is simple in form and easy to solve under time-dependent loading.Then, the general solution of one-dimensional consolidation was obtained by using the separation variable method.The general solution contains only one integral of the loading function, which is easy to operate in complex loading conditions.Finally, this study gave the analytical solutions of the one-dimensional consolidation of soft soil foundation under several common loading types and analyzed the consolidation behavior under one-step loading and cyclic loading.These analytical solutions provide theoretical guidance for determining construction parameters such as the loading mode, loading rate, and loading time in the preloading method for soft soil From these figures, it should be noted that the consolidation degree in each cycle reaches a maximum at the end of unloading and the consolidation degree is the minimum at the beginning of the loading.

Discussion
For the one-dimensional consolidation problem of soft soil foundation under timedependent loading, most of the previous studies obtained numerical solutions [25][26][27][28], while the relevant analytical solutions were few and complicated [24].In this study, a general solution for the one-dimensional consolidation of a soft soil foundation was derived.Compared with existing analytical solutions, the variable substitution method was used to obtain the governing equation, which is simple in form and easy to solve under time-dependent loading.Then, the general solution of one-dimensional consolidation was obtained by using the separation variable method.The general solution contains only one integral of the loading function, which is easy to operate in complex loading conditions.Finally, this study gave the analytical solutions of the one-dimensional consolidation of soft soil foundation under several common loading types and analyzed the consolidation behavior under one-step loading and cyclic loading.These analytical solutions provide theoretical guidance for determining construction parameters such as the loading mode, loading rate, and loading time in the preloading method for soft soil foundation treatment.This study also summarizes the change rule of consolidation degree with time under different loading parameters.
In summary, this study further expands and enriches the analytical theory of the onedimensional consolidation of soil and provides theoretical guidance for the construction process of soft soil foundations under variable loading.

Conclusions
The following conclusions may be drawn from this study: 1.
Based on Barron's theory of equal strain consolidation, a four-element model was used to consider the rheological characteristics of soil, and a set of analytical solutions was developed for consolidation with vertical drains under depth and time-dependent loading.The increase in additional stress is a function depending on both time and depth.It is assumed to vary linearly with depth, and several time functions are considered to represent different loading cases, which include instantaneous loading, one-step loading, multi-step loading, and cyclic loading.

2.
The consolidation rate is accelerated with the decrease in loading time t 1 and the increase in σ t /σ b (the value of the top-to-bottom additional stress ratio).With the decrease both of the modulus of the spring in the Kelvin body and the viscosity coefficient of the independent dashpot, the rheological behavior becomes more and more obvious at the later stage of consolidation.The rate of consolidation becomes faster at an early stage but slower at a later stage, with the increase in the viscosity coefficient of the dashpot in the Kelvin body.

3.
For cyclic loading, the consolidation degree in each cycle reaches a maximum at the end of unloading and reaches the minimum at the beginning of the loading.When the number of cycles increases to a certain value, the variation form of consolidation degree curves will tend to be stable.

Figure 1 .
Figure 1.Schematic diagram for consolidation with a vertical drain.

Figure 1 .
Figure 1.Schematic diagram for consolidation with a vertical drain.

Figure 1 .
Figure 1.Schematic diagram for consolidation with a vertical drain.

Figure 3 .Figure 3 .
Figure 3. Model types: (a) four-element model; (b) Merchant model; (c) Maxwell model; and (d) linear elastic model.Based on the four-element model, the constitutive relationship of clayey soil can be written as follows:

Sustainability 2023 ,
15, x FOR PEER REVIEW 15 of 22  () = [∑ ( * ) ] + Model Validation Li et al. (2013) [29] adopted GDS advanced consolidation test system to test the consolidation process of Xiaoshan soft clay.The diameter and height of the soil sample were 76.2 mm and 20 mm, respectively.The water content, density, and void ratio of the soil sample were 50.22%, 1.68 g/cm 3 , and 1.43, respectively.The consolidation pressure was 200 kPa.The proposed four-element model and the traditional three-element model (i.e., the Merchant model) were applied to simulate the test results.The values of model parameters were set as: E0 = 3.11 MPa, E1 = 4.56 MPa, η0 = 8.51 × 10 5 MPa.s, η1 = 2.80 × 10 4

Figure 5 .
Figure 5.The influence of  on the consolidation behavior under one-step loading.

Figure 5 .
Figure 5.The influence of t 1 on the consolidation behavior under one-step loading.

Sustainability 2023 , 2 Figure 6 .
Figure 6.The influence of   ⁄ on the consolidation behavior under one-step loading.

Figure 7 .
Figure 7.The influence of  on the consolidation behavior under one-step loading.

Sustainability 2023 ,Figure 6 .
Figure 6.The influence of   ⁄ on the consolidation behavior under one-step loading.

Figure 7 .
Figure 7.The influence of  on the consolidation behavior under one-step loading.

Figure 7 .
Figure 7.The influence of E 1 on the consolidation behavior under one-step loading.

Figure 8 .
Figure 8.The influence of η 1 on the consolidation behavior under one-step loading.

Sustainability 2023 , 18 Figure 8 .
Figure 8.The influence of  on the consolidation behavior under one-step loading.

Figure 9 .
Figure 9.The influence of  on the consolidation behavior under one-step loading.

Figure 9 .
Figure 9.The influence of η 0 on the consolidation behavior under one-step loading.

Figure 10 .
Figure 10.The influence of  on the consolidation behavior under cyclic loading.

1 Figure 10 .
Figure 10.The influence of t 1 on the consolidation behavior under cyclic loading.

Figure 11 .
Figure 11.The influence of   ⁄ on the consolidation behavior under cyclic loading.

Figure 12 .
Figure12.The influence of  on the consolidation behavior under cyclic loading.

UFigure 11 .
Figure 11.The influence of σ t /σ b on the consolidation behavior under cyclic loading.

Figure 11 .
Figure 11.The influence of   ⁄ on the consolidation behavior under cyclic loading.

Figure 12 .
Figure 12.The influence of  on the consolidation behavior under cyclic loading.

UFigure 12 .
Figure 12.The influence of E 1 on the consolidation behavior under cyclic loading.

Figure 12 .
Figure 12.The influence of  on the consolidation behavior under cyclic loading.

Figure 13 .
Figure 13.The influence of  on the consolidation behavior under cyclic loading.

Figure 13 .
Figure 13.The influence of η 1 on the consolidation behavior under cyclic loading.

Figure 14 .
Figure 14.The influence of  on the consolidation behavior under cyclic loading.

t 1 =Figure 14 .
Figure 14.The influence of η 0 on the consolidation behavior under cyclic loading.
• • • , and T m (t) is a function of t.It is obvious that Equation (17) satisfies the boundary conditions of Equations ( e λ m1 t + A m2 e λ m2 t + C m (30)

Table 1 .
Solutions for different constitutive models under instantaneous loading.

Table 1 .
Solutions for different constitutive models under instantaneous loading.

Table 2 .
Solutions for different constitutive models under one-step loading.

Table 1 .
Solutions for different constitutive models under instantaneous loading.

Table 2 .
Solutions for different constitutive models under one-step loading.
e λ m1 t + C m where expressions for λ m1 , A m1 and C m are given in

Table 1 .
Solutions for different constitutive models under instantaneous loading.

Table 2 .
Solutions for different constitutive models under one-step loading.

Table 2 .
Solutions for different constitutive models under one-step loading.

Table 1 .
Solutions for different constitutive models under instantaneous loading.

Table 2 .
Solutions for different constitutive models under one-step loading.

Table 3 .
Solutions for different constitutive models under multi-step loading.

Table 3 .
Solutions for different constitutive models under multi-step loading.

Table 3 .
Solutions for different constitutive models under multi-step loading.

Table 3 .
Solutions for different constitutive models under multi-step loading.

Table 4 .
Solutions for different constitutive models under cyclic loading.

Table 4 .
Solutions for different constitutive models under cyclic loading.