One-Dimensional Consolidation of Viscoelastic Soils Incorporating Caputo-Fabrizio Fractional Derivative

: In this paper, the Caputo-Fabrizio fractional derivative is introduced to investigate the one-dimensional consolidation behavior of viscoelastic soils. Using the Caputo-Fabrizio operator, a novel four-element fractional-derivative model is proposed to capture the viscoelastic properties of the soils, and further the one-dimensional consolidation equation is derived to simulate the consolidation behavior of the soils. Using the techniques of eigenfunction expansion and Laplace transform, a series of analytical solutions are derived to calculate the excess pore-water pressure and the average degree of consolidation of the soils. The total vertical stress in the soil is assumed to change linearly with depth, and its distribution patterns are classiﬁed to rectangular pattern, trapezoidal pattern and inverse trapezoidal pattern. Four loading types including instantaneous loading, ramp loading, sinusoidal loading and general cyclic loading are considered. Then, a comparison for several special cases is presented to verify the correctness of the proposed solutions through comparing with existing theories. Moreover, two examples considering ramp and sinusoidal loadings are given to study the consolidation behavior of the viscoelastic soils incorporating the Caputo-Fabrizio fractional derivative.


Introduction
As a branch of mathematics, the fractional calculus deals with the generalization of integrals and derivatives to all real (and even complex) orders [1,2]. People have discovered that the fractional calculus can provide more accurate models than the integer calculus to describe the hereditary and memory characteristics of various complex processes and natural phenomena [3]. During the last two decades, the fractional calculus has been widely applied to model the complex processes in the fields of science, engineering and mathematics [1,[4][5][6][7], and so on. In particular, since the fractional calculus was firstly used to describe the viscoelasticity and rheology in 1980s, the fractional-derivative model has achieved great success in simulating the viscoelastic behavior of real materials [1,8]. However, the application of the fractional-derivative to the one-dimensional consolidation of the viscoelastic soils is scarce.
In the geotechnical community, the soil consolidation is one of the most common issues [9][10][11][12]. After the seminal work of Terzaghi [9], a considerable number of onedimensional consolidation models have been proposed to simulate the consolidation process of the soil [13][14][15][16][17]. For obtaining a simplified solution, the one-dimensional consolidation models put forward in the early days treated the soil as a linear elastic material.
In reality, the soil generally shows obvious rheological properties during the consolidation process. Afterwards, in order to capture the rheological properties of the soil, some viscoelastic models including Maxwell model, Kelvin-Voigt model, Merchant model, four-element model and generalized Kelvin-Voigt model have been introduced to extend the one-dimensional consolidation model of the soil [13,[18][19][20][21]. These consolidation models simulated the rheological properties of the soil using the viscoelastic models with integer derivatives (namely standard viscoelastic model). Recently, the concept of fractional calculus has been introduced to develop the viscoelastic model for describing the rheological phenomena of materials. Gemant [22] put forward a series of fractional-derivative models to study the rheological phenomena. Koeller [2] derived the creep and relaxation functions of the fractional-derivative models. In the geotechnical community, Yin et al. [23,24] proposed a fractional-derivative model to describe the viscoelasticity of the soil samples under triaxial conditions. Zhu et al. [8] adopted a fractional-derivative model to predict the compression properties of the soft clay. Inspired by the above work, Wang et al. [25] introduced a fractional-derivative model to extend the one-dimensional consolidation theory of viscoelastic soil. However, the available fractional-derivative models are usually defined by the Caputo or Riemann-Liouville fractional derivatives with a kernel of power function [26][27][28]. The kernel of power function is investigated to be singular at the end point of the interval [26].
To solve the singularity of the kernel of power function, Caputo and Fabrizio recently put forward a new fractional derivative using the kernel of exponential function [27,28]. This fractional derivative (namely Caputo-Fabrizio fractional derivative) has no singularity and it has been proved to be suitable for some classical physical problems. Using the Caputo-Fabrizio fractional derivative, this paper attempts to establish one-dimensional consolidation model of viscoelastic soils. In order to capture the rheological properties of the soils, the Caputo-Fabrizio operator is adopted to develop a novel four-element fractional-derivative model, and one-dimensional consolidation model is establish to simulate the complex consolidation process of the viscoelastic soils. Using the techniques of eigenfunction expansion and Laplace transform, a series of analytical solutions are derived to calculate the excess pore-water pressure and the average degree of consolidation of the viscoelastic soils. Then, a comparison for several special cases is presented to verify the correctness of the proposed solutions, and two examples are given to investigate the consolidation behavior of the viscoelastic soils.

Fundamental Concept
The Riemann-Liouville and Caputo fractional derivatives, respectively, are defined as [1] d α f (t) where f (t) is a given function; Γ(t) is the Gamma function; α is the fractional order; t is the time; and τ is a integral variable. The main problem of these two derivatives is that the kernel of power function is singular at the end point (i.e., τ = t) of the interval. To avoid this problem, Caputo and Fabrizio recently proposed a new fractional derivative using a kernel of exponential function [27]. This derivative (namely Caputo-Fabrizio fractional derivative) is defined as [27,28] 3 of 16 where ω(α) is a normalized function and it can be reasonably given as ω(α) = 1, for detail see [28]. According to this definition, it is obvious that the new kernel has no singularity. For the case of α → 0 , Caputo and Fabrizio [27] state that Equation (3) reduces to On the other hand, for the case of α → 1 , Equation in addition lim where δ(t) is the Dirac delta function.
Using the Caputo-Fabrizio operator, a fractional-derivative element can be defined as where σ(t), ε(t), E and η stand for the applied stress, corresponding strain, elastic modulus and viscosity coefficient of the fractional-derivative element, respectively; and λ = η/E is a creep time. It is clearly observed that this element will reduce to a linear elastic spring as α = 0, whereas becomes an ideal Newtonian dashpot as α = 1.

Model Formulation
A schematic of the consolidation system consisting of one soil layer is drawn in Figure 1. The thickness of the soil layer is H, and the coefficient of vertical permeability of the soil is k v . The coordinate origin is at the top of the soil layer, and the z-axis is perpendicular to the surface of the soil layer. The time-dependent loading q(t) is applied within a limited area of the surface or uniformly distributed on the whole surface. At any time, total vertical stress in soil layer is assumed to change linearly with depth, and its expression is supposed as [15,17]: where σ(z, t) is total vertical stress in soil layer; ζ is a dimensionless parameter controlling the gradient of total vertical stress in soil layer.
At the initial time, excess pore-water pressure is equal to total vertical stress, thus the initial condition of Equation (9) is given by According to the principle of effective stress [9], the vertical effective stress ( ) in soil layer is given by By substituting Equation (11) into Equation (12) results in ( ) As shown in Figure 1, a novel four-element fractional-derivative model is developed to describe the rheological characteristics of the viscoelastic soils by employing the fractional-derivative element given by Equation (7). This model consists of an independent spring, an independent fractional-derivative element and a fractional-derivative Kelvin-Voigt body, which is made of a fractional-derivative element in parallel connection with a spring. Herein, the fractional-derivative element is represented by a diamond, which has been adopted by many researchers [8,23,24]. On the basis of the hypothesis of 1D consolidation theory [9], the differential equation governing the consolidation process of saturated soils is given by where ε(z, t) and u(z, t) are, respectively, the vertical strain and excess pore-water pressure in the soil layer; and γ w is the unit weight of water. It is assumed that the top surface of the soil layer is fully permeable, while its bottom base is impermeable. Thus the boundary conditions of Equation (9) can be obtained as At the initial time, excess pore-water pressure is equal to total vertical stress, thus the initial condition of Equation (9) is given by According to the principle of effective stress [9], the vertical effective stress σ (z, t) in soil layer is given by By substituting Equation (11) into Equation (12) results in As shown in Figure 1, a novel four-element fractional-derivative model is developed to describe the rheological characteristics of the viscoelastic soils by employing the fractionalderivative element given by Equation (7). This model consists of an independent spring, an independent fractional-derivative element and a fractional-derivative Kelvin-Voigt body, which is made of a fractional-derivative element in parallel connection with a spring. Herein, the fractional-derivative element is represented by a diamond, which has been adopted by many researchers [8,23,24].
The relationship between the stress and strain of independent spring is where σ e , ε e and E 0 are the stress, strain and modulus of the independent spring, respectively. The relationship between the stress and strain of independent fractional-derivative element is where σ v , ε v and η 0 stand for the stress, strain and viscosity coefficient of the independent fractional-derivative element, respectively; and λ 0 = η 0 /E 0 is its creep time.
The relationship between the stress and strain of fractional-derivative Kelvin-Voigt body is where σ ev , ε ev , E 1 and η 1 are the stress, strain, modulus and viscosity coefficient of the fractional-derivative Kelvin-Voigt body, respectively; and λ 1 = η 1 /E 1 is its creep time.
The vertical strain ε z of the soil mass is the sum of the strain components of the four-element fractional-derivative model; that is, The vertical effective stress σ z of the soil mass is equal to the stress components of the four-element fractional-derivative model; that is, Based on Equations (13)- (18), the initial values of all the strain components in the fourelement fractional-derivative model are zero. The Laplace transforms of Equations (14)-(18) produce where L is a Laplace transform operator. The inverse Laplace transform of Equation (19) leads to Equation (20) gives the stress-strain relationship of the four-element fractionalderivative model. It can be seen that the expression form of Equation (20) is mainly governed by the fractional orders α 0 and α 1 . For the conditions that the fractional orders are equal to zero or one, as shown in Table 1, the four-element fractional-derivative model will degenerate into several standard models.

Solution Derivation
The solution to Equation (21), based on eigenfunction expansion, can be expre as  Table 1. Simplifications of the four-element fractional-derivative model.

Solution Derivation
The solution to Equation (21), based on eigenfunction expansion, can be expr as T t is the generalized Fourier coeff that varies with time t . Equation (23) satisfies the boundary conditions given in E tion (10). By substituting Equation (23) into Equation (21), a family of ordinary tegro-differential equations can be obtained as Table 1. Simplifications of the four-element fractional-derivative model.

Solution Derivation
The solution to Equation (21), based on eigenfunction expansion, can be expre as By substituting Equation (23) into Equation (21), a family of ordinary tegro-differential equations can be obtained as By substituting Equations (8), (12) and (20) into Equation (9) leads to the governing equation (21) where c v = k v E 0 /γ w is the vertical coefficient of consolidation; κ = E 1 /E 0 is a modulus ratio; and f q (t) is given by
By substituting Equation (23) into Equation (21), a family of ordinary integro-differential equations can be obtained as Considering the initial condition given in Equation (11), the Laplace transform of Equation (24) yields The inverse Laplace transform into Equation (27) leads to the general solution of Equation (24). The result is where Substituting Equation (28) into Equation (23) produces the general solution of excess pore-water pressure in convolution form. The result is where φ(t) is a dimensionless loading function given by and The average degree of consolidation U p (t), which is defined by excess pore-water pressure, can be calculated from the following formula: Substituting Equations (8) and (33) into Equation (35) yields For instantaneous loading, φ(t) can be expressed as By substituting Equation (37) into Equations (33) and (36), the solutions for excess pore-water pressures and average degree of consolidation can be obtained as follows For ramp loading, φ(t) can be expressed as where t 0 is the construction time of ramp loading. By substituting Equation (40) into Equations (33) and (36), the solution for excess pore-water pressure and average degree of consolidation can be obtained as the following expressions: For sinusoidal loading, φ(t) shown in Figure 2 can be expressed as where T is the period of sinusoidal loading. By substituting Equation (43) into Equations (33) and (36), the solutions for the excess pore-water pressures and average degree of consolidation can be obtained as the following expressions: Using the method of generalized Fourier series, ( ) t φ for a general cyclic loading shown in Figure 3, can be expressed as The Fourier coefficients are given by By substituting Equation (43) into Equations (33) and (36), the solutions for the excess pore-water pressures and average degree of consolidation can be obtained as the following expressions: Using the method of generalized Fourier series, φ(t) for a general cyclic loading shown in Figure 3, can be expressed as By substituting Equation (43) into Equations (33) and (36), the solutions for the excess pore-water pressures and average degree of consolidation can be obtained as the following expressions: Using the method of generalized Fourier series, ( ) t φ for a general cyclic loading shown in Figure 3, can be expressed as The Fourier coefficients are given by By substituting Equation (48) into Equations (33) and (36), the solutions for the excess pore-water pressures and average degree of consolidation can be obtained as the following expressions: The Fourier coefficients are given by a n = 2/T T 0 φ(t) cos(2nπt/T)dt, n = 0, 1, 2 · · · b n = 2/T T 0 φ(t) sin(2nπt/T)dt, n = 1, 2 · · · (48) By substituting Equation (48) into Equations (33) and (36), the solutions for the excess pore-water pressures and average degree of consolidation can be obtained as the following expressions: where

Verification against Some Special Cases
To verify the proposed solution, the following case study is used. The height of soil layer, coefficient of vertical permeability of soil mass and unite weight of water are assumed as follows: H = 6 m, γ w = 10 kN/m 3 , k v = 1 × 10 −7 m/s. The values of the rheological parameter and fractional order for each viscoelastic model are given in Table 2. Considering that external loading is applied instantly and total vertical stress in soil layer is uniform with depth, Figure 4 illustrates the results of excess pore-water pressure obtained from the proposed solutions and existing solutions of Xie et al. [13], Wang et al. [20] and Terzaghi [9]. Figure 4a gives the distribution of excess pore-water pressure with depth, whereas Figure 4b gives the change of excess pore-water pressure with time factor. It can be obviously observed that the results calculated from the proposed solution show good agreement with those from the existing solutions derived from standard four-element model [13], standard Merchant model [20], standard Maxwell model [13] and linear elastic model [9], respectively, for Case I, II, III and IV. Therefore, the results of the proposed solution and some existing theories are identical. Moreover, the proposed solution is more general to predict the consolidation behavior of viscoelastic soils.

Discussion about the Consolidation Behavior
This part presents two examples to discuss the consolidation behavior of viscoelastic soils incorporating the Caputo-Fabrizio fractional derivative. One example is conducted for the ramp loading, and the other one is performed for the sinusoidal loading. In each example, the average degree of consolidation incorporating various parameters is obtained from the proposed solution, and the results are presented as a function of the time

Discussion about the Consolidation Behavior
This part presents two examples to discuss the consolidation behavior of viscoelastic soils incorporating the Caputo-Fabrizio fractional derivative. One example is conducted for the ramp loading, and the other one is performed for the sinusoidal loading. In each example, the average degree of consolidation incorporating various parameters is obtained from the proposed solution, and the results are presented as a function of the time factor.
For the case of ramp loading, the proposed solution for average degree of consolidation, which is given in Equation (42), has seven dimensionless parameters governing the consolidation process. These parameters can be classified into three categories, namely fractional orders (i.e., α 0 and α 1 ), loading parameters (i.e., ζ and T v0 ), and rheological parameters (i.e., κ, ϑ and T λ1 ). α 0 and α 1 reflect the influence of fractional orders of independent fractional-derivative element and that in fractional-derivative Kelvin-Voigt body, which further represent the fluid characteristic of soil mass. Indeed, higher fractional orders indicate that fluid characteristic of soil mass is more obvious. ζ reflects the influence of distribution pattern of total vertical stress with depth. In practical engineering, total vertical stress is supposed to change linearly with depth according to three types of distribution patterns, i.e., rectangular pattern as ζ = 0, trapezoidal pattern as ζ > 0 and inverse trapezoidal pattern as ζ < 0. T v0 (=c v t 0 /H 2 ) reflects the influence of construction time t 0 . When T v0 = 0, ramp loading will become instantaneous loading, and the solution for this condition is given in Equation (39). κ (=E 1 /E 0 ) reflects the influence of modulus ratio between fractional-derivative Kelvin-Voigt body and the independent spring. ϑ (=η 1 /η 0 ) reflects the influence of viscosity coefficient ratio between fractional-derivative Kelvin-Voigt body and the independent fractional-derivative element. T λ1 (=c v λ 1 /H 2 ) reflects the influence of viscosity coefficient of fractional-derivative Kelvin-Voigt body. By changing the values of these parameters one by one, consolidation behavior of viscoelastic soil under ramp loading is investigated in the following discussion. Figure 5 illustrates the influence of fractional orders on consolidation behavior under ramp loading. In Figure 5a, α 0 changes from 0 to 1 while the other parameters remain constant, and the changes of average degree of consolidation (U p ) with time factor are presented. It can be clearly seen that U p -curves under different α 0 intersect each other with the increase in time factor T v . Higher α 0 delivers larger U p in the early stage, but smaller one at a later stage. This indicates that for higher α 0 , excess pore-water pressure dissipates faster in the beginning and then more slowly during consolidation process. Moreover, the ultimate value of U p gradually decreases with the increase in α 0 . In Figure 5b, α 1 changes from 0 to 1, and the other parameters remain constant. In comparison with the results obtained under various α 0 , the results obtained under different α 1 also share similar changing features for U p -curves. The only difference is that for various α 1 , the ultimate value of U p is the same with each other. These influences demonstrate that fractional orders i.e., fluid characteristic of soil mass has a significant impact on consolidation behavior. More fluid soils consolidate initially rapidly but more slowly at a later stage, and even their ultimate degree of consolidation will become smaller, if fluid characteristic is embodied in the independent fractional rheological element.
comparison with the results obtained under various 0 α , the results obtained under different 1 α also share similar changing features for p U -curves. The only difference is that for various 1 α , the ultimate value of p U is the same with each other. These influences demonstrate that fractional orders i.e., fluid characteristic of soil mass has a significant impact on consolidation behavior. More fluid soils consolidate initially rapidly but more slowly at a later stage, and even their ultimate degree of consolidation will become smaller, if fluid characteristic is embodied in the independent fractional rheological element.    Figure 6 shows the influence of loading parameters on consolidation behavior under ramp loading. In Figure 6a, ζ changes from −0.5 to 0.5, and the other parameters remain constant. Herein three U p -curves, respectively, are corresponding to three types of distribution patterns of total vertical stress; that is, inverse trapezoidal pattern (ζ < 0), rectangular pattern (ζ = 0) and trapezoidal pattern (ζ > 0). It can be clearly observed that the distribution pattern of total vertical stress has an important effect on the consolidation behavior. Inverse trapezoidal pattern brings about the most rapid dissipation of excess pore-water pressure; trapezoidal pattern results in the slowest dissipation of excess pore-water pressure; and under rectangular pattern, the dissipation of excess pore-water pressure falls in between them. During the consolidation process, the differences between three U p -curves initially increase and reach a maximum at a given time, and then gradually decrease and finally disappear with the increase in time. Figure 6b illustrates the changes of U p with time factor under various construction times. The construction time factor (i.e., T v0 ) varies from 0 to 1, whereas the other parameters remain constant. It can be obviously observed that larger construction time slow down the dissipation of excess pore-water pressure (i.e., reduce U p ) in the early stage, and then the difference between U p -curves of various construction times decrease gradually as time increases.   varies from 0.1 to 10 one by one, whereas the other parameters remain constant. As observed in Figure 7a, the higher the modulus ratio (i.e., modulus of Kelvin-Voigt body is larger), the faster the dissipation of excess pore-water pressure. Figure 7b indicates that higher viscosity coefficient ratio (i.e., viscosity coefficient of independent fractional-derivative element is smaller) brings about slower dissipation of excess pore-water pressure, and decrease the ultimate consolidation degree as well. Figure 7c shows that p U -curves under different λ1 T intersect each other with the increase in time  Figure 7 illustrates the influence of rheological parameters on consolidation behavior under ramp loading. The dimensionless parameters (i.e., κ, ϑ and T λ1 ) varies from 0.1 to 10 one by one, whereas the other parameters remain constant. As observed in Figure 7a, the higher the modulus ratio (i.e., modulus of Kelvin-Voigt body is larger), the faster the dissipation of excess pore-water pressure. Figure 7b indicates that higher viscosity coefficient ratio (i.e., viscosity coefficient of independent fractional-derivative element is smaller) brings about slower dissipation of excess pore-water pressure, and decrease the ultimate consolidation degree as well. Figure 7c shows that U p -curves under different T λ1 intersect each other with the increase in time factor T v . Higher T λ1 (i.e., the viscosity coefficient of Kelvin-Voigt body is larger) delivers faster dissipation of excess pore-water pressure in the early stage, but more slowly one at a later stage. v0 Figure 7 illustrates the influence of rheological parameters on consolidation behavior under ramp loading. The dimensionless parameters (i.e., κ , ϑ and 1 T λ ) varies from 0.1 to 10 one by one, whereas the other parameters remain constant. As observed in Figure 7a, the higher the modulus ratio (i.e., modulus of Kelvin-Voigt body is larger), the faster the dissipation of excess pore-water pressure. Figure 7b indicates that higher viscosity coefficient ratio (i.e., viscosity coefficient of independent fractional-derivative element is smaller) brings about slower dissipation of excess pore-water pressure, and decrease the ultimate consolidation degree as well. Figure 7c shows that p

U -curves under different λ1
T intersect each other with the increase in time factor v T . Higher λ1 T (i.e., the viscosity coefficient of Kelvin-Voigt body is larger) delivers faster dissipation of excess pore-water pressure in the early stage, but more slowly one at a later stage.   For the case of sinusoidal loading, Equation (45) for average degree of consolidation also has seven governing parameters. All of them are the same as those for the case of ramp loading except the loading parameter T T (=c v T/H 2 ), which reflects the influence of the period (T) of sinusoidal loading on consolidation behavior. By changing the values of these parameters individually while holding the others constant, the consolidation behavior of fractional viscoelastic soils under sinusoidal loading is investigated in the following discussion. Figures 8-10 illustrate the changes of average degree of consolidation (U p ) with the increase in time factor (T v ). It can be observed that the change of U p under sinusoidal loading is a vibrant process. In each period, U p maximizes at the beginning of unloading stage, whereas minimizes at the beginning of the next loading stage. Moreover, the maximum and minimum of U p gradually become stable with the increase in time.
fractional orders result in faster dissipation of excess pore-water pressure in the early stage, but slower one at a later stage. Within each period, the relative variation of p U becomes larger with the increase in fractional orders. Moreover, increasing 0 α (i.e., fractional order of independent fractional rheological element) gradually decreases the ultimate consolidation degree. However, 1 α (i.e., fractional order of Kelvin-Voigt body) has no influence on the ultimate consolidation degree.   Figure 10a, higher modulus ratio results in faster dissipation of excess pore-water pressure in the early stage of consolidation process, and the influence of modulus ratio on the change of p U gradually reduces with the increase in time. Figure 10b demonstrates that higher viscosity coefficient ratio delivers slower dissipation of excess pore-water pressure as well as smaller ultimate consolidation degree. This influence characteristic is the same as the result for the case of ramp loading.. Figure 10c shows that higher λ1 T brings about faster dissipation of excess pore-water pressure in the beginning of consolidation process, and then the difference of p U -curves gradually decreases with the increase in time.   tially go up to a maximum and then gradually decrease as the time increases. Figure 9b shows that higher period (i.e, lower frequency) delivers a larger relative variation of p U within each period, but the period has no influence on the stable value of p U .  T λ individually varying from 0.1 to 10 are taken for the investigations. As observed in Figure 10a, higher modulus ratio results in faster dissipation of excess pore-water pressure in the early stage of consolidation process, and the influence of modulus ratio on the change of p U gradually reduces with the increase in time. Figure 10b demonstrates that higher viscosity coefficient ratio delivers slower dissipation of excess pore-water pressure as well as smaller ultimate consolidation degree. This influence characteristic is the same as the result for the case of ramp loading.. Figure 10c shows that higher λ1 T brings about faster dissipation of excess pore-water pressure in the beginning of consolidation process, and then the difference of p U -curves gradually decreases with the increase in time. T λ1 =1.0  Figure 8 illustrates the influence of fractional orders on consolidation behavior under sinusoidal loading, where α 0 and α 1 individually vary from 0 to 1, whereas the other parameters remain constant. As observed in Figure 8a,b, U p curves under various fractional orders intersect each other with the increase in time factor T v . Higher fractional orders result in faster dissipation of excess pore-water pressure in the early stage, but slower one at a later stage. Within each period, the relative variation of U p becomes larger with the increase in fractional orders. Moreover, increasing α 0 (i.e., fractional order of independent fractional rheological element) gradually decreases the ultimate consolidation degree. However, α 1 (i.e., fractional order of Kelvin-Voigt body) has no influence on the ultimate consolidation degree. Figure 9 shows the influence of loading parameters on consolidation behavior under sinusoidal loading. ζ varying from −0.5 to 0.5 and T T varying from 0.05 to 0.20, respectively, are taken for the investigations in Figure 9a,b. As observed in Figure 9a, the distribution pattern of total vertical stress has a noteworthy effect on the change of U p and its relative variation within each period. Inverse trapezoidal pattern (ζ < 0) delivers the fastest dissipation of excess pore-water pressure as well as the largest relative variation of U p within each period; trapezoidal pattern (ζ > 0) brings about the slowest dissipation of excess pore-water pressure as well as the smallest relative variation of U p within each period; and under rectangular pattern (ζ = 0), the dissipation of excess pore-water pressure and the relative variation of U p fall in between those of two aforementioned patterns. During the consolidation process, the differences between three U p -curves initially go up to a maximum and then gradually decrease as the time increases. Figure 9b shows that higher period (i.e, lower frequency) delivers a larger relative variation of U p within each period, but the period has no influence on the stable value of U p . Figure 10 illustrates the influence of rheological parameters on consolidation behavior under sinusoidal loading. The conditions that κ, ϑ and T λ1 individually varying from 0.1 to 10 are taken for the investigations. As observed in Figure 10a, higher modulus ratio results in faster dissipation of excess pore-water pressure in the early stage of consolidation process, and the influence of modulus ratio on the change of U p gradually reduces with the increase in time. Figure 10b demonstrates that higher viscosity coefficient ratio delivers slower dissipation of excess pore-water pressure as well as smaller ultimate consolidation degree. This influence characteristic is the same as the result for the case of ramp loading. Figure 10c shows that higher T λ1 brings about faster dissipation of excess pore-water pressure in the beginning of consolidation process, and then the difference of U p -curves gradually decreases with the increase in time.

Conclusions
In order to capture the viscoelastic behavior of the soils, this paper has introduced the Caputo-Fabrizio fractional derivative to propose a novel four-element fractional-derivative model, and one-dimensional consolidation model has been further established to predict the consolidation behavior of viscoelastic soils. A series of analytical solutions have been derived to calculate the excess pore-water pressure and the average degree of consolidation. A comparison for several special cases has been given to verify the correctness of the proposed solutions, and two examples have been presented to study the consolidation behavior of the viscoelastic soils. The results show that the proposed solutions agree well with the existing theories. When the fractional orders are equal to zero and/or one, the four-element fractional-derivative model will reduce to several standard models, and particularly it allows a great convenience for practical application because of its generality. The fractional orders, loading parameters and rheological parameters have significant influences on the consolidation behavior. Higher fractional order or viscosity coefficient of Kelvin-Voigt body delivers faster dissipation excess pore-water pressure in the early stage, but slower one at a later stage. Inverse trapezoidal pattern of total vertical stress brings about the fastest dissipation of excess pore-water pressure, and trapezoidal pattern results in the slowest one. The consolidation is accelerated with the increase in modulus ration, while slowed down with the increase in viscosity coefficient ratio as well as construction time. The relative variation of consolidation degree under sinusoidal loading becomes larger as the period increases.