One-Dimensional Consolidation of Multi-Layered Unsaturated Soil with Impeded Drainage Boundaries

: In geotechnical engineering, the consolidation of unsaturated soil is a common issue of great interest. Considering the multi-layered property and impeded drainage boundary condition of the soil stratum in real engineering, this study aimed to develop a general semi-analytical solution for assessing the one-dimensional (1D) consolidation behavior of multi-layered unsaturated soil that is subjected to a general impeded drainage boundary condition and a time-dependent loading. To achieve the ﬁnal solution, the proposed consolidation system is ﬁrstly decoupled and solved in the Laplace domain. Then, the semi-analytical solutions for the excess pore-air pressure and excess pore-water pressures as well as the soil settlement are formulated. The Crump method is employed to provide their ﬁnal results in the time domain. The correctness of the derived solutions was veriﬁed against the available analytical and numerical solutions, and excellent agreements were found for the two comparisons. Moreover, two studied examples are presented to illustrate the 1D consolidation behavior of multi-layered unsaturated soil and the inﬂuences stemming from the impeded drainage parameters are discussed.


Introduction
In geotechnical engineering, the consolidation of unsaturated soil is one of the most important subjects. Because of having both an air and water phase, unsaturated soil shows much more complex properties [1][2][3]. In order to describe the behavior of unsaturated soil undergoing compression and consolidation, Alonso et al. [4] proposed a constitutive model of unsaturated soil on the basis of framework of hardening plasticity with two independent sets of stress variables. Rojas and Chavez [5] developed an elastoplastic framework to account for the volumetric behavior of unsaturated soil. Skutnik et al. [6] discussed the coefficient of consolidation of unsaturated soil on the basis of a set of experimental results. For the consolidation of unsaturated soil, the past few decades have witnessed significant progress in developing various consolidation theories [7][8][9] and the most popular one was propounded by Fredlund and Hasan [10], which used a second-order partial differential system to characterize the dissipating phenomenon of excess pore-air pressure and that of excess pore-water pressure simultaneously. Fredlund and Rahardjo [11] further employed the assumption that all the soil parameters do not change during the consolidation process to develop a simplified form of the previous consolidation system. Afterwards, this consolidation system was expanded to the two-and three-dimensional situations [12,13] and it has also extensively inspired a noticeable number of studies on the consolidation behavior of unsaturated soil [14][15][16]. However, most studies treated both the top surface and bottom base as fully drained or absolutely undrained boundaries for seeking the theoretical solutions more easily. In reality, these boundaries are commonly impeded drainage [17][18][19][20]. To capture this feature, Gray [21] defined a novel impeded drainage boundary by taking advantage of the third kind of boundary condition. This boundary can approach the fully drained one or undrained one by changing the drainage parameter. Using this boundary, Zhou and Zhao [22] proposed a numerical solution for the 1D consolidation system of unsaturated soil on the basis of the differential quadrature method. Wang et al. [23,24] derived semi-analytical solutions for this topic by considering the impeded drainage boundary conditions in one way or the double way. However, these studies are only confined to a simple homogenous soil stratum.
In nature, the soil stratum generally consists of multi-layers due to the complex sedimentation. Under this situation, it seems to be unreasonable to treat this type of soil stratum as a homogenous one. Therefore, it is significant to address the consolidation problem of multi-layered unsaturated soil. To the best knowledge of the authors, however, the studies on the consolidation of multi-layered unsaturated soil are rare. Shan et al. [25] focused on the 1D consolidation system of multi-layered unsaturated soil only catering to the situation that all the soil layers have the same coefficient of air volume change and also have the same coefficient of water volume change. Actually, these coefficients generally vary with increasing depth. Moreover, they only considered ideal single-drainage boundary condition. For the case of the impeded drainage boundary condition, Moradi et al. [26] introduced the well known incremental differential quadrature method (IDQM) to formulate a numerical method for assessing the 1D consolidation behavior of multilayered unsaturated soil. This kind of numerical solution shows a good accuracy and efficiency, but there are still shortcomings in the application to mixed boundary [27,28]. In particular, because of the cumbersome and complex mathematical derivations, it is very difficult to be used by practicing engineers. Therefore, it is crucial to derive an accurate and simple solution for assessing the 1D consolidation behavior of multi-layered unsaturated soil. Compared with the numerical methods, theoretical solutions (analytical and semi-analytical solutions) can be more accurate and simpler. Moreover, they can be also applied to verify the numerical method. To date, however, no effort has been devoted to developing a theoretical solution for assessing the 1D consolidation behavior of the common multi-layer unsaturated soil that caters to a general impeded drainage boundary condition.
This study aimed to propose a general semi-analytical solution for assessing the 1D consolidation behavior of the multi-layered unsaturated soil stratum that is subjected to an impeded drainage boundary condition. The well known governing equations propounded by Fredlund and Hasan [10] are adopted to simulate the consolidation process in each soil layer and formulate the consolidation system of the multi-layered unsaturated soil stratum. To obtain the final solution, the formulated consolidation system is firstly decoupled and solved in the Laplace domain. Then, the semi-analytical solutions for the excess pore-air pressure and excess pore-water pressures as well as the soil settlement are developed. The Crump method is employed to obtain their final results in the time domain. The correctness of the derived solutions is verified with the analytical and numerical solutions from the literature. Moreover, two worked examples are presented to illustrate the 1D consolidation behavior of multi-layered unsaturated soil and the influences stemming from the impeded drainage parameters are discussed.

Governing Equations
Without a loss of generality, as schematically in Figure 1, we consider a 1D consolidation system that consists of several contiguous unsaturated soil layers. Each soil layer is indexed with i (i = 1, 2, . . . , n) and its thickness is h i . Both the top surface and bottom base are treated as impeded drainage to the air and water phases. Assume that a timedependent loading q(t) is uniformly acted on the top of the top soil layer. The flows and the soil settlement are assumed to only occur along the vertical direction. In each soil layer, the coefficient of air permeability and that of water permeability are denoted as k (i) a and k (i) w , respectively. This study adopted all the basic assumptions of 1D consolidation theory for unsaturated soil [1,10], and employed the well known governing equations of Fredlund and Hasan [10] to characterize the 1D consolidation process of each soil layer, that is: where: where the superscript (i) denotes the ith soil layer; u and stand for the derivatives against depth z and time t; and σz is the total normal stress in the vertical direction, respectively.  Figure 1. The one-dimensional consolidation system of multi-layered unsaturated soil under impeded drainage boundaries.
The interactive constants and consolidation parameters in Equation (2) are: The interactive constants and consolidation parameters in Equation (2) are: where m w(i) 1 and m w(i) 2 are the coefficients of the water volume change stemming from a variation in the net normal stress and that stemming from a variation in the matric a0 + u atm and u atm is the atmospheric pressure; g is the acceleration of gravity; R is the universal gas constant; T is the absolute temperature; M is the molecular mass of air; and γ w is the unit weight of water.

Solution Conditions
In each soil layer, the excess pore-air and pore-water pressures are supposed to be uniformly distributed at the initial time, that is: aw are the initial pore-air and pore-water pressures in the ith soil layer, respectively.
Assume that both the top surface and bottom base are impeded drainage to the air and water phases, that is [23,26]: where R ta and R tw are the drainage parameter of the air phase and that of the water phase at the top surface, respectively; and R ba and R bw are the drainage parameters of the air phase and that of the water phase at the bottom base, respectively. The continuity conditions of the excess pore pressures and those of the flow rates at an interface between two adjacent soil layers can be addressed as below: where = 1, 2, · · · , n − 1.
In the matrix form, these boundary and interface continuity conditions can be rewritten as u (1) , where:

Solution Formulations
In the Laplace domain, the excess pore pressure vector can be defined as where s is the Laplace transform parameter. Using the initial conditions (4) to perform the Laplace transform of the system (1) leads to: ,zz − sũ where: Note that the matrix W (i) has the following diagonalized form: with: Reformulate the system (11) using Equation (13) as This system is diagonal for the vector of θ (i) −1ũ (i) and has the general solution: with: where A (i) 1 and A (i) 2 are the two coefficient vectors for the ith soil layer. The Laplace transforms for the boundary conditions (7) and continuity conditions (8) are:ũ Substituting Equation (16) into these equations further produces: (2) . . .
Ψ (2) . . . with Obviously, the coefficient vectors A 2 can be directly solved from Equation (20). Then, we can obtain the excess pore pressures from the inversion of the solution (16) as where j = √ −1; and c is a real constant that exceeds the real part of all the singular points . For a layered unsaturated soil, Equation (22) is commonly difficult to analytically achieve. Herein, the well known Crump's method [29] is employed to approximate its result as where T is a real parameter exceeding one half of the largest time t max .
Once knowing the values of the excess pore-air and pore-water pressures in each soil layer, we can calculate the related volumetric strain ε (i) v from the below expression [9,15]: where m s(i) = m is the coefficient of the soil volume change that is caused by an evolution of the net normal stress; and m is the coefficient of the soil volume change that is caused by an evolution of the matrix suction. This equation takes the following Laplace transform: Integrating this equation against z from z i-1 to z i and summing these integrations of all the soil layers leads to: where S(s) is the Laplace transform of the ground settlement of the multi-layered unsaturated soil; and I is an identity matrix. Similarly, we can adopt the Crump's method [29] to calculate the ground settlement in the time domain as

Verifications
In this part, two consolidation cases of three-layered soil strata are investigated to verify the correctness of the derived semi-analytical solution against the available analytical and numerical solutions. One case is from the work of Shan et al. [25], which has a fully drained top surface and an undrained bottom boundary. For seeking an analytical solution, this case takes the same coefficient of air volume change and that of water volume change for different soil layers. The other case is from the work of Moradi et al. [26], which caters to impeded drainage top and bottom boundaries. In addition, the coefficient of air volume change and that of water volume change are considered to be different in different soil layers. For this case, Moradi et al. [26] proposed an IDQM solution. For these two cases, the specifications of the soil profile and material parameters are given in Table 1 [25,26]. Table 1. Parameters of the two cases of three-layered soil strata adopted from the literature [25,26].
For the first case, the drainage parameters at the top surface and bottom boundary are taken as infinities and zeros (i.e., R ta = R tw = ∞ and R ba = R bw = 0), respectively. This indicates that the top surface is fully drained while the bottom base is undrained to pore-air and pore-water. The applied loading is a constant load of 100 kPa, and at the initial time, both the excess pore-air and pore-water pressures are given as zeros. On the basis of Equation (23), the excess pore pressures are attained and the results are graphically shown in Figure 2a,b in the forms of excess pore pressure distributions. As a comparison, these figures also plot the results of Shan et al. [25] from the analytical solution. For the second case, the top surface is taken as fully drained to pore-air and pore-water, while the bottom base is given as undrained to pore-water but impeded drained to pore-air. The applied loading is an exponential load with the asymptotic value of 100 kPa. The excess pore-air and pore-water pressures are taken as zeros at the initial time. Similarly, the results of excess pore pressures are calculated from Equation (23). For the graphic presentation, the point z = 9.5 m was chosen to be investigated. In Figure 3a,b, the evolutions of excess pore-air and pore-water pressures are depicted. As a comparison, the results of Moradi et al. [26] from the IDQM solution are also given. As seen in Figures 2 and 3, excellent agreements are obtained in these comparisons, which demonstrates the correctness and reliability of the derived semi-analytical solution. Therefore, the proposed solution is feasible to predict the 1D consolidation behavior of multi-layered unsaturated soil stratum with the general impeded drainage boundary conditions. excess pore pressures are calculated from Equation (23). For the graphic presentation, the point z = 9.5 m was chosen to be investigated. In Figure 3a,b, the evolutions of excess poreair and pore-water pressures are depicted. As a comparison, the results of Moradi et al. [26] from the IDQM solution are also given. As seen in Figures 2 and 3, excellent agreements are obtained in these comparisons, which demonstrates the correctness and reliability of the derived semi-analytical solution. Therefore, the proposed solution is feasible to predict the 1D consolidation behavior of multi-layered unsaturated soil stratum with the general impeded drainage boundary conditions. Shan et al. [25] This study

Worked Examples
In this part, two worked examples of three-layered unsaturated soil strata under the impeded drainage boundaries are investigated, and the influences stemming from the drainage parameter of air phase and that of the water phase are discussed. The specifications of the soil profile and material parameters listed in Table 1 for the second case are adopted in the calculations. Two typical types of external loadings including constant and exponential loads are incorporated in the calculation procedure. The excess pore-air and pore-water pressures are given as zeros at the initial time. Herein, at the bottom base, the point z = 9.5 m was chosen to be investigated. In Figure 3a,b, the evolutions of excess poreair and pore-water pressures are depicted. As a comparison, the results of Moradi et al. [26] from the IDQM solution are also given. As seen in Figures 2 and 3, excellent agreements are obtained in these comparisons, which demonstrates the correctness and reliability of the derived semi-analytical solution. Therefore, the proposed solution is feasible to predict the 1D consolidation behavior of multi-layered unsaturated soil stratum with the general impeded drainage boundary conditions.

Worked Examples
In this part, two worked examples of three-layered unsaturated soil strata under the impeded drainage boundaries are investigated, and the influences stemming from the drainage parameter of air phase and that of the water phase are discussed. The specifications of the soil profile and material parameters listed in Table 1 for the second case are adopted in the calculations. Two typical types of external loadings including constant and exponential loads are incorporated in the calculation procedure. The excess pore-air and pore-water pressures are given as zeros at the initial time. Herein, at the bottom base, the

Worked Examples
In this part, two worked examples of three-layered unsaturated soil strata under the impeded drainage boundaries are investigated, and the influences stemming from the drainage parameter of air phase and that of the water phase are discussed. The specifications of the soil profile and material parameters listed in Table 1 for the second case are adopted in the calculations. Two typical types of external loadings including constant and exponential loads are incorporated in the calculation procedure. The excess pore-air and pore-water pressures are given as zeros at the initial time. Herein, at the bottom base, the drainage parameter of the air phase is taken as R ba = 10 and that of the water phase is given as R bw = 5, respectively; while at the top surface, they are changed one by one from 2 to 100. Based on the derived solutions, the excess pore-air and pore-water pressures as well as the ground settlement of these multi-layered unsaturated soil strata are obtained and the results are graphically demonstrated below.
The first worked example takes a changing drainage parameter of air phase at the top surface while that of water phase is taken as R tw = 5. For this situation, Figure 4a,b depict the evolutions of excess pore-air pressure induced by constant and exponential loads, respectively. As can be observed, the action of a constant load immediately results in a great pore-air pressure, which gradually decreases with time passing; while the pore-air pressure first gradually increases as the exponential load increases, and then deceases slowly after the increase in the exponential load becomes slow. In particular, the value of the drainage parameter with respect to the air phase has a pronounced effect on the evolution of pore-air pressure. A larger drainage parameter (i.e., greater drainage efficiency) of air phase delivers a faster dissipation of excess pore-air pressure, especially during the middle and later periods of the consolidation process. Similarly, the evolutions of pore-water pressure are shown in Figure 5a,b. As observed, the evolution of pore-water pressure has two distinct stages. The first one corresponds to the change of pore-air pressure and it is the result of the drainage of pore-air; while the later one is proceeded by the drainage of pore-water. The influence of the drainage parameter above mainly emerges on the first stage, and a larger drainage parameter of air phase leads to a faster dissipation of pore-water pressure. Once the value of this drainage parameter exceeds 100, the curves of pore pressures tend to change as the same way. In addition, for the situation of the constant load, Figure 6a,b demonstrate the isochrones of pore-air pressure and those of pore-water pressure at different times, respectively. As seen, being subjected to a constant load, both the pore-air and pore-water pressures induced in one soil layer are different from those caused in another soil layer. As shown in Table 1, the coefficient of soil volume change in the second layer is the largest one, indicating that the second layer is the softest one. Therefore, being subjected to an external load, this layer may provide a smaller resistance than the other two layers. Under this situation, the air and water phase in the second soil layer have to undertake more load, which means larger pore pressures. It is worthy to note that the drainage parameter (i.e., drainage efficiency) of the air phase shows a profound effect on the distribution of pore-air pressure and that of pore-water pressure. The larger drainage parameter of air phase can speed up their dissipation processes. This effect firstly appears at the position near the top surface, and later spreads far away gradually. The difference resulting from this influence is slightly small in the initial stage, and then it gradually increases with time passing. Later, the induced difference decreases or even disappears. Furthermore, Figure 7a,b show the evolutions of normalized soil settlement under constant and exponential loads, respectively. Being similar to the results in Figure 5a,b, the normalized soil settlement undergoes two distinct developing stages and the influence of the drainage parameter above mainly occurs at the first stage. The larger the drainage parameter of the air phase is, the faster the soil settles. Once the value of this drainage parameter is larger than 100, the curves of the soil settlement also tend to develop as the same pattern. Under this scenario, the influences stemming from the drainage parameter related to the air phase can be ignored, and the top surface approaches full drainage to the air phase.                 The second worked example adopts the changing drainage parameter of the water phase at the top surface whereas that of the air phase is taken as R ta = 5. For this example, the evolutions of excess pore-air pressure under the external excitations of constant load and exponential load are shown in Figure 8a,b, and the evolutions of excess pore-water pressure for the same external excitations are illustrated in Figure 9a,b, respectively. As observed, this drainage parameter has no influence upon the evolution of excess poreair pressure while it shows a pronounced effect on the dissipation of excess pore-water pressure. This influence is different from that caused by the drainage parameter of air phase, as shown in Figures 4 and 5. This phenomenon can be explained that the excess pore-water pressure dissipates significantly after the excess pore-air pressure dissipates completely, as depicted in Figures 5 and 9. Herein, the drainage parameter of the water phase mainly induces a pronounced difference during the later period of the dissipation of excess pore-water pressure. Thus, a change in this drainage parameter will not produce an effect on the evolution of excess pore-air pressure. In general, the larger drainage parameter of the water phase gives rise to a faster dissipation of excess pore-water pressure. When the value of this drainage parameter is larger than 100, its influences on the evolution of excess pore-water pressure can be disregarded. Furthermore, for the situation of the constant load, the distributions of excess pore-air pressures and those of excess pore-water pressures at different times are shown in Figure 10a,b. The results show that the drainage parameter of the water phase does not influence the distributions of excess pore-air pressure at all, but it affects those of excess pore-water pressure significantly. This phenomenon is identical to that in Figures 8 and 9. For the excess pore-water pressure, the larger drainage parameter generally brings about a quicker dissipation. During the consolidation process, the difference caused by this drainage parameter firstly increases to a noticeable value, then gradually decreases or even disappears. Being similar to the phenomena of Figure 6a,b, this drainage parameter initially delivers a noticeable influence at the position near the top surface, and then gradually affects the positions far away. Furthermore, Figure 11a,b demonstrate the changes of the soil settlement induced by the constant and exponential loads, respectively. As seen, the drainage parameter of the water phase mainly produces a noticeable influence during the later period of the soil settling process. Specifically, a larger drainage parameter of the water phase results in a faster settlement rate. When the value of this drainage parameter exceeds 100, its influence in the evolution of the soil settlement can be ignored. Under this scenario, the top surface can be assumed to be a fully drainage boundary for the water phase.
pressure while it shows a pronounced effect on the dissipation of excess pore-water pressure. This influence is different from that caused by the drainage parameter of air phase, as shown in Figures 4 and 5. This phenomenon can be explained that the excess porewater pressure dissipates significantly after the excess pore-air pressure dissipates completely, as depicted in Figures 5 and 9. Herein, the drainage parameter of the water phase mainly induces a pronounced difference during the later period of the dissipation of excess pore-water pressure. Thus, a change in this drainage parameter will not produce an effect on the evolution of excess pore-air pressure. In general, the larger drainage parameter of the water phase gives rise to a faster dissipation of excess pore-water pressure. When the value of this drainage parameter is larger than 100, its influences on the evolution of excess pore-water pressure can be disregarded. Furthermore, for the situation of the constant load, the distributions of excess pore-air pressures and those of excess porewater pressures at different times are shown in Figure 10a,b. The results show that the drainage parameter of the water phase does not influence the distributions of excess poreair pressure at all, but it affects those of excess pore-water pressure significantly. This phenomenon is identical to that in Figures 8 and 9. For the excess pore-water pressure, the larger drainage parameter generally brings about a quicker dissipation. During the consolidation process, the difference caused by this drainage parameter firstly increases to a noticeable value, then gradually decreases or even disappears. Being similar to the phenomena of Figure 6a,b, this drainage parameter initially delivers a noticeable influence at the position near the top surface, and then gradually affects the positions far away. Furthermore, Figure 11a,b demonstrate the changes of the soil settlement induced by the constant and exponential loads, respectively. As seen, the drainage parameter of the water phase mainly produces a noticeable influence during the later period of the soil settling process. Specifically, a larger drainage parameter of the water phase results in a faster settlement rate. When the value of this drainage parameter exceeds 100, its influence in the evolution of the soil settlement can be ignored. Under this scenario, the top surface can be assumed to be a fully drainage boundary for the water phase.

Conclusions
This study incorporates a general impeded drainage boundary condition into 1D consolidation system of multi-layered unsaturated soil. The technique of the Laplace transform and the Crump method are employed to develop a semi-analytical solution for this system. The verifications of the three-layered unsaturated soil conducted in this study exhibit good agreements with the related analytical solution for the traditional boundaries and the numerical incremental differential quadrature method for the impeded drainage boundary condition. The present semi-analytical solution has effectively predicted the evolution of excess pore-air pressure and that of excess pore-water pressures, and it is applicable for the 1D consolidation problems of multi-layered unsaturated soil under various impeded drainage boundaries. The application of an external load to a multi-layered unsaturated soil commonly generates different excess pore-air pressure as well as different excess pore-water pressure in different soil layer. The drainage parameter of the air

Conclusions
This study incorporates a general impeded drainage boundary condition into 1D consolidation system of multi-layered unsaturated soil. The technique of the Laplace transform and the Crump method are employed to develop a semi-analytical solution for this system. The verifications of the three-layered unsaturated soil conducted in this study exhibit good agreements with the related analytical solution for the traditional boundaries and the numerical incremental differential quadrature method for the impeded drainage boundary condition. The present semi-analytical solution has effectively predicted the evolution of excess pore-air pressure and that of excess pore-water pressures, and it is applicable for the 1D consolidation problems of multi-layered unsaturated soil under various impeded drainage boundaries. The application of an external load to a multi-layered unsaturated soil commonly generates different excess pore-air pressure as well as different excess pore-water pressure in different soil layer. The drainage parameter of the air

Conclusions
This study incorporates a general impeded drainage boundary condition into 1D consolidation system of multi-layered unsaturated soil. The technique of the Laplace transform and the Crump method are employed to develop a semi-analytical solution for this system. The verifications of the three-layered unsaturated soil conducted in this study exhibit good agreements with the related analytical solution for the traditional boundaries and the numerical incremental differential quadrature method for the impeded drainage boundary condition. The present semi-analytical solution has effectively predicted the evolution of excess pore-air pressure and that of excess pore-water pressures, and it is applicable for the 1D consolidation problems of multi-layered unsaturated soil under various impeded drainage boundaries. The application of an external load to a multilayered unsaturated soil commonly generates different excess pore-air pressure as well as different excess pore-water pressure in different soil layer. The drainage parameter of the air phase shows great effects on the evolution of excess pore-air pressure and that of excess pore-water pressure as well as the soil settling process; while that of the water phase only affects the evolution of excess pore-water pressure and the soil settling process. These influences initially appear at the position near the boundary and then gradually spread far away. A larger drainage parameter generally expedites the excess pore-pressure dissipations and the soil settling process. Once the drainage parameter is large enough (i.e., larger than 100), the boundary approaches fully drainage.