A General Analytical Solution for Axisymmetric Consolidation of Unsaturated Soil with Impeded Drainage Boundaries

The consolidation of soil is one of the most common phenomena in geotechnical engineering. Previous studies for the axisymmetric consolidation of unsaturated soil have usually idealized the boundary conditions as fully drained and absolutely undrained, but the boundaries of unsaturated soil are actually impeded drainage in most practical situations. In this study, we present a general analytical solution for predicting the axisymmetric consolidation behavior of unsaturated soil that incorporates impeded drainage boundary conditions in both the radial and vertical directions simultaneously. The impeded drainage boundary is modeled using the third kind boundary, and it can also realize fully drained and absolutely undrained ones by changing the drainage parameter. A general analytical solution is developed to predict the excess poreair and pore-water pressures as well as the average degree of consolidation in an unsaturated soil stratum using the common methods of eigenfunction expansion and Laplace transform. The newly developed solution is expressed in the product of the terms of time, depth, and radius, which are derived using Laplace transform, usual Fourier, and FourierBessel series, respectively. The eigenfunctions and eigenvalues are evaluated from the impeded drainage boundaries in both radial and depth dimensions. Then, the correctness of the proposed analytical solution is verified against the existing analytical solution for the case of traditional boundaries and against the finite difference solution for the case of general impeded drainage boundaries, and excellent agreements are obtained. Finally, the axisymmetric consolidation behavior of unsaturated soil involving impeded drainage boundaries is demonstrated and analyzed, and the effects of the drainage parameters are discussed. The results indicate that the larger drainage parameter generally expedites the dissipations of the excess pore pressures and further promotes the soil settling process. As the drainage parameter increases, its influence gradually diminishes and even can be neglected when it is larger than 100. The general analytical solution and findings of this study can help for better understanding the axisymmetric consolidation behavior of the unsaturated soil stratum in the ground improvement project with vertical drains as well as the gas-oil gravity drainage mechanism in the naturally fractured reservoirs.


Introduction
The consolidation of soil is one of the most common phenomena in geotechnical engineering [1,2]. In the twenties of the last century, Terzaghi [3] proposed the fundamental theory of one-dimensional (1D) consolidation, which has been of great importance for the progress of geotechnical engineering. To accelerate the consolidation process, vertical drains such as sand drains have been widely used in the improving sites with thick soil deposits. The advantage of this improvement technique is to offer both radial and vertical drainage paths. One of the usual methods to address the vertical drain system is to model it as the axisymmetric consolidation process of a soil cylinder unit around the vertical drain. On the basis of the fundamental theory of Terzaghi [3], a considerable number of studies have been conducted on the axisymmetric consolidation theories for the saturated soil stratum since the late 1940s [4][5][6][7][8][9][10]. In nature, however, most soil beneath the ground is generally unsaturated owing to both the environmental influences and human activities [1,2]. With the great progress made, previous decades have seen lots of well-developed consolidation theories for unsaturated soil, and the typical ones mainly include those of Blight [11], Scott [12], Barden [13,14], and Fredlund and his coworkers [1,2,[15][16][17]. For the axisymmetric consolidation problem of unsaturated soil, Conte [18] proposed a numerical solution using the finite-element analysis method, Qin et al. [19] put forward a semianalytical solution utilizing modified Bessel functions and Laplace transformation, and Zhou [20] introduced the differential quadrature method for proposing a numerical scheme. A common drawback of these solutions is that they consider the air and water phases to only flow along the radial direction. Recently, to capture both the radial and vertical flow behavior of air and water phases, Ho and his coworkers [21,22] transformed the well-known partial differential system of excess pore-air and pore-water pressures from the 3D Cartesian coordinate to the cylindrical one [17]. Using this system, they proposed a set of analytical solutions to assess the axisymmetric consolidation behavior of unsaturated soil under several typical time-dependent loadings [21] and linearly varying initial conditions [22]. These works give valuable knowledge for assessing the axisymmetric consolidation behavior of unsaturated soil, but they only treated both the radial and vertical boundaries as fully drained and/or absolutely undrained ones (i.e., extreme drainage conditions).
In fact, impeded drainage is usually a noticeable matter of concern in most practical consolidation problems. This phenomenon generally occurs when the drainage medium at a boundary takes a poor permeability so that it cannot discharge the air and water phases from the consolidating soil [23,24]. For example, in the ground improving project with a preloading, an overlying sand cushion may be impeded drainage when the fine particle fouls and clogs the drainage path [25] and an underlying layer that has a relatively small void ratio could impede the discharge of air and water phases [26,27]. In particular, for the vertical drain system, the smear effect around a drain and the leakage between two adjoining drains may also affect the consolidation rate of a soil stratum [28,29]. Therefore, it seems to be more general and appropriate to model the boundaries of a consolidating soil as impeded drainage ones. To consider the impeded effect of a boundary, Gray [30] proposed an impeded drainage boundary utilizing the third kind of boundary condition, which can also realize the fully drained and undrained ones through changing the drainage parameter. Using this boundary, researchers have made significant contributions to the consolidation theories of saturated soil, and the effects of the drainage parameters have been investigated in the last decades [24-26, 30, 31]. Recently, for that of unsaturated soil, Wang et al. [27][28][29] introduced the impeded drainage boundary to extend the existing 1D and 2D consolidation theories of unsaturated soil and proposed a series of semianalytical solutions; Zhao et al. [32] developed a new general analytical solution for the 1D consolidation system of unsaturated soil subjected to impeded drainage boundaries. Moradi et al. [33] presented a numerical solution for the 1D consolidation system of layered unsaturated soil with impeded drainage boundaries. Niu et al. [34] developed a set of solutions for 1D consolidation of unsaturated soil with general boundary conditions subjected to a time-dependent load. Huang and Li [35] provided a semianalytical solution for the 2D consolidation system of unsaturated soil that incorporates vertical impeded drainage boundaries. Li et al. [36] proposed a semianalytical solution for the consolidation of drainage well foundations in unsaturated soils with a radial semipermeable drainage boundary. Huang and Zhao [37] derived a general analytical solution for 2D plane strain consolidation of unsaturated soil incorporating lateral and vertical impeded drainage boundaries. So far, for the axisymmetric consolidation system of unsaturated soil, there seems to be no rigorous solution involving the impeded drainage boundaries. The present work is aimed at filling this blankness.
This study is aimed at proposing a reliable axisymmetric consolidation system of unsaturated soil that incorporates the impeded drainage boundaries in both radial and vertical directions. First, the mathematical model is presented. Then, a general analytical solution is derived to calculate the excess pore pressures and the average degree of consolidation using the methods of eigenfunction expansion and Laplace transform. Afterward, the correctness of the newly developed solution is verified against the available analytical solution and finite difference solution. Finally, typical worked examples considering general impeded drainage boundaries are provided to demonstrate the axisymmetric consolidation behavior of unsaturated soil and to assess the effects of the drainage parameters.

Mathematical Model
2.1. Governing Equations. Working in association with preloading, the vertical drains have been used extensively in a ground improvement project for expediting the consolidation process of the soft soil stratum, and they are usually arranged in a regular pattern, such as square pattern and triangular pattern. Figure 1(a) shows the triangular pattern of vertical drains where a typical influence zone is given. Figure 1(b) illustrates the simplified elevation of a typical drain system that penetrates into an unsaturated soil stratum. This system consists of a typical drain with the radius of r w and its influence zone with the radius of r e . The thickness of this system is h. On the top surface of this system, a constant surcharge is applied uniformly, which could produce excess pore-air and pore-water pressures in the soil stratum. Under this scenario, the air and water phases may simultaneously discharge along both the radial and vertical directions. In these two directions, the air permeability takes the coefficients of k ar and k az and the water permeability gets the coefficients of k wr and k wz , respectively. During the mathematical derivations, this study adopts the assumptions that are extensively used to develop the consolidation theory of unsaturated soil. The prominent ones are listed below [15][16][17][18][19][20][21][22][32][33][34][35]: (1) The soil stratum is homogeneous In order to characterize the axisymmetric consolidation phenomenon of unsaturated soil, Ho and Fatahi [21,22] transformed the well-known partial differential system of excess pore-air and pore-water pressures from the 3D Cartesian coordinate to the cylindrical one. For a constant external loading, it is simplified to [21,22] with the vector and matrices where u a and u w are the excess pore-air and pore-water pressures, respectively; c a vr and c w vr denote the radial consolidation coefficient of the air phase and water phase, respectively; c a vz and c w vz denote the vertical consolidation coefficient of the air phase and water phase, respectively; C a and C w denote the interactive constants of these two phases; u ,t is the first derivative of u against time; u ,r and u ,rr are the first and second derivatives of u against radius, respectively; and u ,zz is the second derivative of u against depth. In addition, the aforementioned consolidation parameters are where m a 1 and m w 1 are the coefficients of air and water volume changing with an increase in net stress, respectively; m a 2 and m w 2 are the coefficients of air and water volume changing with an increase in matric suction, respectively; S r denotes the degree of saturation; n denotes the porosity; u atm is the atmospheric pressure; R is the universal gas constant; Θ is the absolute temperature; M is the molecular mass of air; g is the gravitational constant; and γ w is the unit weight of water.

Solution Conditions.
To address the smearing effect induced by vertical drain and the cross-flow between two adjoining drains, Wang et al. [27,28] and Zhou [20] treated the drain surface and that of its influence zone as impeded drainage boundaries, respectively. Likewise, for extending the consolidation theory of unsaturated soil involving vertical drains, Huang and Li [35] introduced vertical impeded drainage boundaries to characterize the clogging effects of the overlying and underlying soil layers.

Geofluids
In the cylindrical coordinate system, these boundary conditions can be expressed as where a t and a b are the drainage parameter at the top surface and at the bottom base, respectively, and a w and a e are the drainage parameter at the drain surface and at the surface of its influence zone, respectively. Obviously, when one of these drainage parameters is zero, the related boundary turns into absolutely undrained; while it is infinite, the related boundary becomes fully drained.
If the external loading uniformly acts on the soil layer in a very large range, it is reasonable and acceptable to assume that the excess pore-air and pore-water pressures induced at the initial time are uniformly distributed in the whole soil layer [38,39]. This gives where u 0 a and u 0 w denote the initial values of u a and u w , respectively.

Solution Derivation
The eigenfunction expansion method is one of the most elegant techniques to solve the linear partial differential system. This method tries to construct the solution of a linear partial differential system as a sum of infinite generalized Fourier series with a complete set of eigenfunctions and the corresponding coefficients. Taking both the radial and vertical boundary conditions (4) and (5) into account, herein we attempt to solve the system (1) by proposing a solution of the following form: where R i ðrÞ and Z j ðzÞ denote the radial eigenfunction and the vertical eigenfunction, respectively; T ij a ðtÞ and T ij w ðtÞ are the generalized Fourier coefficients that relate to the excess pore-air and pore-water pressures, respectively. Using this solution form to separate the system (1) into where T ij = fT ij a ðtÞ, T ij w ðtÞg T ; μ i and ν j are two unknown positive constants. The last three systems are ordinary differential ones that can be more easily solved. Likewise, making the substitution of Equation (7) in the solution conditions (4), (5), and (6), we have where T ij 0 = T ij ð0Þ. Using these solution conditions to solve the ordinary differential systems (8), (9), and (10) one by one, we obtain with parameters where J κ ðμ i rÞ and Y κ ðμ i rÞ are the Bessel function of the first kind and the second kind, respectively; κð= 0, 1Þ is the order of the Bessel function; and μ i and ν j are the radial and vertical eigenvalues, respectively, which are evaluated from the radial and vertical boundary conditions as These two equations are transcendental ones, and they can be solved by some typical numerical methods, such as the dichotomy method and Newton iteration method.
Note that the eigenfunctions R i ðrÞ and Z j ðzÞ satisfy the orthogonality relationships that are presented in Appendix A (i.e., (A.5) and (A.9)). Thus, we can evaluate T ij 0 from Equation (13) as where Then, the combination of Equations (14), (15), and (16) in the solution form (7) gives the final solutions of excess pore-air and pore-water pressures: Once the excess pore-air and pore-water pressures are known, the average volumetric strain arising in the unsaturated soil is computed using the method of two stress-state variables [21,22]: together with where ε v ðtÞ denotes the average volumetric strain arising in the unsaturated soil, m s 1 = m w 1 + m a 1 is the coefficient of soil volume changing with a variation in net stress, m s 2 = m w 2 + m a 2 . is the coefficient of soil volume changing with a variation in matric suction, and the parameter δ ij is Use Equation (22) to calculate the average degree of consolidation of the unsaturated soil as where U avg denotes the average degree of consolidation that represents the soil settling proportion at any given time during the consolidation process.

Solution Verification
To verify the correctness of the developed analytical solution, a finite difference solution (FDS) for the system (1) is formulated in Appendix B using the finite difference method. Based on this solution and the existing analytical solution (EAS) [22], three typical cases are presented and discussed to verify the correctness of the newly developed analytical solution.
The first two take the traditional drainage boundary conditions [22], and the last one gets the impeded drainage boundary conditions. The calculated parameters of the soil layer are taken as [22] r w = 0:2 m, r e = 1:8 m, h = 5 m, k ar = k az = 10 −9 m/s, k wr = k wz = 10 −10 m/s, m s 1 = −5:64 × 10 −4 kPa −1 , m w 1 = 0:2m s 1 , m s 2 = 0:2m s 1 , m w 2 = 1:8m w 1 , n = 0:50,S r = 80%, u 0 a = 20 kPa, u 0 w = 40 kPa, u atm = 100 kPa, and Θ = 293 K. By using these parameters, the excess pore pressures are calculated from Equation (21) and the results are illustrated in Figure 2 for the point located at r = 0:5ðr w + r e Þ and z = 0:5 h. As comparisons, the results given by EAS for the first two cases and those given by FDS for the last case are also plotted in Figure 2. As seen, the results provided by the newly developed analytical solution are in excellent agreements with those from EAS and FDS, powerfully indicating that the newly developed analytical solution is correct. The advantage of the newly developed analytical solution is that it treats both the radial and vertical boundaries as impeded drainage ones, which are more general in the practical engineering, and the existing analytical solution [22] is a special case of this solution.

Worked Example
In this section, we present two worked examples to demonstrate the axisymmetric consolidation behavior of an unsaturated soil stratum incorporating impeded drainage boundaries. For simplicity, we just treat the external surface of the influence zone and the underlying soil layer as absolutely undrained (i.e., a e = 0 and a b = 0) while considering the drain surface and the top surface as impeded drainage (i.e., a w ≠ 0 and a t ≠ 0). In these worked examples, the drainage parameter of the top surface and that of the drain surface are changed one by one from 0.1 to 1000 to show their effects on the consolidation behavior. The parameters of the soil layer are taken as those of the previous section. Using the newly developed analytical solution, the average excess pore-air and pore-water pressures are calculated for these two examples and the results are depicted in Figures 3 and  4. As observed, the excess pore-air pressure dissipates as a 5 Geofluids typical inverse-S pattern while the excess pore-water pressure decreases as a typical double inverse-S pattern. The drainage parameters show noticeable impacts on the dissipation process of excess pore-air pressure and that of excess pore-water pressure. The larger drainage parameter generally expedites these two dissipation processes. Specifically, for the excess pore-air pressure, the drainage parameter mainly influences the later stage of the dissipation process, as depicted in Figure 3, while for the excess pore-water pressure, the drainage parameter makes pronounced differences during two distinct intervals that are separated by an obvious plateau period, as illustrated in Figure 4. The first interval arises concomitantly with the excess pore-air pres-sure dissipating while the second interval emerges as the dissipating plateau period comes to an end. Figure 5 presents the corresponding changes of the average degree of consolidation, which reflects the variations in the soil settling process. As seen, the average degree of consolidation commonly increases more rapidly when the drainage parameters are larger, indicating that the soil layer settles faster. Being similar to that in Figure 4, the drainage parameter also produces an obvious influence on the variation of the average degree of consolidation during two distinct stages. In addition, as the drainage parameter increases, the induced influence gradually diminishes and even can be ignored once the drainage parameter exceeds  100. This indicates that the corresponding boundary comes close to fully drainage.

Geofluids
and that of the drain surface are investigated. The results indicate that the drainage parameter influences the consolidation process significantly. The larger drainage parameter generally expedites the dissipations of the excess pore pressures and further promotes the soil settling process. As the drainage parameter increases, its influence gradually diminishes and even can be neglected when it is larger than 100 (c) In comparison with the existing solution, the developed analytical solution captures the impeded drainage boundary conditions; thus, it is more general for assessing the axisymmetric consolidation behavior of the unsaturated soil stratum in the ground improvement project with vertical drains as well as the gas-oil gravity drainage mechanism in the naturally fractured reservoirs [40,41]. It is noteworthy that further research is needed to determine the impeded drainage parameters of the radial and vertical boundaries in practical applications