2 D Plane Strain Consolidation Process of Unsaturated Soil with Vertical Impeded Drainage Boundaries

The consolidation process of soil stratum is a common issue in geotechnical engineering. In this paper, the two-dimensional (2D) plane strain consolidation process of unsaturated soil was studied by incorporating vertical impeded drainage boundaries. The eigenfunction expansion and Laplace transform techniques were adopted to transform the partial differential equations for both the air and water phases into two ordinary equations, which can be easily solved. Then, the semi-analytical solutions for the excess pore-pressures and the soil layer settlement were derived in the Laplace domain. The final results in the time domain could be computed by performing the numerical inversion of Laplace transform. Furthermore, two comparisons were presented to verify the accuracy of the proposed semi-analytical solutions. It was found that the semi-analytical solution agreed well with the finite difference solution and the previous analytical solution from the literature. Finally, the 2D plane strain consolidation process of unsaturated soil under different drainage efficiencies of the vertical boundaries was illustrated, and the influences of the air-water permeability ratio, the anisotropic permeability ratio and the spacing-depth ratio were investigated.


Introduction
The consolidation process of soil stratum is a common issue in geotechnical engineering and it has captured a great deal of attention in the geotechnical community.In the mid-1920s, Terzaghi [1] put forward the classical one-dimensional (1D) consolidation process for saturated soil, after which some researchers derived the solutions to 1D consolidation process for saturated soil subjected to different loading and boundary conditions [2][3][4][5].With the problem raised in practical engineering, much attention has been attracted to the consolidation process of unsaturated soil since the 1960s.Scott [6] discussed the consolidation process of unsaturated soil with occluded air bubbles, and Barden [7,8] presented a study of the compacted unsaturated soil.Furthermore, Fredlund and Hasan [9] proposed a 1D consolidation theory by driving two partial differential equations (PDEs) to describe the dissipation processes of excess pore-air and -water pressures.Dakshanamurthy and Fredlund [10] extended this theory to two-dimensional (2D) plane strain consolidation process by referring the concept of 2D heat diffusion.These theories have been widely accepted in geotechnical engineering and have inspired numerous research projects on the consolidation process of unsaturated soil [11][12][13][14][15].
In order to ensure that analytical solutions were available, most of the previous studies treated the top and bottom drainage boundaries of the soil layer as fully drained or undrained [16].However, these boundaries are generally partially drained (i.e., impeded drainage) in most practical consolidation processes.For example, a sand cushion covered on the top surface of the soil layer, which is commonly used to facilitate the drainage, will become partially drained when its drainage capacity is not effective.The underlying layer will exhibit impeded drainage characteristic when its void ratio is relatively small.In order to describe the partially drained effect, Gray [17] defined an impeded drainage boundary using the third type boundary, which can also realize the fully drained and undrained boundaries (i.e., the first and second type boundaries) by changing the values of the boundary parameters.Base on this boundary, a number of solutions have been proposed to solve the 1D consolidation process of saturated soil [16][17][18][19][20] and that of unsaturated soil [21][22][23].These solutions help to understand the impeding effect (i.e., drainage efficiency) of the covered sand cushion and underlying layer on the 1D consolidation process.However, in the study of the 2D plane strain consolidation processes, particularly for unsaturated soil, there appears no solution incorporating the impeded drainage boundaries.The objective of this paper was to develop a general solution to the 2D plane strain consolidation process of unsaturated soil under such impeded drainage boundaries, that is, the top surface and bottom base of the soil layer were considered to be partially drained.
Based on the consolidation equations proposed by Dakshanamurthy and Fredlund [10], this paper attempted to derive semi-analytical solutions to the 2D plane strain consolidation process of unsaturated soil by incorporating the vertical impeded drainage boundaries.To obtain the final solutions, the eigenfunction expansion and Laplace transform techniques were used to transform PDEs of both air and water phases into two ordinary equations.Then, the solutions for the excess pore-pressures and the soil layer settlement were derived in the Laplace domain, and the final results can be computed by performing the numerical inversion of Laplace transform.Furthermore, two comparisons against the finite difference solution and the previous analytical solution were made to verify the accuracy of the proposed solution.Finally, some computations were conducted to study the 2D plane strain consolidation process of unsaturated soil incorporating vertical impeded drainage boundaries.

Governing Equations
In geotechnical practice, vertical drain assisted with preloading is one of the most commonly used techniques to improve the soft soil.The vertical drain can facilitate the consolidation process by providing a shorter drainage path in the horizontal direction.Figure 1 illustrates a referential consolidation system of a homogeneous unsaturated soil layer with two installed vertical drains, in which the vertical boundaries are presumed to be impeded drainage to the air and water phases.This system extends the problem considered by Ho et al. [14,15].
Processes 2018, 6, x FOR PEER REVIEW 2 of 21 However, these boundaries are generally partially drained (i.e., impeded drainage) in most practical consolidation processes.For example, a sand cushion covered on the top surface of the soil layer, which is commonly used to facilitate the drainage, will become partially drained when its drainage capacity is not effective.The underlying layer will exhibit impeded drainage characteristic when its void ratio is relatively small.In order to describe the partially drained effect, Gray [17] defined an impeded drainage boundary using the third type boundary, which can also realize the fully drained and undrained boundaries (i.e., the first and second type boundaries) by changing the values of the boundary parameters.Base on this boundary, a number of solutions have been proposed to solve the 1D consolidation process of saturated soil [16][17][18][19][20] and that of unsaturated soil [21][22][23].These solutions help to understand the impeding effect (i.e., drainage efficiency) of the covered sand cushion and underlying layer on the 1D consolidation process.However, in the study of the 2D plane strain consolidation processes, particularly for unsaturated soil, there appears no solution incorporating the impeded drainage boundaries.The objective of this paper was to develop a general solution to the 2D plane strain consolidation process of unsaturated soil under such impeded drainage boundaries, that is, the top surface and bottom base of the soil layer were considered to be partially drained.
Based on the consolidation equations proposed by Dakshanamurthy and Fredlund [10], this paper attempted to derive semi-analytical solutions to the 2D plane strain consolidation process of unsaturated soil by incorporating the vertical impeded drainage boundaries.To obtain the final solutions, the eigenfunction expansion and Laplace transform techniques were used to transform PDEs of both air and water phases into two ordinary equations.Then, the solutions for the excess pore-pressures and the soil layer settlement were derived in the Laplace domain, and the final results can be computed by performing the numerical inversion of Laplace transform.Furthermore, two comparisons against the finite difference solution and the previous analytical solution were made to verify the accuracy of the proposed solution.Finally, some computations were conducted to study the 2D plane strain consolidation process of unsaturated soil incorporating vertical impeded drainage boundaries.

Governing Equations
In geotechnical practice, vertical drain assisted with preloading is one of the most commonly used techniques to improve the soft soil.The vertical drain can facilitate the consolidation process by providing a shorter drainage path in the horizontal direction.Figure 1 illustrates a referential consolidation system of a homogeneous unsaturated soil layer with two installed vertical drains, in which the vertical boundaries are presumed to be impeded drainage to the air and water phases.This system extends the problem considered by Ho et al. [14,15].The thickness of the soil layer was H, and the internal spacing of vertical drains was L. In the soil layer, the coefficients of permeability with respect to the air and water in the x-direction were k ax and k wx , and those in z-direction were k az and k wz , respectively.On and under the soil layer were two impedance layers, and their thicknesses were H 0 and H 1 , respectively.These two impedance layers were considered to be permeable in z-direction only.In the upper impedance layer, the coefficients of permeability with respect to the air and water were k a0 and k w0 , respectively; those in the down impedance layer were k a1 and k w1 , respectively.The external load q 0 was assumed to be constant and uniformly distributed on the top of the upper impedance layer.
For a constant external load instantaneously applied, the governing equations for 2D plane strain consolidation process of unsaturated soil can be expressed as follows [10,14,15]: where u a and u w are the excess pore-air and -water pressures, respectively (kPa); C a and C w are the interactive constants associated with the air and water phases, respectively; c a vx and c w vx are the consolidation coefficients with respect to the air and water phases in x-direction, respectively (m 2 /s); and c a vz and c w vz are the consolidation coefficients with respect to the air and water phases in z-direction, respectively (m 2 /s).These consolidation parameters can be expressed as follows: where m a 1 and m w 1 are the coefficients of air and water volume change with respect to the change of the net normal stress, respectively (kPa −1 ); m a 2 and m w 2 are the coefficients of air and water volume change with respect to the change of suction, respectively (kPa −1 ); n is the soil porosity; S r is the degree of saturation; u atm is the atmospheric pressure (kPa); Θ is the absolute temperature (K); R is the universal air constant (8.314J•mol −1 •K −1 ); g is the gravitational acceleration (9.8 m/s 2 ); M is the molecular mass of air phase (0.029 kg/mol); and γ w is the unite weight of water (9.8 kN/m 3 ).

Initial and Boundary Conditions
It was assumed that immediately after the application of the external load, the initial excess pore-air and -water pressures were distributed uniformly throughout the entire stratum.That is: where u 0 a and u 0 w are the values of the initial excess pore-air and -water pressures, respectively.
Processes 2019, 7, 5 4 of 20 The horizontal (i.e., left and right) boundaries, which were defined by two vertical drains, were all considered to permeable to the air and water phases.That is: The vertical (i.e., top and bottom) boundaries, which were separately defined by sand cushion and underlying layer, were considered to be partially permeable to the air and water phases.That is: where R ta = k az0 H/k az H 0 and R tw = k wz0 H/k wz H 0 , and they are two dimensionless characteristic factors describing the drainage efficiency of the top surface to the air and water phases, respectively; and they are two dimensionless characteristic factors describing the drainage efficiency of the bottom base to the air and water phases, respectively.When it is assumed that R ta = R tw = ∞ and R ba = R bw = ∞, the vertical boundaries are double drainage (i.e., top-base drainage [14,15]).For this scenario, both the top and bottom surfaces were permeable to the air and water phases.When it is assumed that R ta = R tw = ∞ and R ba = R bw = 0, the vertical boundaries are single drainage (i.e., top drainage [14,15]).For this scenario, the top surface was permeable to the air and water phase, whereas the bottom base was impermeable to the air and water phases.

Solution Formulation and Verification
In general, unsaturated soil is complex in natural property and engineering behavior due to the lack of homogeneity, the complicate texture, the interaction of multiphase and the nonlinear feature.In order to handle the mathematical derivation easily, some essential assumptions were made along with those used in the 2D plane strain consolidation theory of unsaturated soil [10,14,15].
The prominent assumptions are listed as follows: (1) The soil stratum was considered to be homogeneous; (2) The soil particle and pore water were incompressible; (3) The flows of the air and water phases were continuous and independent, and simultaneously, they only took place in two directions (i.e., the xand z-directions); (4) During the consolidation process, the deformation occurring was small and the coefficients of the volume change for the soil stratum remained constant; (5) The consolidation parameters with respect to the air phase (i.e., C a , c a vx and c a vz ) and those with respect to the water phase (i.e., C w , c w vx and c w vz ) remained constant during the consolidation process; and (6) The effect of air diffusing through water, air dissolving in water, the movement of water vapor and temperature change were ignored.
The above assumptions were not completely accurate for all cases.The consolidation parameters may change because of the variations in the soil properties, such as the permeability coefficients, degree of saturation and porosity, and moreover, the modulus for the soil and the water phase were nonlinear.However, in order to obtain the solution to the consolidation equations of unsaturated soil more easily, it may be acceptable to assume that these soil properties remain constant during the transient process for a small stress increment.Moreover, these constant properties have been commonly adopted to study the consolidation process of unsaturated soil in some existing analytical and semi-analytical methods proposed by Lo et al. [11,14,15], Zhou et al. [12], Shan et al. [13] and Wang et al. [21][22][23].

Solution Formulation
Equations (1a) and (1b) can be rewritten as follows: where κ a = k ax /k az and κ w = k wx /k wz , and they are two anisotropic permeability ratios of the air and water phases, respectively.Through defining two alternative variables φ 1 and φ 2 , Equations (6a) and (6b) can be transformed into: where: The details of derivation process and meanings of c 21 and c 12 can be found in Appendix A. It is noteworthy that Equations (7a) and (7b) will become two uncoupled PDEs when κ a = κ w , and their solutions can be obtained more easily than coupled ones.However, for the sake of generality the solutions are derived below for the coupled case of Equations (7a) and (7b).
Substituting Equation (4) into Equation (8) gives: These boundary conditions above are homogeneous; thus, the general solutions of Equations (7a) and (7b) can be expressed as follows: where ϕ i 1 and ϕ i 2 are the generalized Fourier coefficients of φ 1 and φ 2 with respect to x, respectively; and µ i = iπ/L.
Implementing the Laplace transform on Equation (A22) gives: By substituting Equations (15a) and (15b) into Equation ( 16), the solutions for the excess pore-air and -water pressures can be obtained in the Laplace domain.The results are: where α 1 ~α6 are presented in Equation (A26).
According to the method of two stress-state variables, the volumetric strain of unsaturated soil induced by a constant external load instantaneously applied during the 2D plane strain consolidation process is [14]: where ε v is the volumetric strain; m s 1 = m w 1 + m a 1 , is the coefficient of the soil volume change with respect to the change in net normal stress; m s 2 = m w 2 + m a 2 is the coefficient of the soil volume change with respect to the change in suction.
Applying the Laplace transform to Equation (18) gives: Integrating Equation ( 19) along both the xand z-directions, we can obtain the average settlement of the soil layer in the Laplace domain (denoted as W) as follows: Processes 2019, 7, 5 7 of 20

Solution Verification
Using Crump's method (see Appendix C) [24] to implement the numerical inversion of Laplace transform (NILT) on Equations (17a), (17b) and ( 20), the semi-analytical solutions for the excess pore-air and -water pressures and the soil layer settlement were obtained in the time domain.In this part, two comparisons against the previous analytical solution [14] and the finite difference solution were performed to verify the accuracy of the semi-analytical solution.The details of finite difference solution to Equations (1a) and (1b) with vertical impeded drainage boundaries are presented in Appendix D. The soil properties adopted in these comparisons were assumed as follows [14] approach infinity, the vertical impeded drainage boundaries can be degenerated into those of the top-base drainage system; when R t = ∞ and R t = 0, the vertical impeded drainage boundaries can be simplified into those of the top drainage system.Compared with the results from Ho et al. [14], the variations of the excess pore-air and -water pressures were computed using the semi-analytical solution.Because of limited space, only the isotropic permeability condition for the air and water phases (i.e., κ a = κ w = 1) was considered.For the graphic presentation, the point of investigation was located at x = 1 m and z = 4 m. Figure 2 shows the variations of the excess pore-air and -water pressures under the top-base drainage system, and those under the top drainage system are illustrated in Figure 3.In these two figures, the semi-analytical solution and the previous analytical solution were individually abbreviated to SAS and PAS, respectively.It can be observed that the results of the semi-analytical solution with R b = R t = ∞ were the same as those of previous analytical solution for the top-base drainage system, while the results of the semi-analytical solution with R b = 0 and R t = ∞ were identical to those of the previous analytical solution for the top drainage system.
under the top-base drainage system, and those under the top drainage system are illustrated in Figure 3.In these two figures, the semi-analytical solution and the previous analytical solution were individually abbreviated to SAS and PAS, respectively.It can be observed that the results of the semi-analytical solution with b t R R = =∞ were the same as those of previous analytical solution for the top-base drainage system, while the results of the semi-analytical solution with b 0 R = and t R = ∞ were identical to those of the previous analytical solution for the top drainage system.In order to verify the accuracy of the semi-analytical solution under vertical impeded drainage boundary, another comparison against the finite difference solution was performed.The excess pore-air and -water pressures under the vertical boundaries with R b = R t = 10 were computed using both the semi-analytical solution and the finite difference solution.For the graphic presentation, the point of investigation also was located at x = 1 m and z = 4 m. Figure 4 illustrates the variations of the excess pore-air and -water pressures, where finite difference solution is abbreviated to finite difference solution (FDS).It can be found that the excess pore-air and -water pressures calculated from the semi-analytical solution were in line with those of the finite difference solution.Therefore, these comparisons showed that the proposed semi-analytical solution incorporating the vertical impeded drainage boundaries was correct and it is also feasible for the 2D plane strain consolidation processes of unsaturated soil with top-base drainage and top drainage boundaries.

2D Consolidation Process of Unsaturated Soil with Vertical Impeded Drainage Boundaries
Based on the proposed semi-analytical solution above, some computations were performed to illustrate the 2D plane strain consolidation process of unsaturated soil with the vertical impeded drainage boundaries.The variations of the excess pore-air and -water pressures and the average settlement of the soil layer were presented, and the influences of the drainage efficiency of vertical boundary, the air-water permeability ratio, the anisotropic permeability ratio and the spacing-depth ratio were investigated.Because of limited space, the drainage efficiencies of top surface and bottom base for air phase were only considered to be the same as those for water phases, that is, The parameters of the soil property were assumed to be the same as those used in the previous section.

Consolidation Process with Different Drainage Efficiencies
The vertical drainage boundaries can be classified into symmetrical and asymmetrical ones according to the drainage efficiencies of top and bottom boundaries (i.e., t t b = R R , the vertical drainage boundaries are symmetrical; they are asymmetrical when t b R R ≠ .For these two types of vertical boundaries, the excess-pore pressures and the soil layer settlement were computed using the proposed semi-analytical solutions considering different drainage efficiencies.
For the graphic presentation, the point of investigation on the excess-pore pressures was located at x = 1 m and z = 2.5 m.Figures 5 and 6, respectively, show the variations of the excess pore-pressures and the soil layer settlement under the symmetrical impeded drainage boundary, in which t  Therefore, these comparisons showed that the proposed semi-analytical solution incorporating the vertical impeded drainage boundaries was correct and it is also feasible for the 2D plane strain consolidation processes of unsaturated soil with top-base drainage and top drainage boundaries.

2D Consolidation Process of Unsaturated Soil with Vertical Impeded Drainage Boundaries
Based on the proposed semi-analytical solution above, some computations were performed to illustrate the 2D plane strain consolidation process of unsaturated soil with the vertical impeded drainage boundaries.The variations of the excess pore-air and -water pressures and the average settlement of the soil layer were presented, and the influences of the drainage efficiency of vertical boundary, the air-water permeability ratio, the anisotropic permeability ratio and the spacing-depth ratio were investigated.Because of limited space, the drainage efficiencies of top surface and bottom base for air phase were only considered to be the same as those for water phases, that is, The parameters of the soil property were assumed to be the same as those used in the previous section.

Consolidation Process with Different Drainage Efficiencies
The vertical drainage boundaries can be classified into symmetrical and asymmetrical ones according to the drainage efficiencies of top and bottom boundaries (i.e., R t and R b ).When R t = R b , the vertical drainage boundaries are symmetrical; they are asymmetrical when R t = R b .For these two types of vertical boundaries, the excess-pore pressures and the soil layer settlement were computed using the proposed semi-analytical solutions considering different drainage efficiencies.For the graphic presentation, the point of investigation on the excess-pore pressures was located at x = 1 m and z = 2.5 m.Figures 5 and 6, respectively, show the variations of the excess pore-pressures and the soil layer settlement under the symmetrical impeded drainage boundary, in which R t = R b vary from 0.2 to 100.
these two types of vertical boundaries, the excess-pore pressures and the soil layer settlement were computed using the proposed semi-analytical solutions considering different drainage efficiencies.For the graphic presentation, the point of investigation on the excess-pore pressures was located at x = 1 m and z = 2.5 m.Figures 5 and 6, respectively, show the variations of the excess pore-pressures and the soil layer settlement under the symmetrical impeded drainage boundary, in which t R = b R vary from 0.2 to 100.It can be observed that the drainage efficiency has a significant influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.Larger t R = b R generally resulted in a greater dissipation rate of the excess pore-pressures and a larger settlement rate of the soil layer.Figure 5a graphically indicates that the normalized curves of the excess pore-air pressures gradually shifted to the left with the increase of t R = b R .Figure 5b illustrates the dissipation pattern of water phase and it exhibits the typical double inverse S curves, in which the influence of the drainage efficiency can be separated into two stages: the first one emerged during the dissipation process of the excess pore-air water, and the second one was caused by the dissipation process of the excess pore-air water.Between these two stages, the plateau period emerged on the dissipation curve of the excess pore-air water.Similarly, the influence of the drainage efficiency on the soil layer settlement can also be divided into two stages based on the typical double inverse S curves, as shown in Figure 6.In addition, it can be also found that the differences of the excess pore-pressures and the soil layer settlement induced by different drainage efficiencies decreased with the increase of t R were close to 100, the differences of the excess pore-pressures and the soil layer settlement were negligible.This indicates that the corresponding boundary can be considered to be permeable to air and water phases.Figures 7 and 8, respectively, illustrate the variations of the excess pore-pressures and the soil layer settlement under the asymmetrical vertical impeded drainage boundary, in which b R = 0 or ∞ and t R = 1 or 25.. Being the same as those in Figures 5 and 6, a bigger t R delivered a larger dissipation rate of the excess pore-pressure, and further resulted in a larger settlement rate.The dissipation patterns of the excess pore-air pressure are typical inverse S curves, whereas those of the excess pore-water pressure and the soil layer settlement are typical double inverse S ones with two obvious affected stages.Compared with the results under the symmetrical vertical impeded drainage boundary, Figures 7a and 7b indicate that the excess pore-pressures under the vertical boundaries of b 0 R = dissipated more slowly and consequently resulted in a smaller settlement rate (see Figure 8a).In contrast, Figures 7c and 7d show that the excess pore-pressures under the vertical boundaries of b R = ∞ dissipated faster and consequently induced a larger settlement rate It can be observed that the drainage efficiency has a significant influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.Larger R t = R b generally resulted in a greater dissipation rate of the excess pore-pressures and a larger settlement rate of the soil layer.Figure 5a graphically indicates that the normalized curves of the excess pore-air pressures gradually shifted to the left with the increase of R t = R b .Figure 5b illustrates the dissipation pattern of water phase and it exhibits the typical double inverse S curves, in which the influence of the drainage efficiency can be separated into two stages: the first one emerged during the dissipation process of the excess pore-air water, and the second one was caused by the dissipation process of the excess pore-air water.Between these two stages, the plateau period emerged on the dissipation curve of the excess pore-air water.Similarly, the influence of the drainage efficiency on the soil layer settlement can also be divided into two stages based on the typical double inverse S curves, as shown in Figure 6.In addition, it can be also found that the differences of the excess pore-pressures and the soil layer settlement induced by different drainage efficiencies decreased with the increase of R t = R b .When R t = R b were close to 100, the differences of the excess pore-pressures and the soil layer settlement were negligible.This indicates that the corresponding boundary can be considered to be permeable to air and water phases.
Figures 7 and 8, respectively, illustrate the variations of the excess pore-pressures and the soil layer settlement under the asymmetrical vertical impeded drainage boundary, in which R b = 0 or ∞ and R t = 1 or 25.Being the same as those in Figures 5 and 6, a bigger R t delivered a larger dissipation rate of the excess pore-pressure, and further resulted in a larger settlement rate.The dissipation patterns of the excess pore-air pressure are typical inverse S curves, whereas those of the excess pore-water pressure and the soil layer settlement are typical double inverse S ones with two obvious affected stages.Compared with the results under the symmetrical vertical impeded drainage boundary, Figure 7a,b indicate that the excess pore-pressures under the vertical boundaries of R b = 0 dissipated more slowly and consequently resulted in a smaller settlement rate (see Figure 8a).In contrast, Figure 7c,d show that the excess pore-pressures under the vertical boundaries of R b =∞ dissipated faster and consequently induced a larger settlement rate (see Figure 8b).This phenomenon happened because of the different bottom boundaries.The bottom boundary with R b = 0 was impermeable and it might impede the drainage of the air and water phases; the bottom boundary with R b = ∞ was permeable and could accelerate the drainage process of the air and water phases.

Consolidation Process with Different Air-Water Permeability Ratios
The air-water permeability ratio was defined as the ratio between the vertical permeability coefficient of air phase and that of water phase, i.e., az wz / k k .Moreover, the anisotropic permeability conditions of these phases were assumed to be the same, i.e., a κ = w κ .Figure 9

Consolidation Process with Different Air-Water Permeability Ratios
The air-water permeability ratio was defined as the ratio between the vertical permeability coefficient of air phase and that of water phase, i.e., az wz / k k .Moreover, the anisotropic permeability conditions of these phases were assumed to be the same, i.e., a κ = w κ .Figure 9

Consolidation Process with Different Air-Water Permeability Ratios
The air-water permeability ratio was defined as the ratio between the vertical permeability coefficient of air phase and that of water phase, i.e., k az /k wz .Moreover, the anisotropic permeability conditions of these phases were assumed to be the same, i.e., κ a =κ w .Figure 9 shows the changes of the excess pore-pressures under different values of k az /k wz , in which the vertical impeded drainage boundaries are given by R b = R t = 5 and R b = R t = ∞, respectively.The corresponding variations of the soil layer settlement are shown in Figure 10.It can be observed that the air-water permeability ratio had a great influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.Larger k az /k wz delivered faster dissipation of the excess pore-air pressure at the later stage, and resulted in faster dissipation of the excess pore-water pressure at the intermediate stage.After that, the excess pore-air pressure was fully dissipated, the plateau period emerged on the normalized curve of the excess pore-water pressure, and then, different curves tended to gradually convergence into a single one (see Figure 9b).A similar characteristic can be found in the normalized curves of the soil layer settlement (see Figure 10).Moreover, the bigger k az /k wz can obviously prolong the plateau periods above.In comparison with the results obtained when R b = R t = ∞, the excess pore-pressures and soil layer settlement obtained when R b = R t = 5 were smaller; that is, it needed more time to complete the consolidation process.This phenomenon happened as a result of the impeding effects of the top and bottom drainage boundaries given by R b = R t = 5, and it was identical to the results under 1D consolidation with the impeded drainage boundaries [21].
Processes 2018, 6, x FOR PEER REVIEW 12 of 21 fully dissipated, the plateau period emerged on the normalized curve of the excess pore-water pressure, and then, different curves tended to gradually convergence into a single one (see Figure 9b).A similar characteristic can be found in the normalized curves of the soil layer settlement (see Figure 10).Moreover, the bigger az wz / k k can obviously prolong the plateau periods above.In comparison with the results obtained when b t R R = = ∞, the excess pore-pressures and soil layer settlement obtained when b t R R = 5 = were smaller; that is, it needed more time to complete the consolidation process.This phenomenon happened as a result of the impeding effects of the top and bottom drainage boundaries given by b t R R = 5 = , and it was identical to the results under 1D consolidation with the impeded drainage boundaries [21].

Consolidation Process with Different Anisotropic Permeability Ratios
In this part, for simplicity, the anisotropic permeability ratio of air phase was only considered to be the same as that of water phase, i.e. obvious that the anisotropic permeability ratio had a great effect on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The dissipation with a larger a κ = w κ tended to progress relatively faster than that with a smaller one.The reason is that the vertical drains installed in the soil layer provide a shorter drainage path in the horizontal direction, and they allow the excess pore-pressures dissipate faster when a κ = w κ is larger.Being the same as the results in Figures 5 and 6, the normalized curves of the excess pore-air pressures were a group of fully dissipated, the plateau period emerged on the normalized curve of the excess pore-water pressure, and then, different curves tended to gradually convergence into a single one (see Figure 9b).A similar characteristic can be found in the normalized curves of the soil layer settlement (see Figure 10).Moreover, the bigger az wz / k k can obviously prolong the plateau periods above.In comparison with the results obtained when b t R R = = ∞, the excess pore-pressures and soil layer settlement obtained when b t R R = 5 = were smaller; that is, it needed more time to complete the consolidation process.This phenomenon happened as a result of the impeding effects of the top and bottom drainage boundaries given by b t R R = 5 = , and it was identical to the results under 1D consolidation with the impeded drainage boundaries [21].

Consolidation Process with Different Anisotropic Permeability Ratios
In this part, for simplicity, the anisotropic permeability ratio of air phase was only considered to be the same as that of water phase, i.e. obvious that the anisotropic permeability ratio had a great effect on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The dissipation with a larger a κ = w κ tended to progress relatively faster than that with a smaller one.The reason is that the vertical drains installed in the soil layer provide a shorter drainage path in the horizontal direction, and they allow the excess pore-pressures dissipate faster when a κ = w κ is larger.Being the same as the results in Figures 5 and 6, the normalized curves of the excess pore-air pressures were a group of typical inverse S ones with single influence stage while those of the excess pore-water pressures and

Consolidation Process with Different Anisotropic Permeability Ratios
In this part, for simplicity, the anisotropic permeability ratio of air phase was only considered to be the same as that of water phase, i.e., κ a = κ w .The vertical impeded drainage boundaries are given by R b = R t = 5 and R b = R t = ∞.Figures 11 and 12 illustrate the variations of the excess pore-pressures and the soil layer settlement under different values of κ a = κ w , respectively.It was obvious that the anisotropic permeability ratio had a great effect on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The dissipation with a larger κ a = κ w tended to progress relatively faster than that with a smaller one.The reason is that the vertical drains installed in the soil layer provide a shorter drainage path in the horizontal direction, and they allow the excess pore-pressures dissipate faster when κ a = κ w is larger.Being the same as the results in Figures 5 and 6, the normalized curves of the excess pore-air pressures were a group of typical inverse S ones with single influence stage while those of the excess pore-water pressures and the soil layer settlements were typical double inverse S curves with two obvious influence stages (see Figures 11  and 12).These curves shifted to the left with an increase of κ a = κ w .Being the same as those observed in Figures 9 and 10, under a given value of κ a = κ w , the vertical boundary condition of R b = R t = 5 induced slower dissipations of the excess pore-pressures and slower development of the soil layer settlement than the condition of R b = R t = ∞.For these two boundary conditions, it can be observed that the larger κ a = κ w was, the smaller the difference between these results was.That is, the drainage efficiencies of the vertical impeded boundaries affected the 2D consolidation behavior more slightly when κ a = κ w was larger (e.g., 4).The reason is that for the scenario of larger κ a = κ w , most of the excess pore-pressures were dissipated horizontally, and the vertical drainage effect was small and could even be neglected.conditions, it can be observed that the larger a κ = w κ was, the smaller the difference between these results was.That is, the drainage efficiencies of the vertical impeded boundaries affected the 2D consolidation behavior more slightly when a κ = w κ was larger (e.g., 4).The reason is that for the scenario of larger a κ = w κ , most of the excess pore-pressures were dissipated horizontally, and the vertical drainage effect was small and could even be neglected.

Consolidation Process with Different Spacing-Depth Ratios
For simplicity, only the isotropic permeability condition for air and water phases was considered in this part, and the vertical impeded drainage boundaries were assumed to be symmetrical.Figures 13 and 14, respectively, illustrate the changes in the excess pore-pressures and soil layer settlement under different spacing-depth ratios (i.e., / L H delivered a faster dissipation of the excess-pore pressures due to a shorter horizontal drainage path.As a result, the soil layer settlement develops more quickly.Being the same as the results in Figures 5 and 6, there was a single influence stage on the dissipation pattern of the excess pore-air pressure while there were two obvious influence stages on the dissipation of the excess pore-water pressure and on the development curve of the soil layer settlement (see Figures 13 and conditions, it can be observed that the larger a κ = w κ was, the smaller the difference between these results was.That is, the drainage efficiencies of the vertical impeded boundaries affected the 2D consolidation behavior more slightly when a κ = w κ was larger (e.g., 4).The reason is that for the scenario of larger a κ = w κ , most of the excess pore-pressures were dissipated horizontally, and the vertical drainage effect was small and could even be neglected.

Consolidation Process with Different Spacing-Depth Ratios
For simplicity, only the isotropic permeability condition for air and water phases was considered in this part, and the vertical impeded drainage boundaries were assumed to be symmetrical.Figures 13 and 14, respectively, illustrate the changes in the excess pore-pressures and soil layer settlement under different spacing-depth ratios (i.e., / L H delivered a faster dissipation of the excess-pore pressures due to a shorter horizontal drainage path.As a result, the soil layer settlement develops more quickly.Being the same as the results in Figures 5 and 6, there was a single influence stage on the dissipation pattern of the excess pore-air pressure while there were two obvious influence stages on the dissipation of the excess pore-water pressure and on the development curve of the soil layer settlement (see Figures 13 and

Consolidation Process with Different Spacing-Depth Ratios
For simplicity, only the isotropic permeability condition for air and water phases was considered in this part, and the vertical impeded drainage boundaries were assumed to be symmetrical.Figures 13  and 14, respectively, illustrate the changes in the excess pore-pressures and soil layer settlement under different spacing-depth ratios (i.e., L/H), in which the vertical impeded drainage boundaries are given by R b = R t = 5 and R b = R t = ∞.It can be observed that a smaller L/H delivered a faster dissipation of the excess-pore pressures due to a shorter horizontal drainage path.As a result, the soil layer settlement develops more quickly.Being the same as the results in Figures 5 and 6, there was a single influence stage on the dissipation pattern of the excess pore-air pressure while there were two obvious influence stages on the dissipation of the excess pore-water pressure and on the development curve of the soil layer settlement (see Figures 13 and 14).Under a given L/H, the excess pore-pressures under the vertical boundary condition of R b = R t = 5 dissipated more slowly that those under the vertical boundary condition of R b = R t = ∞.The phenomenon was similar to that observed in Figures 9  and 10.For these two boundary conditions, it was obvious that the difference between their results decreased with the increase of L/H.That is, the drainage efficiencies of vertical impeded boundaries had fewer effects on the 2D consolidation behavior when the horizontal drainage path was shorter (e.g., L/H = 0.5).The reason is that most of the excess pore-pressures were dissipated horizontally for the scenario of shorter horizontal drainage path.R R = = ∞.The phenomenon was similar to that observed in Figures 9 and 10.For these two boundary conditions, it was obvious that the difference between their results decreased with the increase of / L H .That is, the drainage efficiencies of vertical impeded boundaries had fewer effects on the 2D consolidation behavior when the horizontal drainage path was shorter (e.g., / L H = 0.5).
The reason is that most of the excess pore-pressures were dissipated horizontally for the scenario of shorter horizontal drainage path.

Conclusions
In this paper, a semi-analytical solution to 2D plane strain consolidation process of unsaturated soil incorporating vertical impeded drainage boundaries was obtained using eigenfunction expansion and Laplace transform techniques.By performing the inversion of Laplace transform, some computations were performed to study the consolidation processes under different drainage efficiencies of the vertical boundaries, air-water permeability ratios, anisotropic permeability ratios and spacing-depth ratios.As a result of the analysis the following conclusions can be drawn.
a) The present semi-analytical solution was validated to be consistent with the finite difference solution and the previous analytical solution in the literature, and it is more general in the drainage boundary compared with the existing analytical solution.b) The drainage efficiency of the vertical boundary had a great influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The larger drainage efficiency delivered the faster dissipation of the excess pore-pressures and the faster development of the soil layer settlement.When the drainage efficiency was large enough (i.e., t 100 R ≥ or b 100 R ≥ ), the corresponding boundary can be considered to be fully drained.R R = = ∞.The phenomenon was similar to that observed in Figures 9 and 10.For these two boundary conditions, it was obvious that the difference between their results decreased with the increase of / L H .That is, the drainage efficiencies of vertical impeded boundaries had fewer effects on the 2D consolidation behavior when the horizontal drainage path was shorter (e.g., / L H = 0.5).
The reason is that most of the excess pore-pressures were dissipated horizontally for the scenario of shorter horizontal drainage path.

Conclusions
In this paper, a semi-analytical solution to 2D plane strain consolidation process of unsaturated soil incorporating vertical impeded drainage boundaries was obtained using eigenfunction expansion and Laplace transform techniques.By performing the inversion of Laplace transform, some computations were performed to study the consolidation processes under different drainage efficiencies of the vertical boundaries, air-water permeability ratios, anisotropic permeability ratios and spacing-depth ratios.As a result of the analysis the following conclusions can be drawn.
a) The present semi-analytical solution was validated to be consistent with the finite difference solution and the previous analytical solution in the literature, and it is more general in the drainage boundary compared with the existing analytical solution.b) The drainage efficiency of the vertical boundary had a great influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The larger drainage efficiency delivered the faster dissipation of the excess pore-pressures and the faster development of the soil layer settlement.When the drainage efficiency was large enough (i.e., t 100 R ≥ or b 100 R ≥ ), the corresponding boundary can be considered to be fully drained.

Conclusions
In this paper, a semi-analytical solution to 2D plane strain consolidation process of unsaturated soil incorporating vertical impeded drainage boundaries was obtained using eigenfunction expansion and Laplace transform techniques.By performing the inversion of Laplace transform, some computations were performed to study the consolidation processes under different drainage efficiencies of the vertical boundaries, air-water permeability ratios, anisotropic permeability ratios and spacing-depth ratios.As a result of the analysis the following conclusions can be drawn.
(a) The present semi-analytical solution was validated to be consistent with the finite difference solution and the previous analytical solution in the literature, and it is more general in the drainage boundary compared with the existing analytical solution.(b) The drainage efficiency of the vertical boundary had a great influence on the dissipation processes of the excess pore-pressures and the development of the soil layer settlement.The larger drainage efficiency delivered the faster dissipation of the excess pore-pressures and the faster development of the soil layer settlement.When the drainage efficiency was large enough (i.e., R t ≥ 100 or R b ≥ 100), the corresponding boundary can be considered to be fully drained.(c) The consolidation process with a larger air-water permeability ratio tended to progress more quickly than that with a smaller one.Either a larger anisotropic permeability ratio (i.e., larger horizontal permeability coefficient) or a smaller spacing-depth ratio delivered faster dissipation of the excess pore pressures and faster development of the soil layer settlement.In comparison with the results under top-base drainage boundary, the impeded drainage boundary with limited drainage efficiency slowed down the consolidation process significantly. where: where: ( )  The initial and boundary conditions can be obtained as follows: (i) The initial excess pore-air and pore-water pressures ( i = 1 ∼ Nx + 1, j = 1 ∼ Nz + 1):  Based on Equations (A32)-(A43), a computer program can be developed to provide the numerical results of the excess pore-air and pore-water pressures under vertical impeded drainage boundaries.In order to make the computer scheme converge, the time interval and space internal should be reasonably chosen and controlled.For the domain as shown in Figure A1, the mesh sizes should satisfy the following stability conditions [29,30]: x + β a z ≤ 0.5, β w x + β w z ≤ 0.5 (A44)

Figure 2 .Figure 3 .
Figure 2. Excess pore-pressures under top-base drainage boundary calculated from semi-analytical solution (SAS) and previous analytical solution (PAS): (a) air phase and (b) water phase.

Figure 2 .
Figure 2. Excess pore-pressures under top-base drainage boundary calculated from semi-analytical solution (SAS) and previous analytical solution (PAS): (a) air phase and (b) water phase.

Figure 4 .
Figure 4. Excess pore-pressures under impeded drainage boundary calculated from semi-analytical solution (SAS) and finite difference solution (FDS): (a) air phase and (b) water phase.

Figure 4 .
Figure 4. Excess pore-pressures under impeded drainage boundary calculated from semi-analytical solution (SAS) and finite difference solution (FDS): (a) air phase and (b) water phase.

Figure 5 . 21 Figure 5 .
Figure 5. Variations of excess pore-pressures under symmetrical vertical impeded drainage boundary with different drainage efficiencies: (a) air phase and (b) water phase.

1 H=5mFigure 6 .
Figure 6.Variations of soil layer settlement under symmetrical vertical impeded drainage boundary with different drainage efficiencies.

Figure 6 .
Figure 6.Variations of soil layer settlement under symmetrical vertical impeded drainage boundary with different drainage efficiencies.

Figure 7 .
Figure 7. Variations of excess pore-pressures under asymmetrical vertical impeded drainage boundary with different drainage efficiencies: (a,b) b 0 R = and (c,d) b R = ∞ .

Figure 8 .
Figure 8. Variations of soil layer settlement under asymmetrical vertical impeded drainage boundary with different drainage efficiencies: (a) b 0 R = and (b) b R = ∞ .
shows the changes of the excess pore-pressures under different values of az wz / k k , in which the vertical impeded drainage boundaries are given by b , respectively.The corresponding variations of the soil layer settlement are shown in Figure 10.It can be observed that

Figure 8 .
Figure 8. Variations of soil layer settlement under asymmetrical vertical impeded drainage boundary with different drainage efficiencies: (a) b 0 R = and (b) b R = ∞ .
shows the changes of the excess pore-pressures under different values of az wz / k k , in which the vertical impeded drainage boundaries are given by b , respectively.The corresponding variations of the soil layer settlement are shown in Figure10.It can be observed that the air-water permeability ratio had a great influence on the dissipation processes of the excess

Figure 8 .
Figure 8. Variations of soil layer settlement under asymmetrical vertical impeded drainage boundary with different drainage efficiencies: (a) R b = 0 and (b) R b = ∞.

Figure 9 .
Figure 9. Variations of excess pore-pressures under vertical impeded drainage boundary with different air-water permeability ratios: (a) air phase and (b) water phase.

Figure 10 .
Figure 10.Variations of soil layer settlement under vertical impeded drainage boundary with different air-water permeability ratios.
, a κ = w κ .The vertical impeded drainage boundaries are given by b t 5 R R = = and b t R R = = ∞.Figures 11 and 12 illustrate the variations of the excess pore-pressures and the soil layer settlement under different values of a κ = w κ , respectively.It was

Figure 9 .
Figure 9. Variations of excess pore-pressures under vertical impeded drainage boundary with different air-water permeability ratios: (a) air phase and (b) water phase.

Figure 9 .
Figure 9. Variations of excess pore-pressures under vertical impeded drainage boundary with different air-water permeability ratios: (a) air phase and (b) water phase.

Figure 10 .
Figure 10.Variations of soil layer settlement under vertical impeded drainage boundary with different air-water permeability ratios.
, a κ = w κ .The vertical impeded drainage boundaries are given by b t 5 R R = = and b t R R = = ∞.Figures 11 and 12 illustrate the variations of the excess pore-pressures and the soil layer settlement under different values of a κ = w κ , respectively.It was

Figure 10 .
Figure 10.Variations of soil layer settlement under vertical impeded drainage boundary with different air-water permeability ratios.

Processes 2018, 6 , 5 R
x FOR PEER REVIEW 13 of 21 the soil layer settlements were typical double inverse S curves with two obvious influence stages (see Figures11 and 12).These curves shifted to the left with an increase of a κ = w κ .Being the same as those observed in Figures 9 and 10, under a given value of a κ = w κ , the vertical boundary condition of b R t = = induced slower dissipations of the excess pore-pressures and slower development of the soil layer settlement than the condition of b t R R = = ∞.For these two boundary

Figure 11 .
Figure 11.Variations of excess pore-pressures under vertical impeded drainage boundary with different anisotropic permeability ratios: (a) air phase and (b) water phase.

1 Figure 12 .
Figure 12.Variations of soil layer settlement under vertical impeded drainage boundary with different anisotropic permeability ratios.

/
L H ), in which the vertical impeded drainage boundaries are given by b t 5 R R = = and b t R R = = ∞.It can be observed that a smaller

Figure 11 . 5 R
Figure 11.Variations of excess pore-pressures under vertical impeded drainage boundary with different anisotropic permeability ratios: (a) air phase and (b) water phase.

Figure 11 .
Figure 11.Variations of excess pore-pressures under vertical impeded drainage boundary with different anisotropic permeability ratios: (a) air phase and (b) water phase.

1 Figure 12 .
Figure 12.Variations of soil layer settlement under vertical impeded drainage boundary with different anisotropic permeability ratios.

/
L H ), in which the vertical impeded drainage boundaries are given by b t 5 R R = = and b t R R = = ∞.It can be observed that a smaller

Figure
Figure 12.Variations of soil layer settlement under vertical impeded drainage boundary with different anisotropic permeability ratios.

12 .
Figure 12.Variations of soil layer settlement under vertical impeded drainage boundary with different anisotropic permeability ratios.

Processes 2018, 6 , 5 =
x FOR PEER REVIEW 14 of 2114).Under a given / L H , the excess pore-pressures under the vertical boundary condition of dissipated more slowly that those under the vertical boundary condition of b t

Figure 13 .
Figure 13.Variations of excess pore-pressures under vertical impeded drainage boundary with different spacing-depth ratios: (a) air phase and (b) water phase.

Figure 14 .
Figure 14.Variations of soil layer settlement under vertical impeded drainage boundary with different spacing-depth ratios.

Figure 13 . 5 =
Figure 13.Variations of excess pore-pressures under vertical impeded drainage boundary with different spacing-depth ratios: (a) air phase and (b) water phase.

Figure 13 .
Figure 13.Variations of excess pore-pressures under vertical impeded drainage boundary with different spacing-depth ratios: (a) air phase and (b) water phase.

Figure 14 .
Figure 14.Variations of soil layer settlement under vertical impeded drainage boundary with different spacing-depth ratios.

Figure 14 .
Figure 14.Variations of soil layer settlement under vertical impeded drainage boundary with different spacing-depth ratios.

Figure A1 .
Figure A1.Forward time and central space (FTCS) difference scheme for the governing equations: (a) air phase and (b) water phase.

Figure A1 .
Figure A1.Forward time and central space (FTCS) difference scheme for the governing equations: (a) air phase and (b) water phase.