## 3 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_{\rm 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]:  
\[\sigma(z,t)=\Big{(}1+\zeta\frac{z}{H}\Big{)}q(t) \tag{8}\]  
where \(\sigma(z,t)\) is total vertical stress in soil layer; \(\zeta\) is a dimensionless parameter controlling the gradient of total vertical stress in soil layer. Figure 1: One-dimensional consolidation model of the viscoelastic soils. On the basis of the hypothesis of 1D consolidation theory [9], the differential equation governing the consolidation process of saturated soils is given by  
\[\frac{k_{\mathrm{v}}}{\gamma_{\mathrm{w}}}\frac{\partial^{2}u(z,t)}{\partial z^ {2}}=-\frac{\partial\epsilon(z,t)}{\partial t} \tag{9}\]  
where \(\epsilon(z,t)\) and \(u(z,t)\) are, respectively, the vertical strain and excess pore-water pressure in the soil layer; and \(\gamma_{\mathrm{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  
\[u(z,t)\big{|}_{z=0}=0,\quad\frac{\partial u(z,t)}{\partial z}\Big{|}_{z=H}=0 \tag{10}\]  
At the initial time, excess pore-water pressure is equal to total vertical stress, thus the initial condition of Equation (9) is given by  
\[u(z,0)=\sigma(z,0)=\Big{(}1+\zeta\frac{z}{H}\Big{)}q(0) \tag{11}\]  
According to the principle of effective stress [9], the vertical effective stress \(\sigma^{\prime}(z,t)\) in soil layer is given by  
\[\sigma^{\prime}(z,t)=\sigma(z,t)-u(z,t) \tag{12}\]  
By substituting Equation (11) into Equation (12) results in  
\[\sigma^{\prime}(z,0)=0 \tag{13}\]  
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]. The relationship between the stress and strain of independent spring is  
\[\sigma_{\mathrm{e}}=E_{0}\epsilon_{\mathrm{e}} \tag{14}\]  
where \(\sigma_{\mathrm{e}},\epsilon_{\mathrm{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  
\[\sigma_{\mathrm{v}}=\eta_{0}\frac{d^{\alpha}\sigma_{\mathrm{e}v}}{d^{\alpha }\eta_{0}}=\frac{E_{0}}{1-\alpha_{0}}\int_{0}^{t}\frac{d\epsilon_{\mathrm{e}v }}{d\tau}e^{-\frac{\alpha_{0}}{1-\alpha_{0}}(\frac{t-\tau}{\lambda_{1}})}d \tau,\quad(0\leq\alpha_{0}\leq 1) \tag{15}\]  
where \(\sigma_{\mathrm{v}},\epsilon_{\mathrm{v}}\) and \(\eta_{0}\) stand for the stress, strain and viscosity coefficient of the independent fractional-derivative element, respectively; and \(\lambda_{0}=\eta_{0}/E_{0}\) is its creep time. The relationship between the stress and strain of fractional-derivative Kelvin-Voigt body is  
\[\sigma(z,t)=\Big{(}1+\zeta\frac{z}{H}\Big{)}q(t) \tag{8}\]0  
where \(\sigma_{\mathrm{e}v}\), \(\epsilon_{\mathrm{e}v}\), \(E_{1}\) and \(\eta_{1}\) are the stress, strain, modulus and viscosity coefficient of the fractional-derivative Kelvin-Voigt body, respectively; and \(\lambda_{1}=\eta_{1}/E_{1}\) is its creep time.