Consolidation of Viscoelastic Soil With Vertical Drains for Continuous Drainage Boundary Conditions Incorporating a Fractional Derivative Model

In geotechnical engineering, vertical drainage is the most economical method for accelerating the consolidation of a large area of soft ground. In this study, we analyze the viscoelasticity of the soil and the actual drainage conditions on the top surface of the soil, and then we introduce continuous drainage boundary conditions and adopt a fractional derivative model to describe the viscoelasticity of the soil. With the use of a viscoelasticity model, the governing partial differential equation for vertical drains under continuous drainage boundary conditions is obtained. With the application of the Crump numerical inversion method, the consolidation solution for vertical drains is also obtained. Further, the rationality of the proposed solution is verified by several examples. Moreover, some examples are provided to discuss the influence of interface drainage parameters on the top surface of soil and the viscoelasticity parameters of soil on the consolidation behavior of vertical drains. The proposed method can be applied in the fields of transport engineering to predict the consolidation settlement of a foundation reinforced by vertical drains.

In geotechnical engineering, vertical drainage is the most economical method for accelerating the consolidation of a large area of soft ground. In this study, we analyze the viscoelasticity of the soil and the actual drainage conditions on the top surface of the soil, and then we introduce continuous drainage boundary conditions and adopt a fractional derivative model to describe the viscoelasticity of the soil. With the use of a viscoelasticity model, the governing partial differential equation for vertical drains under continuous drainage boundary conditions is obtained. With the application of the Crump numerical inversion method, the consolidation solution for vertical drains is also obtained. Further, the rationality of the proposed solution is verified by several examples. Moreover, some examples are provided to discuss the influence of interface drainage parameters on the top surface of soil and the viscoelasticity parameters of soil on the consolidation behavior of vertical drains. The proposed method can be applied in the fields of transport engineering to predict the consolidation settlement of a foundation reinforced by vertical drains.

INTRODUCTION
Consolidation of vertical drains is the most economical and effective method for reinforcement treatment of a large area of soft soil (Indraratna et al., 2010(Indraratna et al., , 2012. In order to theoretically study the consolidation of vertical drains, Carrillo (1942) first proposed the division of a combination of seepage problems in vertical drains into radial seepage and vertical seepage for separate analysis. Subsequently, Barron (1942), Horne (1964), Yoshikuni and Nakanodo (1974), Yoshikuni (1979), and Xie et al. (1993) divided the vertical drain consolidation problem into two ideal conditions, free strain and equal strain, to analyze and answer the problems separately. They established a linear elastic consolidation theory for vertical drains. However, a large amount of scientific research and engineering measurement data has shown that soil has viscoelasticity properties and that its consolidation process cannot be completed in a short time. The soil consolidation process often takes years or even decades to complete (Karlsson et al., 2016).
For the purpose of analyzing viscoelastic soil, many scholars have introduced viscoelasticity models into the consolidation theory of vertical drains, such as the classic Kelvin, Maxwell, Merchant and the generalized Kelvin. Zhao (1988) and Liu et al. (1998Liu et al. ( , 2005 established the consolidation differential equation for integer-order viscoelastic vertical drains and obtained numerous excellent results. In 1940, Blair (1947) and Gerasimov (1948) proposed a fractional derivative model to describe viscoelasticity behaviors of a linear elasticity model and an integer-order viscoelasticity model. Torvik (1986, 2012) verified the rationality of a fractional derivative rheological model by combining a molecular theory with a rheological theory based on fractional derivatives. Subsequently, the fractional derivative viscoelasticity model began to appear in geotechnical engineering. Sun and Wei (2007) and Li-jun et al. (2011) studied the creep of soft soil with a fractional calculus model. Wang et al. (2017) used the fractional derivative model to describe the 1D consolidation of viscoelastic soil. Qing-zi et al. (2016) experimentally verified that the application of the fractional derivative model is effective in describing creep deformation of soft soils. Huang and Li (2019) studied the consolidation characteristics of fractional derivative viscoelastic vertical drains for variable loads based on the fractional-derivative Merchant (FM) model. Because of the studies mentioned above, the fractional derivative model will be used more frequently in the study of vertical drains consolidation in the future.
Traditional boundary conditions are mainly Terzaghi (1943) drainage (permeab le and impervious). Gray (1945) proposed a semi-permeable boundary. Schiffmann and Stein (1970) studied the effect of permeability coefficient on the 1D consolidation of soil. However, the Terzaghi drainage boundary has nothing to do with time, and the boundary conditions contradict the initial conditions. Although the semi-permeable boundary reflects the change process in the permeable boundary, its boundary condition is a function that complicates the solution process. Based on this, Mei et al. (2013) proposed a continuous drainage boundary that was simple in form and could reflect changes in drainage boundary pore pressure with time to describe the actual boundary drainage during soil consolidation.
This study is conducted based on the research results mentioned above and on the FM model. We have introduced the continuous drainage boundary conditions proposed by Mei and have derived a solution to the fractional derivative viscoelastic vertical drain consolidation problem for continuous drainage boundary conditions. Finally, some special cases are given to study the effect of interface drainage parameters and soil fractional viscoelasticity parameters on consolidation properties.

Continuous Drainage Boundary
According to Mei et al. (2013), the expression of the continuous drainage boundary is where q 0 is the initial load and b is the interface drainage parameter, which can be obtained through experiments. t is the time. From Equation (1), it can be easily found that when t → 0, the pore-water pressure on the top surface of soil results in u (0, 0) = q 0 , satisfying the initial condition. When t → ∞, the porewater pressure on the top surface of soil approaches 0. Therefore, Equation (1) solves the initial condition problem, and it reflects the change in pore pressure at the boundary of soil with time.

Consolidation Model of Vertical Drains
In order to facilitate mathematical derivation, the following basic assumptions are adopted in this study, as shown below: (1) Darcy's law applies to seepage of water in soil.
(2) The amount of inflowing water from soil at any depth is equivalent to the amount of outflowing water in vertical drains. (3) External load is infinite and evenly distributed.
(4) Equal strain condition is established; that is, the soil and the vertical drains at the same depth only undergo vertical deformation, and the vertical deformation is equal. (5) The compressibility of smeared area is the same as the compressibility of undisturbed soil. (6) Permeability coefficients in the vertical direction of the smeared area are equal to those of the undisturbed soil, and the horizontal permeability coefficient of the smeared area is less than that of the undisturbed soil.
The consolidation system is composed of vertical drains, smeared area, and natural soil, as shown in Figure 1. From the figure, q 0 is the pre-compression load on the top surface of the soil. r w , r s , and r e are vertical drains radius, smear area radius, and influence radius, respectively. r and z are the radial and vertical coordinates, respectively. k v is the vertical permeability coefficient of the soil in the smeared area and the natural area. k s and k h are the horizontal permeability coefficient of the smeared soil and the natural soil, respectively. k w is the vertical permeability of the vertical drains. γ w is the average severity of water. ε v is the strain of soil in vertical direction. u s (r, z, t) is the excess pore-water pressure in the smeared area. u n (r, z, t) is the excess pore-water pressure in the natural soil.ū (z, t) is the average pore-water pressure in the soil. u w (z, t) is the pore-water pressure in the vertical drains. h is the thickness of the soil. Based on the above assumptions, the basic equation for vertical drain consolidation is as follows (Barron, 1942;Tang and Onitsuka, 2000): The average pore-water pressure of the soil can be expressed as follows: The top boundary in this model is considered to continuous drainage and the bottom boundary is considered to be undrained. At the same time, according to symmetry, no seepage occurs at the boundary of the soil-affected area. Thus, boundary conditions are as follows:ū where b 1 and b 2 represent the interface drainage parameters on the top surface of the soil outside the vertical drains and the interface drainage parameters on the top surface of the vertical drains, respectively. Initial conditions of the consolidation model are as follows:ū

Fractional Derivative Model
In order to determine the rheological consolidation properties of vertical drains more accurately, the theory of vertical drain consolidation adopted the FM model (Huang and Li, 2019), as shown in Figure 2.
The constitutive equations for the FM model can be expressed as where σ is the total stress, ε is the total strain, E 1 is the elastic modulus of spring 1, E 2 is the elastic modulus of spring 2, λ=η/E 2 is the viscosity coefficient of the bomb element, ε 1 and ε 2 are strain of spring 1 and strain of spring 2, respectively, and RL D t α is the fractional-derivative operator of type α of function f (t), which is defined as is the gamma function, and τ is an integral variable. The Laplace transform for Equation (9) can be obtained as follows: This study introduces the FM model to describe the effective stress-strain relationship in the soil. Using the effective stress principle, the effective stress of soil σ ′ v under loading is Obviously, the initial conditions of the soil are By substituting Equation (12) into Equation (11) and performing the inverse Laplace transform on Equation (11), the expression changes as follows: where ⊗ denotes convolution. Furthermore, where δ (t) is the Dirac function. When α = 1, the FM model degenerates to a standard Merchant (SM) model. When α = 0, the FM model degenerates to a linear elasticity model.

Governing Equations and Solutions
The consolidation differential Equation of the fractional viscoelastic vertical drains for continuous drainage boundary conditions is the same as the derivation Equation of (Huang and Li, 2019): where c v = k v E 1 /γ w , c h = k h E 1 /γ w , and n(=r e /r w ) represents the ratio of vertical drains radius to influence radius, and E α,α (t) is a Mittag-Leffler function.
The following dimensionless parameters and variables are defined as Substituting the above parameters and variables into Equations (16) and (17) results in The continuous drainage initial conditions and boundary conditions can be rewritten as Since the continuous drainage boundary conditions are nonhomogeneous, the finite Fourier sine transform is very suitable for solving such a problem. The expression of the finite Fourier sine transform of f (x) is (Sneddon, 1995).
where 0 ≤ z ≤ l. The inverse transformation is defined as The above finite Fourier sine transform has the following basic properties: In order to facilitate the application of finite Fourier sine transform, we applied the principle of mirroring and found that the consolidation problem under the conditions of the single drainage boundary with soil thickness of h is equivalent to the problem in a double sided continuous drainage boundary with soil thickness of 2h. The problem in the double-sided continuous drainage boundary is shown in Figure 3.
The bottom boundary conditions of the equivalent consolidation model are the same as top boundary conditions, which areū Frontiers in Materials | www.frontiersin.org where c 1 = h 2 b 1 /c v and c 2 = h 2 b 2 /c v are the dimensionless interface drainage parameters on the top surface of the foundation soil outside the vertical drains and the dimensionless interface drainage parameters on the top surface of the vertical drains, respectively. Substituting Equations (27) and (28) Implementing the Laplace transform of T v for Equation (18) and using the given initial condition of Equation (29) generate Correspondingly, the boundary conditions of Equations (27) and (28) can be transformed intõ Performing the finite Fourier sine transform onṽ w ξ , p results inṼ where M = mπ/2m = 1, 2, · · · . The inverse transformation of Equation (34) is The basic properties of the finite Fourier sine transform are The finite Fourier sine transform is performed on both ends of Equation (31), and Equations (32) are used to obtain the following: where D m = 8G(n 2 − 1)/ M 2 n 2 . By solvingṼ m in Equation (37), the expression ofṼ m can be obtained by, Performing the inverse Laplace transform onṼ m results in Performing the inverse finite Fourier sine transform on Equation (40) leads to the solution for pore-water pressure in the vertical drains, which is where T m is the inverse Laplace transform ofT m .
where j = √ −1, and β is a real number. Generally, for the case of 0 < α < 1, Equation (44) is difficult to have an analytical solution. The Laplace transform numerical inversion Crump method (Crump, 1976;Koeller, 1984) is used to obtain the solution for T m , which is where T is a real parameter exceeding one-half of the longest time. By substituting Equation (41) into Equation (19), the average pore-water pressure of the soil layer can be obtained as Based on Equation (44), the average degree of consolidation U can be obtained as

Degeneration to the Terzaghi Drainage Boundary
In this section, the rationality of the consolidation solution derived from this study is verified by comparison with two existing theories. The first case involves the degeneration of the solution in this study to the consolidation solution for the Terzaghi drainage boundary condition (i.e., c 1 → ∞, c 2 → ∞). After degeneration, Equation (42) can be simplified tõ Based on Equation (45), the following expression can be obtained: Obviously, the degradation solution in this study is exactly the same as that in the result obtained by Huang and Li (2019) based on the fractional derivative viscoelastic vertical drain consolidation model for the completely pervious boundary condition. As shown in Figure 4A, the calculation result of the solution in Huang and Li (2019) is consistent with the calculation result of the solution in this study.

Degeneration to SM Model
In the second case, the solution in this study is degenerated into the solution with the SM model for the completely pervious boundary (i.e., c 1 → ∞, c 2 → ∞, α= 1) and compared with the consolidation solution of Liu et al. (1998), which is based on the integer-order viscoelasticity model for vertical drains. In order to facilitate the comparison, it is assumed that the parameters of the vertical drains are the same as those in the study of Liu et al. (1998). The calculation parameters are h/(2r w ) = 100, n = 15, s = 1.3, k h /k s = 2, k h /k v = 1.5, k h /k w = 0.00015, and T λ = 0.01κ. κ is changed from 0 to ∞. The time parameter is where T h is the same as that given by Liu et al. (1998). From Figure 4B, it can be observed that the consolidation solution in this study degenerates into an integer-order viscoelasticity model for a completely pervious boundary, which is the same as the U-curve calculated by Liu et al. (1998). Therefore, the solution derived in this study is consistent with the existing theory.

PARAMETER ANALYSIS
From the deduced fractional-derivative viscoelastic vertical drain consolidation solution, it can be easily found that the consolidation characteristics are related to the interface drainage parameters on the top surface of the soil outside vertical  drains c 1 and at the top of vertical drains c 2 . In addition, the viscoelastic parameters of the soil also affect soil consolidation. This study describes the compilation of the program based on the solution mentioned above, mainly calculation and analysis  of the parameters c 1 , c 2 , α, κ, and T λ , and discusses the general law of the consolidation characteristics of fractional derivative viscoelastic vertical drains. In the calculation, the design parameters of the vertical drains are consistent with those given in Huang and Li (2019) and are as follows: h/(2r w ) = 100, n = 15, s = 1.3, k h /k s = 2, k h /k v = 1.5, k h /k w = 0.00015, T λ = 0.1, κ = 10, and α = 0.5. Figure 5 shows the influence of equivalent change in the drainage parameters of the soil top surface interface inside and outside the vertical drains on the U-curves. Figure 6 illustrates the impact of interface drainage parameters on the top surface of the soil outside vertical drains c 1 on the U-curves when interface drainage parameters at the top of vertical drains c 2 remain constant. Figure 7 depicts the impact of the interface drainage parameters c 2 on U-curves when the interface drainage parameters c 1 remain constant. It can be seen from Figures 5-7 that the change in the interface drainage parameters has a crucial impact on U-curves. The pore-water pressure of different interface drainage parameters  dissipates differently with time. The greater the interface drainage parameter, the faster the dissipation of the excess pore water pressure. When the interface drainage parameter is increased to 10 4 , the U-curve approaches a curve. This indicates that when the interface drainage parameter is >10 4 , the drainage boundary approaches the complete drainage boundary. Therefore, the continuous drainage boundary can describe the drainage boundary between the completely drained boundary and the completely undrained boundary, which is more in line with engineering reality. Figure 8 shows the distribution curve of the pore-water pressure v w in the vertical drains and the average pore-water pressurev along the depth range when the interface drainage parameters c 1 and c 2 change equivalently. From the figure, it can be observed that under the same consolidation time, v w andv both decrease with the increase in the interface drainage  parameters of the soil. When the interface drainage parameters are the same, as the consolidation time increases from T h = 1 to T h = 10, the pore-water pressure gradually dissipates. Figure 9 shows the distribution curves of v w andv along the depth range when the interface drainage parameters c 2 remain unchanged and the interface drainage parameters c 1 changes. As shown in Figure 9B, it is clearly seen that the dissipation ofv is affected by the interface drainage parameters c 1 . However, after a certain depth,v is no longer affected by the interface drainage parameters c 1 . Moreover, v w relative tov is less affected by the interface drainage parameters c 1 . Figure 10 shows the distribution curves of v w andv along the depth range when the interface drainage parameters c 1 remain constant and the interface drainage parameters c 2 changes. Figure 10 shows the significant effect of the interface drainage parameters c 2 on pore-water pressures v w andv. Furthermore, in the initial stage, the pore-water pressure v w in the vertical drains is more affected by the interface drainage parameters c 2 than the average pore-water pressurev. Figure 11 illustrates the effect of the fractional order α on Ucurves. In this figure, the other parameters remain unchanged except for the fractional order α. Obviously, Figure 11 shows that in the initial stage, U-curves are very close in value, indicating that pore water pressure is not sensitive to fractional order α. However, in the later stage, the lesser the fractional orderα, the earlier soil consolidation completed, revealing that α has a great impact on the dissipation rate of porewater pressure. This phenomenon is caused by the fractional derivative element. In the later stage of consolidation, when α = 0, the soil is a linear elastic body, and the consolidation speed is fast; when α gradually increases to 1, the soil becomes more viscous; and the effective stress increases slowly, which causes the pore-water pressure to dissipate more slowly. Figure 12 illustrates the effect of the modulus ratio κ on U-curves. In this figure, the other parameters remain unchanged except for κ. Since the modulus E 1 of the independent spring is set to 1, the influence of the independent spring modulus E 1 is eliminated; and the influence of spring 2 is reflected independently. Obviously, the figure shows that the average degree of consolidation decreases when the soil layer modulus ratio increases. Figure 13 illustrates the effect of the creep time factor T λ on U-curves. In this figure, the other parameters remain unchanged except for T λ . The creep time factor is T λ = c v λ/h 2 , which is related to the viscosity coefficient of the soil. From the figure, it is clearly seen that in the early stage, the greater T λ is the faster excess pore-water pressure dissipates. However, in the late stage, the greater T λ , the slower the dissipation rate of the pore-water pressure. This phenomenon is in a good agreement with the 1D consolidation of viscoelastic soil (Xie et al., 2008).

CONCLUSIONS
This study combined continuous drainage boundary conditions and the fractional-derivative model to establish the governing differential equations of vertical drains. After a series of derivation, the consolidated solution was obtained, and the accuracy of the solution was verified by comparison with existing theories. Finally, the solution is used to study the effect of the parameters of interface drainage and viscoelasticity on the behavior of consolidation. Some important conclusions are summarized as follows: (1) Continuous drainage boundary can describe the change in pore-water pressure at the boundary of soil with time.
In addition, the greater the interface drainage parameter, the faster the pore-water pressure dissipates. When the interface drainage parameter is >10 4 , the drainage boundary approaches the complete drainage boundary.
(2) The interface drainage parameters on the top surface of the soil at the boundary outside the vertical drains mainly affect the dissipation rate of average pore-water pressure in the soil layer, when other parameters remain consistent. On the other hand, the interface drainage parameters on the top surface of the vertical drains mainly influence the dissipation rate of the pore-water pressure in the vertical drains, when other parameters remain consistent. (3) Fractional order is a crucial factor affecting the consolidation of viscoelastic vertical drains. In the later stage of consolidation, the greater the fractional order, the more viscous the soil and the slower the progress of soil consolidation. (4) In the initial stage, the greater the creep time factors, the faster the pore-water pressure dissipates; but in the later stage, the situation is just the opposite: the modulus ratio will slow down the dissipation of pore-water pressure during the entire consolidation process.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.