Nonlinear Analysis of a Single Vertical Drain Under Vacuum Preloading Based on Axisymmetric Biot’s Consolidation Theory

Abstract

This study incorporates a nonlinear seepage relationship into Biot’s consolidation theory and simulates the consolidation of a single vertical drain under vacuum preloading using the finite element method. The model, simplified via the equal-strain assumption, is validated against theoretical predictions. Under the axisymmetric Biot’s framework, consolidation behavior is analyzed in detail. The results show that in the early stages of consolidation, excess pore water pressure in the vicinity of the prefabricated vertical drain (PVD) does not fully dissipate and may even increase, indicating the occurrence of the Mandel–Cryer effect. As the consolidation process advances, the consolidation front gradually extends outward, and the void ratio near the PVD decreases rapidly, leading to the formation of a clogging zone. In contrast, the reduction in the void ratio in the non-clogging region is relatively slow. The progressive development of the clogging zone significantly impedes the overall consolidation rate. Furthermore, this study explores the influence of key parameters—including the compression index, permeability coefficient, well diameter ratio, smear effect, and well resistance—on the formation of the clogging zone and the Mandel–Cryer effect.

1. Introduction

Vacuum preloading is a widely employed method to accelerate the consolidation of soft soil foundations. Substantial progress has been made in modeling this process through vertical drain consolidation theory, which is fundamentally an axisymmetric problem in soil mechanics [1]. Vacuum preloading is also extensively applied in engineering practice, with its theoretical foundation tracing back to Terzaghi’s classic one-dimensional consolidation theory proposed in 1925. This theory established the first mathematical model for foundation consolidation analysis by quantifying the relationship between pore water pressure dissipation and soil compression, laying the cornerstone for subsequent research. However, the one-dimensional drainage assumption is inadequate for describing the coupled effects of multi-dimensional seepage and lateral deformation inherent in vacuum preloading. To address this limitation, Carrillo introduced a superposition principle for pore pressure under multi-directional seepage in 1942 [2], fostering the development of the early framework for vertical drain consolidation theory. Barron further developed axisymmetric consolidation differential equations in 1948 [3], confirming through comparisons between free-strain and equal-strain assumptions that the equal-strain solution is more suitable for engineering practice. Biot’s consolidation theory further broke through dimensional limitations, achieving the first characterization of total stress’s spatiotemporal evolution via a coupled seepage–displacement field model. With the deepening of the theoretical framework, Horne extended Barron’s equation to radial–axial two-way seepage scenarios in 1964 [4]; Yoshikuni introduced well resistance effects to improve free-strain equations between 1974 and 1979 [5,6]; Hansbo constructed an equal-strain analytical solution incorporating well resistance and smear effects in 1981 [7]; Onoue optimized free-strain theory using radial orthogonal formulations in 1988 [8]; and Xie systematically corrected the deviation in volumetric strain distribution in 1989 [9], forming the current standardized framework for sand drain foundation design.

Subsequent studies have incorporated nonlinear soil behavior into vertical drain consolidation analysis. Basak pioneered the integration of nonlinear soil properties into vertical drain consolidation analysis, focusing on the radial consolidation of ideal sand drains under instantaneous loading scenarios [10]. Lekha considered time-dependent variations in pore water pressure and load increment ratios, deriving closed-form solutions for sand drain consolidation under variable loading conditions [11]. Indraratna, while neglecting well resistance, accounted for smear effects and introduced the e-lgp and e-lgk models to determine soil compressibility and permeability coefficients, and derived approximate solutions for radial nonlinear consolidation using Hansbo’s approach [12]. Zhou argued that these models inadequately describe highly compressible soft soils, thus employing more general expressions for compressibility and permeability to derive nonlinear radial consolidation solutions [13]. Guo investigated nonlinear vertical drain consolidation by considering coupled radial–vertical seepage. He analyzed linear, parabolic, and constant variations in permeability coefficients within the smear zone and provided analytical solutions for foundation consolidation [14]. Han linearized soil compression curves in double logarithmic coordinates lg(1 + e)-lgp and lg(1 + e)-lgk, and obtained approximate solutions for nonlinear consolidation under vacuum loading with depth-dependent attenuation, incorporating well resistance [15]. Gao derived analytical solutions for vertical drain consolidation based on soil nonlinearity and parabolic radial permeability variations in the smear zone, validating the results through degenerate analysis [16]. Li, ignoring well resistance, assumed pore water flow followed Slepicka’s exponential seepage model and obtained analytical solutions for vertical drain consolidation [17]. Lu provided analytical solutions for vertical drain foundations, considering radial–vertical well resistance coupling, variable loading, and double-layered soil conditions [18,19,20].

Schiffman was among the first to explore nonlinear variations in permeability for soft soil and investigated nonlinear consolidation under variable loading [21]. Davis assumed proportional changes in soil compressibility and permeability during consolidation, deriving analytical solutions for one-dimensional nonlinear consolidation under instantaneous loading [22]. Mesri introduced the e-lgk and e-lgp models into nonlinear consolidation analysis and solved them numerically using the finite difference method [23]. Lekha applied these models to equal-strain nonlinear consolidation in vertical drain foundations, successfully deriving analytical solutions for equal-strain nonlinear consolidation [24]. Recent studies have critically examined practical limitations, particularly drainage clogging, a prevalent issue that significantly impairs consolidation efficiency. Comprehensive reviews have synthesized the mechanisms of clogging, including filter clogging, particle migration, and nonuniform consolidation, as well as evaluation metrics and potential mitigation strategies [25]. Experimental investigations have quantified severe permeability reductions (ranging from 50 to 100 times) near drainage boundaries, revealing how nonuniform consolidation occurs rapidly (within 24 to 48 h) and substantially reduces drainage rates by up to 90% [26]. Microstructural studies have further shown that clogging zones, extending from 0.2 to 0.4 times the drainage radius, form due to soil matrix compression under high vacuum gradients. This process induces permeability anisotropy and smear effects through particle reorientation [27].

This study incorporates the Modified Cam-Clay model to investigate nonlinear consolidation of vertical drain foundations under vacuum preloading, considering the variation in permeability coefficient k with void ratio e. The parameters λ and κ of the Modified Cam-Clay model are obtained from one-dimensional consolidation tests. By reducing the proposed model into a poroelastic equal-strain formulation and comparing it with the nonlinear equal-strain analytical solutions considering the e-lgk and e-lgp relationships, the effectiveness of the numerical simulation in this study is verified. Distinct from prior studies, the simulation successfully captures the spatiotemporal evolution of clogging zones and the incipient Mandel–Cryer effect through excess pore pressure transients, quantified by the effectiveness factor. Finally, the influences of key parameters such as the well diameter ratio, compression index, and permeability index on nonlinear consolidation are investigated.

2. Analytical Solution for Nonlinear Radial Consolidation of Vertical Drains Under Equal-Strain Conditions

The analysis model for vertical drains with vacuum preloading is presented in Figure 1. Barron derived the governing equation for radial consolidation under equal-strain conditions, expressed as

(1)${\displaystyle \frac{E_{s}k}{{\gamma }_w}}({\displaystyle \frac{{\partial }^{2}u}{\partial r^2}}+{\displaystyle \frac{\partial u}{\partial r}}{\displaystyle \frac{1}{r}})={\displaystyle \frac{\partial \overline{u}}{\partial t}}$

The governing equation is subject to the following boundary and initial conditions:

For 0 < t and r = r_e (impervious surface), ∂u/∂r = 0 (no seepage on impervious surface);

For 0 < t and r = r_w (permeable surface), u = −80 kPa (vacuum preloading pressure);

At t = 0 and 0 < r ≤ r_e, u = 0 (initial conditions).

Following Xie’s theoretical framework, integrating Equation (1) with boundary condition 1 yields

(2)${\displaystyle \frac{\partial u}{\partial r}}={\displaystyle \frac{{\gamma }_w}{2E_{s}k}}\left({\displaystyle \frac{r_e^2-r^2}{r}}\right){\displaystyle \frac{\partial \overline{u}}{\partial t}}$

Further integration using boundary condition 2 leads to

(3)$\overline{u}+P={\displaystyle \frac{{\gamma }_w}{2E_{s}k}}{\displaystyle \frac{\partial \overline{u}}{\partial t}}(r_e^{2}ln{\displaystyle \frac{r}{r_w}}-{\displaystyle \frac{r^2-r_w^2}{2}})$

where P denotes the absolute vacuum pressure.

From the definition of the average pore pressure at any depth $\overline{u}={\int }_{r_e}^{r_w}2\pi \mathit{rudr}/\pi (r_e^2-r_w^2)$, the relationship between the average pore water pressure and time is given by

(4)$\overline{u}+P={\displaystyle \frac{{\gamma }_w}{E_{s}k}}{F}_n{\displaystyle \frac{\partial \overline{u}}{\partial t}}$

Considering the nonlinear characteristics of consolidation, both E_s and k are functions of the soil’s consolidation state. Specifically, E_s represents the current compression modulus and k, the current permeability coefficient. $F_n={\displaystyle \frac{{r_e}^2}{2}}\left[(lnn-{\displaystyle \frac{3}{4}}){\displaystyle \frac{n^2}{n^2-1}}+{\displaystyle \frac{1}{n^2-1}}(1-{\displaystyle \frac{1}{1-4n^2}})\right]$, $n=r_e/r_w$.

Laboratory consolidation tests reveal that the e-lgp′ relationship for normal consolidated soils follows

(5)$e=e_0-C_c\lg ({\displaystyle \frac{p^{\prime }}{p_0^{\prime }}})$

where e is the current void ratio, e₀ is the initial void ratio, C_c is the compression index, p′₀ is the initial average effective stress, and $p^{\prime }=({\sigma }_x^{\prime }+{\sigma }_y^{\prime }+{\sigma }_z^{\prime })/3$ is the average effective stress.

Similarly, Hvorslev’s linear e-lgk relationship is expressed as [28]

(6)$e=e_0+C_k\lg ({\displaystyle \frac{k}{k_0}})$

where k is the current permeability, k₀ is the initial permeability, and C_k is the permeability index.

The relationship between vertical effective stress and mean effective stress is given by

(7)$d{\sigma }_z^{\prime }={\displaystyle \frac{3}{1+2K_0}}{\mathrm{dp}}^{\prime }$

The vertical strain can be expressed through the void ratio relationship. Combining this with Equation (5) yields

(8)$d{\epsilon }_z=-{\displaystyle \frac{de}{1+e_0}}={\displaystyle \frac{C_c}{(1+e_0)p^{\prime }ln10}}{dp}^{\prime }$

The compression modulus is derived by combining Equations (7) and (8):

(9)$E_s={\displaystyle \frac{d{\sigma }_z^{\prime }}{d{\epsilon }_z}}={\displaystyle \frac{3(1+e_0)p^{\prime }ln10}{C_c(1+2K_0)}}$

Consequently, the relationship between the current compression modulus E_s and the initial compression modulus E_s0 is

(10)${\displaystyle \frac{E_s}{E_{s0}}}={\displaystyle \frac{p^{\prime }}{p_0^{\prime }}}$

Substituting Equation (6) into Equation (5) yields the relationship between the current permeability coefficient k and the initial permeability k₀:

(11)${\displaystyle \frac{k}{k_0}}={\left({\displaystyle \frac{p^{\prime }}{p_0^{\prime }}}\right)}^{{\displaystyle \frac{C_c}{C_k}}}$

By substituting Equations (10) and (11) into Equation (4), the nonlinear relationship between the average pore water pressure $\overline{u}$ and time under equal-strain conditions is obtained:

(12)$\overline{u}+P={\displaystyle \frac{{\gamma }_w}{{E_s}_0{k}_0}}{F}_n{\left({\displaystyle \frac{p_0^{\prime }}{p\prime }}\right)}^{1-{\displaystyle \frac{C_c}{C_k}}}{\displaystyle \frac{\partial \overline{u}}{\partial t}}$

During vacuum preloading, the average effective stress increases from $p_0^{\prime }$ to $p_0^{\prime }+P$. The consolidation degree U based on effective stress is defined as

(13)$U={\displaystyle \frac{p^{\prime }}{p_0^{\prime }+P}}={\displaystyle \frac{p_0^{\prime }-\overline{u}}{p_0^{\prime }+P}}$

Substituting Equation (13) into Equation (12) results in the time-dependent consolidation degree:

(14)${\left(1-U\right)U}^{1-{\displaystyle \frac{C_c}{C_k}}}={\displaystyle \frac{{\gamma }_w}{{E_s}_0{k}_0}}{F}_n{\left({\displaystyle \frac{p_0^{\prime }}{p_0^{\prime }+P}}\right)}^{1-{\displaystyle \frac{C_c}{C_k}}}{\displaystyle \frac{\partial U}{\partial t}}$

3. Nonlinear Finite Element Model for Vertical Drains Based on Axisymmetric Biot’s Consolidation Theory

3.1. Finite Element Model

Based on the Terzaghi–Rendulic consolidation theory, Barron categorized the consolidation behavior of foundations with vertical drains into two idealized deformation patterns: the free strain condition and the equal strain condition. Xie noted that these assumptions represent extreme deformation cases, whereas the actual behavior of soft soils during consolidation typically falls between these two extremes [9].

Therefore, in this study, the finite element software Abaqus (version·2020), utilizing Biot’s consolidation theory, is employed to simulate the axisymmetric nonlinear consolidation of vertical drains. A representative model is established, as illustrated in Figure 2, where r_e = 0.5 m, r_w = 0.05 m, and h = 1 m.

(1). Meshing

According to Xie’s study on finite element meshing for foundations with vertical drains [29], to reduce errors and achieve accurate excess pore water pressure results, mesh refinement is necessary at the cross-section of the soil adjacent to the sand drains and near vertical drainage boundaries. Therefore, the current model adopts the meshing strategy illustrated in Figure 2. The total number of elements and nodes is determined by 45 seeds in the vertical direction and 20 seeds in the horizontal direction, resulting in a grid comprising 900 elements and 945 nodes. In the vertical direction, the influence of vertical drains and gravity is accounted for by refining the mesh at both the top and bottom boundaries. To accurately capture the large stress gradients near the bottom of the drain, a bidirectional biased seeding approach is employed. This results in a finer mesh near the boundaries and a coarser mesh in the middle. In the horizontal direction, significant variations occur near the prefabricated vertical drain (PVD), leading to the division of the domain into four regions. Each region is uniformly subdivided into five grids, with the mesh density being highest near the PVD and gradually coarsening with increasing distance, following a near-dense-to-far-sparse distribution. To improve accuracy and efficiency, a structured mesh is used, employing four-node axisymmetric quadrilateral pore pressure elements (CAX4P).

Regarding mesh density, we conducted a comparative analysis of five different grid schemes, as shown in Table 1.

Comparing the vertical arrangement schemes Cases 1-2-5, it can be found that changes in vertical element density have little effect on the void ratio but have an impact on consolidation time and maximum vertical effective stress. The vertical arrangement of Case 5 is already consistent in terms of maximum vertical effective stress, with only a 1.2% deviation in consolidation time. Comparing the radial arrangement schemes Cases 3-4-5, it is observed that radial arrangement affects all three parameters. In the case of the low-density Case 3, the differences in various values are significant, while Case 4 and Case 5 are very close in terms of consolidation time and void ratio.

Considering both directions comprehensively, Case 5 is the optimal choice in terms of the efficiency of seed arrangement and the validity of calculation.

(2). Analyze Step Setting

Following the Initial step, a Geostatic step is defined to establish an in situ stress equilibrium. Subsequently, a Soils-type consolidation step is performed for a total simulation time of 10⁸ s. The initial time increment for the Soils step is selected according to Equation (15).

(15)$\Delta t\ge {\displaystyle \frac{{\gamma }_w}{6Ek}}(\Delta L{)}^2$

where ΔL is the size of the soil element near the drainage surface, γ_w is the unit weight of water, E is the initial elastic modulus of the soil, and k is the permeability coefficient. For this model, ΔL = 1.24 × 10⁻² m, γ_w = 10⁴ N/m³, and the initial permeability coefficient at the bottom is taken as k = 4.73 × 10⁻¹⁰ m/s. Substituting these values yields Δt ≥ 1838.63 s. Therefore, an initial time increment of 5000 s is adopted in the simulation.

(3). Loads and Boundary Conditions

In the Initial step, the void ratio is set according to stress conditions at the top and bottom boundaries, with initial void ratios of 1.403 and 1.315, respectively. Linear variation with depth is assumed.

Using the Modified Cam-Clay model, convergence fails if the initial mean effective stress at the soil surface is zero. To avoid this, an initial vertical stress of 8 kPa is applied at the top during the Geostatic step to achieve stress equilibrium.

Additionally, a body force of −27 kN/m³ simulates the soil’s self-weight. It is assumed that the phreatic surface coincides with the ground surface, resulting in an initial vertical effective stress increase from −8 kPa at the top to −25 kPa at the bottom. The coefficient of earth pressure at rest is K₀ = μ/(1 − μ) = 0.4286, where μ, Poisson’s ratio, is 0.3.

At the beginning of the Consolidation step, a negative pore pressure of −80 kPa is applied to the left boundary to simulate vacuum preloading.

The summarized boundary conditions of the numerical model are presented in Table 2.

3.2. Soil Parameters

The soil parameters used in this analysis were obtained from soft clay excavated at a construction site in Wenzhou, China. The geological strata in this region are primarily composed of deep silty soft clay, making it representative of typical soft clay deposits in China. The results from the one-dimensional consolidation test conducted on the soil sample are shown in Figure 3, where the e-lgp curve exhibits three characteristic segments: (1) a primary loading path (OCF), (2) unloading path (FGH), and (3) reloading path (HIJ). Notably, Point E identifies the preconsolidation pressure, which is 40 kPa.

In this study, the Modified Cam-Clay model is employed to simulate the nonlinear behavior of the soil [30]. The cap yield surface is implemented with its associated flow rule, employing plastic volumetric strain as the hardening parameter. The parameters for the Modified Cam-Clay model corresponding to this soil are summarized in Table 3.

The permeability coefficient of the soil is calculated using Equation (4), where the permeability index C_k = 0.6.

3.3. Finite Element Model for Nonlinear Consolidation Under Equal-Strain Conditions

To validate the accuracy of the finite element model, the finite element model described above is simplified into a model based on the equal-strain assumption, which is then compared with the analytical solution presented in Section 2. This model reduction is implemented as follows:

First, in Abaqus, the Clay Plastic section within the soil material definition is removed. Then, the value of κ in the Log Bulk Modulus section under Porous Elastic is replaced with the λ value previously defined in the Clay Plastic section. This procedure effectively simplifies the constitutive model from the Modified Cam-Clay model to the nonlinear porous elastic model described in Section 2.

Next, two assumptions are introduced:

(1). the soil undergoes vertical deformation only, with no lateral deformation;

(2). the deformation follows an equal-strain pattern.

To enforce the first assumption, all horizontal displacements in the model are constrained. For the second assumption, a rigid body is placed atop the soil, and a Hard Contact interaction is defined between the rigid body and the soil surface.

With these adjustments, the finite element model is simplified into a Barron-type nonlinear consolidation model under equal-strain conditions, consistent with that described in Section 2.

3.4. Finite Element Model Validation

Based on the equal-strain nonlinear pore water pressure–time relationship derived from Equation (12), a numerical solution is obtained using an iterative method implemented via MATLAB (version·R2022b).

To simulate the range of pore pressure from 0 to −80 kPa induced by vacuum preloading with a step size of du = 0.01 Pa, the pore pressure domain is discretized into N = 8 × 10⁶ steps, and the expression for the pore pressure is u(i) = idu.

The initial stress state of the soil is computed as follows:

(16)${\sigma }_{z0}^{\prime }=8 \mathrm{kPa}$

(17)${\sigma }_z^{\prime }=-du+{\sigma }_{z0}^{\prime }$

(18)${\sigma }_x^{\prime }={\sigma }_y^{\prime }={\sigma }_z^{\prime }\times {\displaystyle \frac{\mu }{1-\mu }}$

Using Equations (16)–(18), the initial average stress p′₀ is calculated by the following:

(19)$p_0^{\prime }={\displaystyle \frac{1}{3}}\left({\sigma }_{z0}^{\prime }+{\displaystyle \frac{2\mu {\sigma }_{z0}^{\prime }}{1-\mu }}\right)$

The average stress at any point is expressed as follows:

(20)$p^{\prime }={\displaystyle \frac{{\sigma }_x^{\prime }+{\sigma }_y^{\prime }+{\sigma }_z^{\prime }}{3}}$

The relationship between the void ratio and the average stress is given by the following:

(21)$e=e_0-\lambda (ln{p}^{\prime }-ln{p}_0^{\prime })$

Using Equation (21), the modulus of elasticity E is related to the porosity ratio and stress state:

(22)$E={\displaystyle \frac{3p^{\prime }(1+e_0\left)\right(1-2\mu )}{\kappa }}$

Then, the constrained modulus of compression E_s is expressed as

(23)$E_s={\displaystyle \frac{E(1-\mu )}{(1+\mu \left)\right(1-2\mu )}}$

The permeability coefficient, derived from Equation (6) can be expressed as

(24)$k=k_0\times 10^{-\left({\displaystyle \frac{e_0-e}{C_k}}\right)}$

Using Equation (12), set $\alpha ={\displaystyle \frac{F_n{\gamma }_w}{k{E}_s}}$, where C_k = 0.6. The relationship between pore pressure and time is derived as dt as follows:

(25)$dt=-{\displaystyle \frac{\alpha }{idu+80000}}{\left({\displaystyle \frac{p_0^{\prime }}{p^{\prime }}}\right)}^{1-{\displaystyle \frac{C_c}{C_k}}}du$

Therefore, the expression of the degree of consolidation characterized by pore pressure is obtained:

(26)$U(i)={\displaystyle \frac{idu}{80000}}$

Based on Equations (25) and (26), the computed consolidation degree for the equal-strain model as a function of time is illustrated in Figure 4. The numerical simulation results from the finite element model are consistent with the theoretical analytical solution, as quantitatively demonstrated by the coherence analysis of the two curves: the mean coherence value is 0.990 with a variance of 0.004, confirming that the finite element method in this study is reasonable and effective.

4. Result Analysis

This section compares the differences and similarities between the nonlinear and linear models. To ensure consistency in the stress state between the two models, an initial consolidation stress of 8 kPa is applied in both models via the Load module. As a result, the effective stress field in the linear model ranges from 8 kPa at the top to 25 kPa at the bottom. Based on the e–ln p relationship, the corresponding void ratios are calculated as 1.403 at the top and 1.315 at the bottom, with both the stress and the void ratio assumed to vary linearly with depth.

In the Property module, the linear model uses the permeability coefficient k, corresponding to the average void ratio at the initial state. This value is assumed to remain constant despite changes in the void ratio. The elastic modulus in the linear model is calculated using Equation (27), where both the initial void ratio and the mean effective stress are determined as the average of the final values at the top and bottom of the nonlinear model:

(27)$E={\displaystyle \frac{3(1-2v\left)\right(1+e_0)p^{\prime }}{\kappa }}$

The parameters used in both models are summarized in Table 4.

4.1. Vertical Stress Spatial Distribution and Settlement Deformation Process

Figure 5 presents a comparison of the spatial distributions of deformation and vertical effective stress at the late stage of consolidation (t = 3.0 × 10⁷ s) for both the nonlinear and linear models. In the linear model, the soil exhibits uniform deformation characteristics during the late stages of consolidation; in other words, vertical deformation at the same elevation is nearly unchanged across the domain. Similarly, vertical effective stress remains almost constant along the horizontal planes. This uniformity arises from the homogeneous distribution of soil properties—such as the coefficient of consolidation—throughout the soil layer. At any given depth, the initial stress state and boundary conditions are the same. Furthermore, due to the path-independent nature of linear elastic problems, the final stress state is unaffected by the consolidation path. As a result, both the stress distribution and the degree of consolidation tend to be consistent at each depth level.

In contrast, the nonlinear model displays significant spatial nonuniformity in deformation. Settlement is markedly greater on the right side of the top surface compared to the left. Additionally, vertical effective stress decreases with increasing distance from the PVD, and a local maximum stress is observed near the bottom adjacent to the drainage boundary. This behavior is primarily attributed to pore pressure dissipation near the PVD, which elevates effective stress and, consequently, the compression modulus. The high-modulus region leads to stress concentration. Moreover, asynchronous consolidation progress across regions causes stress redistribution, further amplifying stress near the drainage boundary. In both models, vertical effective stress increases with depth due to the effect of gravity; however, the nonlinear model demonstrates significantly more complex mechanical behavior.

This complexity fundamentally stems from nonlinear material effects that differentiate constitutive responses. In the linear model, initial stress disequilibrium persists until pore pressure propagates uniformly throughout the domain, after which both stress and strain converge to a symmetric equilibrium configuration. Crucially, the linear model reaches a unique final state independent of loading history. Conversely, the nonlinear model exhibits path-dependent evolution, where both the initial conditions and transient consolidation processes govern the mechanical response.

To further investigate the mechanisms underlying stress development in the two models, a dimensionless position parameter is defined as R = (r − r_w)/(r_e − r_w). At a depth of 0.5 m, two nodes are selected for analysis: one located at R = 1/8 (near the PVD) and one at R = 1 (far from the PVD). The evolution of their stress paths is illustrated in Figure 6. In the nonlinear model, the generalized shear stress q at both locations initially remains constant, followed by a decrease and subsequent increase throughout the consolidation process. However, the onset of the decrease occurs later at the point farther from the PVD. At the final time step, the values of q are 29.7 kPa at R = 1 and 39.0 kPa at R = 1/8. Regarding the evolution of the mean effective stress p’, both positions exhibit a similar trend, but consistently higher values are observed near the PVD, reaching 97.3 kPa at R = 1/8 and 77.7 kPa at R = 1. This nonuniform stress distribution leads to greater changes in the void ratio (i.e., the difference between the initial and final values) near the PVD, resulting in significant differential settlement.

In contrast, the linear model shows a moderate rise–fall trend in generalized shear stress q near the PVD, with a small peak value. The mean effective stress development at R = 1 is similar to that in the nonlinear model; however, the curve at R = 1/8 exhibits a steady increase without significant fluctuations. At the final time step, both positions converge to consistent stress values: 45.7 kPa for q and 57.0 kPa for p′. The uniform stress distribution in the linear model results in consistent changes in the void ratio throughout the soil mass, leading to uniform settlement.

Comparing the stress field responses of the two models reveals that the final generalized shear stress in the linear model is approximately 1.2 to 1.5 times greater than in the nonlinear model. Conversely, the final mean effective stress in the nonlinear model is approximately 1.3 to 1.7 times higher than in the linear model, with the disparity more pronounced near the PVD. This divergence stems from the stress-dependent behavior in the nonlinear model, where both permeability and compressibility vary dynamically. In high-stress regions, decreased permeability induces localized clogging, causing a greater load to be borne by the nearby soil, thereby amplifying void ratio changes. In contrast, the use of constant parameters in the linear model leads to uniform stress transfer, failing to capture the complex behavior exhibited by real soils.

Figure 7 illustrates the relationship between surface settlement and time at the top of the soil for various distances from the PVD under both modeling approaches. The settlement curves reveal several key observations. Throughout the consolidation process, the linear model consistently shows greater settlement during the early stages and reaches its final (stabilized) value earlier than the nonlinear model. Near the PVD, the final settlement values of both models are nearly identical (Figure 7a,b). However, as the distance from the PVD increases, the final settlement of the nonlinear model increases—from 10.3 cm in Figure 7a to 11.5 cm in Figure 7d—whereas the final settlement in the linear model remains nearly constant at approximately 10.3 cm.

At t = 1.2 × 10⁵ s, the settlement at R = 1/8 is 1.04 cm, while at R = 1, it is only 0.25 cm. By t = 8.4 × 10⁶ s, the settlements at both positions converge to nearly identical values of 6.93 cm. At the end of consolidation (t = 5.0 × 10⁷ s), the settlement reaches 10.81 cm at R = 1/8 and 12.33 cm at R = 1. These results suggest that the nonlinear model produces a distinct settlement pattern: it shows greater settlement near the PVD during the early consolidation stage, while the far side experiences more significant settlement in the later stage. This uneven settlement progression ultimately results in slight upward bulging (heaving) near the drainage boundary.

Table 5 consolidates strain data at element centers (depth: 0.05 m) for normalized radial positions R = 1/10 and R = 9/10, captured at two consolidation stages: the initial consolidation phase (t = 10,000 s) and final consolidation completion (t = 6.16 × 10⁷ s). At the initial stage, all strain components at R = 1/10 markedly exceed those at R = 9/10, reflecting intensified radial deformation near the PVD. Upon full consolidation, volumetric strains ε_v converge to comparable magnitudes (−0.1592 at R = 1/10 vs. −0.1555 at R = 9/10). However, greater settlement manifests at R = 9/10—a paradox resolved by contrasting axial strains ε_z. During the late stage of consolidation, this apparent contradiction arises from fundamentally distinct strain states: proximal zones (R = 1/10) experience radial compression (ε_r < 0), whereas distal regions (R = 9/10) develop tensile radial strains (ε_r > 0). Although ε_v = ε_r + ε_z + ε_θ yields similar magnitudes, The dominant ε_z component (−0.1522 at R = 9/10 vs. −0.0957 at R = 1/10) signifies a fundamental shift in the soil deformation mechanism, wherein axial strain governs settlement behavior in distal regions (R = 9/10), while radial drainage dominates proximal zones (R = 1/10). This fundamental strain mechanism governs the observed spatiotemporal settlement evolution.

4.2. Pore Pressure and Mandel–Cryer Effect

Figure 8 illustrates the evolution of pore water pressure over time at different locations in both models. Due to gravitational effects, the initial pore pressure at the bottom of the model is higher than that at the top. As shown in Figure 8a, in the linear model, during the early stage of consolidation, the pore pressure at the bottom near the PVD (R = 1/8) decreases more rapidly than that at the top. Around t = 3.0 × 10⁴ s, the bottom pressure drops to the same level as the top, after which both decrease at an equal rate until stabilizing at −80 kPa. As the distance from the PVD increases (from R = 1/8 to 1), the time required for the top and bottom pore pressures to equalize becomes progressively longer.

In contrast, in the nonlinear model, the bottom pore pressure not only decreases to match that at the top but continues to drop below it before both stabilize at −80 kPa. Moreover, the evolution of pore pressure is more complex in the nonlinear case. Notably, during the early stage of consolidation, the pore pressure at the top initially rises before beginning to dissipate. Meanwhile, at the bottom near the PVD (Figure 8a), the pore pressure initially drops rapidly, then encounters resistance to further dissipation, and ultimately aligns with the top pore pressure. Figure 8d shows that at locations far from the PVD, the initial pore pressure at the top remains nearly constant and close to zero.

Under single-drain consolidation conditions, the degree of consolidation is typically characterized using the time factor $T_v={\displaystyle \frac{kE}{{\gamma }_w{d}_e^2}}t$. However, in the nonlinear model, both the permeability coefficient and the modulus vary with space and time, making this time factor inapplicable for accurately describing the degree of consolidation.

After the application of vacuum loading, the pore water pressure within the soil does not decrease immediately; instead, it increases until it reaches a peak, after which it begins to decline. This phenomenon, first identified by Mandel and Cryer, is referred to as the Mandel–Cryer effect, or the stress transfer effect. To characterize the timing and impact of the Mandel–Cryer effect, a dimensionless pore pressure, denoted as η_t, is defined.

(28)${\eta }_t={\displaystyle \frac{u_t-u_{min}}{u_{max}-u_{min}}}$

In this context, a Mandel–Cryer effect factor M is introduced. This factor is derived by applying a logarithmic transformation to the dimensionless pore pressure η_t. To ensure the validity of the logarithmic function’s domain, an infinitesimal constant $\epsilon ={10}^{-2}$ is added.

(29)$M =-\lg ({\eta }_t+\epsilon )$

Figure 9 illustrates the variation in the Mandel–Cryer effect factor corresponding to pore pressure during the early stage of consolidation in both the nonlinear and linear models. The peak times of the Mandel–Cryer effect at different locations are summarized in Table 6. The timing of the peak value of this effect factor varies with location. In soils near the PVD, the peak of the Mandel–Cryer effect occurs early in the consolidation process (Figure 9a). As the distance from the PVD increases, the peak time of the effect factor is progressively delayed. Comparing the peak times at the top and bottom at the same radial position, it is observed that the peak always occurs earlier at the bottom than at the top. As shown in Figure 9b, the linear model produces the Mandel–Cryer effect only at the top region near the PVD (R < 1/2), with a relatively early peak time. The effect concludes during the early consolidation stage, around t = 2.3 × 10⁵ s.

The engineering implications of the Mandel–Cryer effect warrant rigorous examination. During initial consolidation, this phenomenon induces critical excess pore pressure rises of 2.7 kPa at the top and 3.3 kPa at the bottom. Significantly, these pressures constitute 33.7% of the initial vertical stress at the top (8 kPa) and 13.2% at the bottom (25 kPa), transiently reducing the effective stresses to 5.3 kPa (top) and 21.7 kPa (bottom). Notably, the top zone exhibited exceptionally high susceptibility, with pore pressure rises reaching over one-third of the initial stress state. This significant, albeit temporary, reduction in effective stress can potentially compromise the shear strength and stability of the soil. Critically, the elevated pore pressures, particularly near the PVD, may induce hydraulic fracturing. Furthermore, the risk of localized soil instability in highly sensitive clays experiencing this transient stress state should be considered. Although subsequent stress evolution diminishes the relative impact when pore pressure peaks, the Mandel–Cryer effect remains non-negligible—particularly in top regions during early consolidation—due to its pronounced impact on effective stress and associated risks.

4.3. Spatiotemporal Variation in Void Ratio and Clogging Effect

The clogging effect refers to the substantial reduction in soil permeability near the PVD during vacuum preloading, primarily caused by a decrease in the surrounding soil’s void ratio. This reduction leads to clogging and a sharp decline in drainage capacity. Experimental observations [26] indicate that severe clogging reduces drainage rates by 90% and permeability by 50–100 times, sometimes causing complete blockage. Given the gradual permeability evolution governed by the e-lgk relationship, the compression-induced low-permeability region is referred to as the clogging zone for clarity in subsequent discussions. The onset of clogging is defined as a two-thirds reduction (3-fold decrease) from the initial maximum permeability coefficient (6.64 × 10⁻¹⁰ cm·s⁻¹). While this definition of the clogging zone is inherently bounded by idealized threshold criteria that may not fully capture the continuous, spatially heterogeneous nature of permeability degradation in real soils, it remains a pragmatic and optimized approach for quantifying clogging effects within computational frameworks. Base on the empirical e-lgk relationship defined by Equation (6), where each void ratio e uniquely determines permeability k, the critical void ratio for clogging in the nonlinear model is e_c = 1.03 (corresponding permeability = 2.21 × 10⁻¹⁰ cm·s⁻¹). In contrast, because permeability remains constant in the linear model, the clogging effect does not develop.

Figure 10 depicts the spatial distribution of the void ratio in the nonlinear model at four consolidation stages: early stage (t = 5.0 × 10⁶ s), mid-stage (t = 1.0 × 10⁷ s), late stage (t = 2.0 × 10⁷ s), and completion (t = 6.16 × 10⁷ s). As shown in Figure 10a, the clogging phenomenon (defined as regions where the void ratio e ≤ 1.03, indicated in black) is not prominent during the early stage of consolidation, with the clogging zone appearing only in the upper portion of the PVD. As shown in Figure 10b, the clogging zone begins to extend downward along the PVD and eventually covers its entire length, with a slight increase in its radial extent. During the middle-to-late stages of consolidation, as illustrated in Figure 10c,d, the clogging zone further expands radially, and its thickness continues to increase. By the late stage of consolidation, the maximum thickness of the clogging zone reaches 3/80.

Figure 11 presents the temporal evolution of the void ratio at different locations for both models. As shown in Figure 11a, in the nonlinear model, the void ratio at the top decreases rapidly to match the value at the bottom, after which the rate of decrease slows. This trend indicates that in the early stage of consolidation (t < 3.0 × 10⁵ s), the formation of the clogging zone significantly reduces the consolidation rate and prolongs the overall process. This effect is particularly evident near the PVD, where the clogging zone alters local consolidation behavior.

As shown in Figure 11, in the linear model, the void ratio at the top of the model consistently remains higher than that at the bottom throughout the consolidation process. The void ratio at the top stabilizes earlier, and the stabilization time increases with distance from the PVD. In contrast, in the nonlinear model, the void ratio at the top—particularly near the PVD—decreases much more rapidly than at the bottom, even falling below the bottom value during the early consolidation stage. This results in a denser soil state near the PVD, which aligns with the clogging effect discussed earlier.

A comparison between the nonlinear and linear models reveals distinct differences in the development of top-layer void ratios. Near the PVD (Figure 11a), the linear model shows a noticeably faster rate of reduction than the nonlinear model. Additionally, the time required to reach a stable void ratio is significantly shorter in the linear model. For bottom-layer void ratio development, the linear model generally displays lower values than the nonlinear model within the range R ≤ 1/4 (Figure 11a,b). In the zones where R = 1/2 and R = 1 (Figure 11c,d), both models exhibit similar initial development trends. However, between approximately t = 3.0 × 10⁵ s and t = 8.0 × 10⁵ s, the linear model consistently maintains lower void ratios than the nonlinear model.

As the distance from the PVD increases (from R = 1/8 to 1), the difference in stabilization times for the top and bottom void ratios between the two models gradually diminishes. By the end of consolidation, the void ratios at the bottom in the nonlinear model are 1.06, 1.08, 1.11, and 1.12, indicating that the final void ratio increases with distance from the PVD. In contrast, the final void ratios in the linear model remain nearly uniform (top = 1.17; bottom = 1.07). The fact that the top-layer void ratio in the linear model is higher than in the nonlinear model suggests that the nonlinear model achieves a denser state at the top, indicating better consolidation performance.

In the nonlinear model, the top-layer void ratio curve near the PVD (Figure 11a) exhibits a sharp decline in the early stage of consolidation, followed by gradual stabilization. This indicates that the formation of the clogging zone hinders further consolidation. Additionally, the void ratio at the bottom is significantly higher than at the top upon completion of consolidation, and this difference increases with distance from the PVD. This behavior is primarily attributed to the fully fixed boundary condition at the model base, which restricts soil movement and limits volumetric deformation. As permeability decreases along the drainage path, excess pore pressure caused by external loading cannot dissipate efficiently, further slowing the consolidation rate. Therefore, soil consolidation behavior is influenced not only by soil properties and external loads but also by boundary conditions.

4.4. Degree of Consolidation

The average degree of consolidation of the soil U is defined as follows:

(30)$U={\displaystyle \frac{V_0-V}{V_0-V_1}}$

where V is the total soil volume at a given time during the consolidation process, V₀ is the initial total volume at the beginning of consolidation, and V₁ is the final total volume upon completion of consolidation.

Based on the finite element simulation results showing the evolution of the degree of consolidation over time, as illustrated in Figure 12, a comparison of the two models reveals that the nonlinear model exhibits a slower increase in the degree of consolidation throughout the process. This slower response indicates that the nonlinear e–lgk and e–lgp relationships introduce substantial resistance to consolidation. In contrast, the linear model simplifies soil mechanical behavior but fails to capture complex phenomena that occur in real consolidation, such as the Mandel–Cryer effect and the clogging effect.

The nonlinear model more accurately reflects the intricacies of soil behavior during consolidation, offering a more realistic representation of the time-dependent changes in pore pressure, deformation, and permeability.

5. Analysis of the Clogging Effect, Mandel–Cryer Effect, and Influencing Factors on the Consolidation Rate in the Nonlinear Model

5.1. Influence of Compression Index λ on the Consolidation Process

Figure 13 illustrates the evolution of the Mandel–Cryer effect factor with varying compression indices (λ = 0.1, 0.15, 0.2, and 0.25). The peak times of the Mandel–Cryer effect at different locations are summarized in Table 7. As illustrated, an increase in λ enhanced the compressibility of the soil, leading to a notable extension in the overall consolidation time. However, the peak time of the Mandel–Cryer effect does not increase significantly. Near the PVD (R = 1/8), the peak times are t = 1.0 × 10⁴ s, 1.5 × 10⁴ s, 2.0 × 10⁴ s, and 2.0 × 10⁴ s, respectively. For λ = 0.25, the peak occurs at the first analysis step, indicating a true peak time less than 2.0 × 10⁴ s. This implies that the Mandel–Cryer effect peak time near the PVD initially increases with λ, reaching a maximum at λ = 0.2, then decreases thereafter.

Conversely, at the location farthest from the PVD (R = 1), the peak times for the Mandel–Cryer effect are t = 3.9 × 10⁵ s, 5.8 × 10⁵ s, 5.2 × 10⁵ s, and 7.0 × 10⁵ s, respectively. Although the general trend suggests a gradual increase, overall, the peak time shows an increasing trend; however, at λ = 0.2, the peak time is unexpectedly reduced.

Overall, the peak of the Mandel–Cryer effect progressively shifts outward in the radial direction, away from the PVD. This pattern aligns with the observed settlement behavior: areas near the PVD experience faster initial settlement and higher consolidation rates, while more distant regions undergo delayed settlement and reduced rates of consolidation.

Figure 14 illustrates the evolution of the clogging zone in soils with different compression indices (λ = 0.1, 0.15, 0.2, and 0.25). As shown in Figure 14a, when the compression index is relatively low (λ = 0.1), no clogging zone is observed. However, as illustrated in Figure 14b–d, increasing the compression index (λ = 0.15 to 0.25) significantly shortens the times for the onset of clogging. The times needed for the clogging zone to develop downward and cover half the height of the PVD are t = 8.0 × 10⁶ s, 5.2 × 10⁵ s, and 8.0 × 10⁴ s, respectively. The corresponding times required for the clogging zone to develop downward and cover the bottom of the PVD are t = 1.5 × 10⁷ s, 2.6 × 10⁶ s, and 1.3 × 10⁶ s. Additionally, the spatial extent of the affected area expands significantly with an increasing λ.

Notably, as the compression index increases (Figure 14b,c), the clogging zone rapidly expands. At a higher compression index (Figure 14d, black curve), the clogging zone nearly spans the entire soil domain. This phenomenon is primarily attributed to the increased compressibility of the soil at higher compression indices. Under the same stress level, a higher λ results in a lower void ratio, which further reduces the permeability coefficient. Consequently, the formation of the clogging zone occurs much earlier, and its spatial influence becomes substantially more extensive.

Engineering practice must account for the positive correlation between soil compressibility and clogging severity, whereby higher compression indices accelerate clogging zone development while expanding its spatial influence (Figure 14). Consequently, soils exhibiting elevated compressibility require proactive anti-clogging measures.

Figure 15 presents consolidation curves for varying compression indices (λ = 0.1, 0.15, 0.2, and 0.25). As illustrated, higher λ values significantly prolong the consolidation process, with completion times increasing from 2.8 × 10⁷ s (λ = 0.1) to 1.6 × 10⁸ s (λ = 0.25). This progression demonstrates that soil compressibility exerts a substantial influence on consolidation duration, where a 150% increase in λ corresponds to a 450% extension in consolidation time.

5.2. Influence of Permeability Coefficient C_k on the Consolidation Process

Figure 16 illustrates the variation in the Mandel–Cryer effect factor in soils with different permeability coefficients (C_k = 0.4, 0.6, 0.8, and 1.0). The peak times of the Mandel–Cryer effect at different locations are summarized in Table 8. At a location far from the PVD (R = 1), the peak times of the Mandel–Cryer effect occur at t = 5.2 × 10⁵ s, 5.8 × 10⁵ s, 5.8 × 10⁵ s, and 3.9 × 10⁵ s, respectively. The overall trend exhibits a non-monotonic behavior, where the peak time initially increases and then decreases, reaching a maximum at C_k = 0.6–0.8.

Conversely, at a location near the PVD (R = 1/8), the peak times are t = 2.0 × 10⁴ s, 1.5 × 10⁴ s, 1.0 × 10⁴ s, and 1.0 × 10⁴ s, respectively, showing a consistently decreasing trend as C_k increases. This indicates that near the PVD, a higher permeability coefficient accelerates the dissipation of excess pore pressure, causing the Mandel–Cryer peak to occur earlier.

The non-monotonic behavior of the peak time farther from the PVD may be attributed to slower drainage at low permeability, which delays the Mandel–Cryer peak. As permeability increases, drainage efficiency improves and the peak time shortens. However, when permeability becomes excessively high, the pore pressure gradients tend to flatten, resulting in a smoother spatial distribution and, consequently, a reduced Mandel–Cryer peak time.

Figure 17 illustrates the development of the clogging zone in soils with varying permeability coefficients (C_k = 0.4, 0.6, 0.8, and 1.0). As the permeability coefficient increases, the initiation of clogging zone formation occurs progressively earlier. The times required for the clogging zone to develop and cover half the height of the PVD are t = 1.0 × 10⁷ s, 6.0 × 10⁶ s, 4.5 × 10⁶ s, and 3.9 × 10⁶ s, respectively. Likewise, the times to fully cover the PVD are t = 2.0 × 10⁷ s, 1.0 × 10⁷ s, 7.8 × 10⁶ s, and 6.9 × 10⁶ s. However, the rate of time reduction diminishes with increasing permeability.

This trend can be attributed to the fact that a higher permeability coefficient facilitates more efficient drainage, thereby accelerating the formation of the clogging zone. Simultaneously, the enhanced dissipation of excess pore pressure leads to a gradual increase in effective stress within the soil, further promoting earlier development of the clogging zone.

Nevertheless, although permeability significantly influences the timing of clogging zone formation, it does not substantially affect its final spatial configuration, as evidenced by the black curves in Figure 16. In other words, the permeability coefficient primarily affects the temporal characteristics of the consolidation process, while exerting minimal influence on the ultimate spatial extent of the clogging zone.

Figure 18 depicts consolidation curves under varying permeability coefficients (C_k = 0.4, 0.6, 0.8, and 1.0). The consolidation process accelerates substantially with increased permeability, evidenced by completion times decreasing from 1.0 × 10⁸ s (C_k = 0.4) to 4.0 × 10⁷ s (C_k = 1.0). This inverse relationship demonstrates that permeability governs consolidation efficiency, where a 150% increase in C_k reduces consolidation time by approximately 60%. However, the relative reduction in consolidation time diminishes progressively with higher C_k values, indicating a nonlinear dependency between permeability enhancement and time-saving.

5.3. Influence of the Well Diameter Ratio r_e/r_w on the Consolidation Process

Figure 19 illustrates the evolution of the Mandel–Cryer effect in soils with varying well diameter ratios (r_e/r_w = 5, 10, 15, 20). The peak times of the Mandel–Cryer effect at different locations are summarized in Table 9. By comparing points with the same ratio of distance (relative position) from the PVD across different models, the peak times at the locations far from the PVD (R = 1) are observed to be t = 8.1 × 10⁴ s, t = 5.4 × 10⁵ s, t = 1.4 × 10⁶ s, and t = 2.6 × 10⁶ s, respectively. At the locations near the PVD (R = 1/8), the peak times are t = 5.0 × 10³ s, t = 1.0 × 10⁴ s, t = 3.0 × 10⁴ s, and t = 6.0 × 10⁴ s, respectively. These results indicate that the peak time increases with the well diameter ratio for the same relative position as the distance from the PVD increases.

By comparing spatial nodes at identical absolute positions from the PVD across different models (Figure 19a, R = 1 and r_e/r_w = 5; Figure 19b, R = 1/2 and r_e/r_w = 10; and Figure 19d, R = 1/4 and r_e/r_w = 20), the peak times of the effect are observed to be t = 8.1 × 10⁴ s, t = 1.9 × 10⁵ s, and t = 2.2 × 10⁵ s, respectively. Similarly, from Figure 19a (R = 1/2 and r_e/r_w = 5), Figure 19b (R = 1/4 and r_e/r_w = 10), and Figure 19d (R = 1/8 and r_e/r_w = 20), the corresponding peak times are t = 2.8 × 10⁴ s, t = 3.8 × 10⁴ s, and t = 6.0 × 10⁴ s. These results suggest that the propagation speed of the Mandel–Cryer effect decreases as the well diameter ratio increases, even when the absolute positions remain unchanged. This phenomenon can be attributed to the fact that while the absolute positions are consistent, a larger well diameter ratio leads to a slower consolidation process, thereby causing a delay in the peak time of the Mandel–Cryer effect.

Figure 20 illustrates the evolution of the clogging effect in soils with varying well diameter ratios (r_e/r_w = 5, 10, 15, and 20). As the well diameter ratio increases, the onset of clogging zone formation is progressively delayed. The times required for the clogging zone to extend across half of the PVD are t = 2.1 × 10⁶ s, 3.0 × 10⁶ s, 6.5 × 10⁶ s, and 1.2 × 10⁷ s, respectively. Similarly, the times required for the zone to fully extend across the entire PVD are t = 5.8 × 10⁶ s, 8.4 × 10⁶ s, 2.0 × 10⁷ s, and 3.4 × 10⁷ s, respectively. Overall, the formation time exhibits a clear linear positive correlation with the well diameter ratio.

A larger well diameter ratio leads to more uniform drainage conditions and smoother pore pressure dissipation, which in turn delay the initiation of the clogging zone. Consequently, the appearance of the clogging zone is linearly postponed as the well diameter ratio increases. Moreover, the final consolidation time also increases accordingly, as a larger well diameter ratio results in an extended drainage path, thereby slowing the overall consolidation rate.

However, despite its significant influence on the timing of clogging zone formation and the duration of consolidation, the well diameter ratio does not substantially affect the final shape or spatial extent of the clogging zone.

Figure 21 presents consolidation curves under varying well diameter ratios (r_e/r_w = 5, 10, 15, and 20). Consolidation duration exhibits strong geometric dependence, with completion times increasing from 1.3 × 10⁷ s (r_e/r_w = 5) to 2.3 × 10⁸ s (r_e/r_w = 20). This demonstrates that drain spacing governs consolidation efficiency, where a 300% increase in r_e/r_w extends consolidation time by 1690%, consistent with theoretical inverse-square relationships in radial drainage systems.

5.4. Influence of Smear Effect and Well Resistance on the Consolidation Process

For the model considering the smearing effect, the radius of the smeared zone is 150 mm, and the height of the smeared zone is 1000 mm. The ratio of the permeability coefficient of the undisturbed soil to that of the smeared zone (k_h/k_s) is assumed to be 2.0. As for the model considering the well resistance effect, the attenuation value of the well resistance is 50%, which means that the vacuum preloading suction linearly attenuates from −80 kPa at the top to −40 kPa.

Figure 22 illustrates the evolution of the Mandel–Cryer effect in soils influenced by the smear effect and well resistance. The peak times of the Mandel–Cryer effect at different locations are summarized in Table 10. On the side far from the PVD (R = 1), the time points when the Mandel–Cryer effect reaches its peak are t = 5.8 × 10⁵ s, t = 7.8 × 10⁵ s, and t = 5.5 × 10⁵ s, respectively. On the side close to the PVD (R = 1/8), the time points when the Mandel–Cryer effect reaches its peak are t = 1.5 × 10⁴ s, t = 3.0 × 10⁴ s, and t = 1.5 × 10⁴ s, respectively. It can be observed that the original model and the model considering well resistance exhibit consistent responses to the Mandel–Cryer effect. For the model considering the smear effect, on the sides close to the PVD (R = 1/8, 1/4), the times for the Mandel–Cryer effect to reach its peak are approximately twice those of the original model. However, as the distance from the PVD increases, this gap narrows, eventually reaching approximately 1.35 times.

Figure 23 illustrates the evolution of the clogging effect in soils influenced by the smear effect and well resistance. It can be observed that in the model considering the smear effect, the clogging zone develops very rapidly. By t = 2.0 × 10⁴ s, it has already fully covered the PVD, and its final extent surpasses the smear zone. This indicates that the presence of the smear effect is highly sensitive to the response of the clogging zone. In contrast, the model accounting for well resistance shows a much slower and more limited development of the clogging zone. At t = 1.1 × 10⁷ s, it covers only approximately 1/4 of the PVD, and the final morphology of the clogging zone does not even cover half of the PVD. Its consolidation completion time is relatively close to that of the original model, which suggests that the presence of well resistance is not conducive to the formation of the clogging zone but exerts a slight impeding effect on consolidation progress.

Figure 24 presents consolidation curves for the model with the smear effect, model with well resistance, and original model. When comparing the model considering the smear effect, the model considering well resistance, and the original model, their consolidation completion times are 9.2 × 10⁷ s, 6.3 × 10⁷ s, and 6.1 × 10⁷ s, respectively. The presence of well resistance prolongs the consolidation time by approximately 3%, while the presence of the smear effect extends it by about 51%. These results indicate that although well resistance slightly reduces the efficiency of vacuum preloading, its associated clogging zone remains relatively limited and does not significantly affect the overall consolidation time. In contrast, the smear effect notably accelerates the development of the clogging zone, resulting in a substantially greater impact on consolidation behavior.

6. Limitations of This Study

This study has corresponding limitations, as follows. Firstly, the numerical model established in this study was not validated by any physical on-site measurements or laboratory data, and its results represent theoretical analyses. Consequently, any conclusions regarding practical applications must be treated with considerable caution. Secondly, the Modified Cam-Clay model employed is primarily applicable to normally consolidated clays, with limited suitability for other soil types such as structured or cemented clays. This restricts the generalizability of this study’s findings to soils outside this classification.

7. Conclusions

This study investigates the consolidation behavior of soil under vacuum preloading using an axisymmetric nonlinear consolidation framework. The following conclusions are drawn:

A nonlinear finite element model for a single vertical drain under vacuum preloading is established based on the axisymmetric Biot’s consolidation theory. By degenerating the general consolidation equation, an equal-strain nonlinear consolidation equation is derived. A comparison between the consolidation curves generated by this model and the analytical solution of the equal-strain nonlinear formulation demonstrates strong agreement, thereby validating the accuracy and reliability of the proposed approach.

The nonlinear model effectively captures the settlement pattern during consolidation, which features rapid settlement near the PVD in the early stage and accelerated settlement at more distant locations in the later stage, leading to the formation of a bulging zone adjacent to the drain. The model successfully represents the development of a densified zone with reduced permeability and comprehensively accounts for both gravity and nonlinear seepage behavior.

The finite element model successfully simulates the development of a densified zone with reduced permeability using void ratio data, providing empirical support for the decline in consolidation efficiency observed in practice. The results reveal the evolution of the Mandel–Cryer effect in the early consolidation stage: radially, the closer the location is to the PVD, the earlier the peak pore pressure occurs; vertically, the peak appears earlier at positions nearer to the bottom. As consolidation progresses outward, rapid void ratio reduction near the PVD leads to clogging, which significantly impedes overall consolidation. In contrast, the Mandel–Cryer effect in the linear model is weak, appearing only briefly near the top of the PVD at the start, and it fails to produce a densified zone due to its restrictive constitutive assumptions.

A parameter sensitivity analysis is conducted to evaluate the influence of the compression index λ, permeability coefficient C_k, well radius ratio r_e/r_w, smear effect, and well resistance on soil consolidation behavior. The results indicate that the Mandel–Cryer peak time is relatively insensitive to λ but increases initially and then decreases with a rising C_k. As r_e/r_w increases, the peak time at both relative and absolute positions within the soil increases. Additionally, the smear effect prolongs the Mandel–Cryer peak time, while well resistance has essentially no impact on it. The development of a densified zone with reduced permeability is highly sensitive to λ, with larger values significantly enhancing its development. However, variations in C_k and r_e/r_w influence only the timing, not the final spatial configuration, of the densified zone. Furthermore, the smear effect significantly expands this zone, whereas well resistance delays its emergence.

Future research should focus on the following areas:

(a). Experimental validation to further refine model accuracy, particularly in relation to clogging development and spatial settlement patterns.

(b). Exploring the relationship between MCC yielding/hardening and consolidation-induced stress anomalies, with an emphasis on stress redistribution mechanisms under axisymmetric drainage conditions.

Author Contributions

Conceptualization, X.P. and D.L.; methodology, X.P. and D.L.; software, D.L.; formal analysis, D.L. and J.F.; writing—original draft preparation, D.L. and X.P.; writing—reviewing and editing, X.P. and X.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author(s).

Acknowledgments

This project was supported by the Engineering Research Center of the Ministry of Education for Renewable Energy Infrastructure Construction Technology.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Axisymmetric consolidation model of a single vertical drain foundation: (a) analysis scheme; (b) unit body.

Figure 2 Finite element mesh configuration: (a) overall dimensions; (b) lateral arrangement.

Figure 3 Void ratio versus vertical stress from the one-dimensional consolidation test.

Figure 4 Comparison of the degree of consolidation between the analytical solution and the numerical solution for the nonlinear equal-strain model.

Figure 5 Vertical stress distribution and deformation profile during consolidation (t = 3.0 × 10⁷ s): (a) nonlinear model; (b) linear model.

Figure 6 Stress development curves: (a) nonlinear model; (b) linear model.

Figure 7 Surface settlement curves at different locations: (a) R = 1/8; (b) R = 1/4; (c) R = 1/2; (d) R = 1.

Figure 8 Pore pressure variation at different locations: (a) R = 1/8; (b) R = 1/4; (c) R = 1/2; (d) R = 1.

Figure 9 Variation in the Mandel–Cryer effect factor at the early stage of consolidation: (a) nonlinear model; (b) linear model.

Figure 10 Development of the clogging zone at different consolidation stages: (a) early stage (t = 5.0 × 10⁶ s); (b) middle stage (t = 1.0 × 10⁷ s); (c) late stage (t = 2.0 × 10⁷ s); (d) final stage (t = 6.16 × 10⁷ s).

Figure 11 Void ratio evolution at different positions: (a) R = 1/8; (b) R = 1/4; (c) R = 1/2; (d) R = 1.

Figure 12 Consolidation curves for the nonlinear and linear models.

Figure 13 Mandel–Cryer effect factor for different compression indices: (a) λ = 0.1; (b) λ = 0.15; (c) λ = 0.2; (d) λ = 0.25.

Figure 14 Development of the clogging zone for different compression indices: (a) λ = 0.1; (b) λ = 0.15; (c) λ = 0.2; (d) λ = 0.25.

Figure 15 Consolidation curves for different compression indices λ.

Figure 16 Mandel–Cryer effect for different permeability coefficients: (a) C_k = 0.4; (b) C_k = 0.6; (c) C_k = 0.8; (d) C_k = 1.0.

Figure 17 Development of the clogging zone for different permeability coefficients: (a) C_k = 0.4; (b) C_k = 0.6; (c) C_k = 0.8; (d) C_k = 1.0.

Figure 18 Consolidation curves for different permeability coefficients C_k.

Figure 19 Mandel–Cryer effect under different well diameter ratios: (a) r_e/r_w = 5; (b) r_e/r_w = 10; (c) r_e/r_w = 15; (d) r_e/r_w = 20.

Figure 20 Development of the clogging effect for different well diameter ratios: (a) r_e/r_w = 5; (b) r_e/r_w = 10; (c) r_e/r_w = 15; (d) r_e/r_w = 20.

Figure 21 Consolidation curves for different well radius ratios r_e/r_w.

Figure 22 Mandel–Cryer effect considering smear effect and well resistance: (a) model with smear effect; (b) model with well resistance.

Figure 23 Development of clogging effect in (a) model with smear effect; (b) model with well resistance; (c) original model.

Figure 24 Consolidation curves for model with smear effect, model with well resistance, and original model.

Table 1

Meshing schemes.

Meshing Schemes	Vertical	Horizontal	Consolidation Time(s)	σ′max-z (kPa)	emin
Case 1	10 uniform elements	As in Figure 2b	6.09 × 107	206.8	0.982
Case 2	100 uniform elements	As in Figure 2b	6.23 × 107	201.5	0.981
Case 3	As in Figure 2a	8 uniform elements	6.02 × 107	174.9	1.017
Case 4	As in Figure 2a	80 uniform elements	6.15 × 107	208.0	0.972
Case 5	As in Figure 2a	As in Figure 2b	6.16 × 107	201.5	0.981

Table 2

Boundary conditions applied in the numerical model.

Location	Type	Magnitude
Top Boundary	Pressure	−8 kPa
Bottom Boundary	Displacement/Rotation	Fully Fixed
Left Boundary	Displacement/Rotation	Fixed in X-direction
Left Boundary	Pore pressure	−80 kPa
Right Boundary	Displacement/Rotation	Fixed in X-direction

Table 3

Parameters of Modified Cam-Clay model.

Parameters	Value
γw (kN/m3)	10
Poisson’s ratio, μ	0.3
Log bulk modulus, κ	0.0261
Log bulk plastic modulus, λ	0.1565
Stress ratio, M	1
Initial yield surface size, a0 (kPa)	0
Wet yield surface size, β	1
Flow stress rate, k	1

Table 4

Nonlinear and linear model parameters.

Parameters	Nonlinear	Linear
Log Bulk Modulus, κ	0.0261	/
Log Bulk Plastic Modulus, λ	0.1565	/
Elastic Modulus, E (MPa)	/	0.58
Poisson’s Ratio, μ	0.3	0.3
Initial Yield Surface Size, a0 (kPa)	0	/
Top/Bottom Initial Void Ratio, e0	1.403/1.315	1.403/1.315
$\mathrm{Top}/\mathrm{Bottom} \mathrm{Initial} \mathrm{Vertical} \mathrm{Effective} \mathrm{Stress}, {\sigma }_v^{\prime }$ (kPa)	8/25	8/25
Top/Bottom Initial Permeability Coefficient, k	6.64 × 10−10/4.73 × 10−10	5.60 × 10−10/5.60 × 10−10

Table 5

Strain component evolution at radial positions during consolidation.

Time	Location	Type	Value
1.0 × 104 s	R = 0.1/0.9	εr	0.013/0.0004
	R = 0.1/0.9	εz	−0.014/−0.0003
	R = 0.1/0.9	εθ	−0.012/−6.5 × 10−5
	R = 0.1/0.9	εv	−0.013/3.5 × 10−5
6.16 × 107 s	R = 0.1/0.9	εr	−0.0413/0.0075
	R = 0.1/0.9	εz	−0.0957/−0.1522
	R = 0.1/0.9	εθ	−0.0222/−0.0108
	R = 0.1/0.9	εv	−0.1592/−0.1555

Table 6

Summary of peak times of the Mandel–Cryer effect in nonlinear and linear models.

Nonlinear		Linear
Location	Peak Time	Location	Peak Time
Top, R = 1/8	1.5 × 104 s	Top, R = 1/8	6.8 × 103 s
Top, R = 1/4	3.9 × 104 s	Top, R = 1/4	2.3 × 104 s
Top, R = 1/2	1.8 × 105 s	Top, R = 1/2	/
Top, R = 1	5.8 × 105 s	Top, R = 1	/
Bottom, R = 1/8	5.0 × 103 s	Bottom, R = 1/8	/
Bottom, R = 1/4	5.0 × 103 s	Bottom, R = 1/4	/
Bottom, R = 1/2	2.8 × 104 s	Bottom, R = 1/2	/
Bottom, R = 1	1.2 × 105 s	Bottom, R = 1	/

Table 7

Summary of peak times of the Mandel–Cryer effect under different compression indices λ.

Model	Location	Peak Time	Model	Location	Peak Time
λ = 0.1	R = 1/8	1.0 × 104 s	λ = 0.2	R = 1/8	2.0 × 104 s
λ = 0.1	R = 1/4	2.8 × 104 s	λ = 0.2	R = 1/4	5.5 × 104 s
λ = 0.1	R = 1/2	1.2 × 105 s	λ = 0.2	R = 1/2	1.6 × 105 s
λ = 0.1	R = 1	3.9 × 105 s	λ = 0.2	R = 1	5.2 × 105 s
λ = 0.15	R = 1/8	1.5 × 104 s	λ = 0.25	R = 1/8	2.0 × 104 s
λ = 0.15	R = 1/4	3.9 × 104 s	λ = 0.25	R = 1/4	6.0 × 104 s
λ = 0.15	R = 1/2	1.8 × 105 s	λ = 0.25	R = 1/2	2.2 × 105 s
λ = 0.15	R = 1	5.8 × 105 s	λ = 0.25	R = 1	7.0 × 105 s

Table 8

Summary of peak times of the Mandel–Cryer effect under different permeability coefficients C_k.

Model	Location	Peak Time	Model	Location	Peak Time
Ck = 0.4	R = 1/8	2.0 × 104 s	Ck = 0.8	R = 1/8	1.0 × 104 s
Ck = 0.4	R = 1/4	5.5 × 104 s	Ck = 0.8	R = 1/4	3.9 × 104 s
Ck = 0.4	R = 1/2	1.6 × 105 s	Ck = 0.8	R = 1/2	1.8 × 105 s
Ck = 0.4	R = 1	5.2 × 105 s	Ck = 0.8	R = 1	5.8 × 105 s
Ck = 0.6	R = 1/8	1.5 × 104 s	Ck = 1.0	R = 1/8	1.0 × 104 s
Ck = 0.6	R = 1/4	3.9 × 104 s	Ck = 1.0	R = 1/4	3.9 × 104 s
Ck = 0.6	R = 1/2	1.8 × 105 s	Ck = 1.0	R = 1/2	1.8 × 105 s
Ck = 0.6	R = 1	5.8 × 105 s	Ck = 1.0	R = 1	3.9 × 105 s

Table 9

Summary of peak times of the Mandel–Cryer effect under different well radius ratios r_e/r_w.

Model	Location	Peak Time	Model	Location	Peak Time
re/rw = 5	R = 1/8	5.0 × 103 s	re/rw = 15	R = 1/8	3.0 × 104 s
re/rw = 5	R = 1/4	1.0 × 104 s	re/rw = 15	R = 1/4	1.1 × 105 s
re/rw = 5	R = 1/2	2.8 × 104 s	re/rw = 15	R = 1/2	4.8 × 105 s
re/rw = 5	R = 1	8.1 × 104 s	re/rw = 15	R = 1	1.4 × 106 s
re/rw = 10	R = 1/8	1.0 × 104 s	re/rw = 20	R = 1/8	6.0 × 104 s
re/rw = 10	R = 1/4	3.9 × 104 s	re/rw = 20	R = 1/4	2.2 × 105 s
re/rw = 10	R = 1/2	1.9 × 105 s	re/rw = 20	R = 1/2	8.6 × 105 s
re/rw = 10	R = 1	5.4 × 105 s	re/rw = 20	R = 1	2.6 × 106 s

Table 10

Summary of peak times of the Mandel–Cryer effect in the original model, and models considering the smear effect and well resistance.

Model	Location	Peak Time
original	R = 1/8	1.5 × 104 s
original	R = 1/4	3.9 × 104 s
original	R = 1/2	1.8 × 105 s
original	R = 1	5.8 × 105 s
smear effect	R = 1/8	3.0 × 104 s
smear effect	R = 1/4	7.8 × 104 s
smear effect	R = 1/2	2.4 × 105 s
smear effect	R = 1	7.8 × 105 s
well resistance	R = 1/8	1.5 × 104 s
well resistance	R = 1/4	3.9 × 104 s
well resistance	R = 1/2	1.8 × 105 s
well resistance	R = 1	5.8 × 105 s

References

1. Xie, L.; Liang, Z.; Feng, G.; Li, Y.; Wu, T. Improved analytical solution for air–boosted vacuum consolidation of saturated soil using eigenfunction expansion method. Symmetry; 2022; 14, 1757. [DOI: https://dx.doi.org/10.3390/sym14091757]

2. Carrillo, N. Simple two and three dimensional case in the theory of consolidation of soils. J. Math. Phys.; 1942; 21, pp. 1-5. [DOI: https://dx.doi.org/10.1002/sapm19422111]

3. Barron, R.A. Consolidation of fine-grained soils by drain wells by drain wells. Trans. Am. Soc. Civ. Eng.; 1948; 113, pp. 718-742. [DOI: https://dx.doi.org/10.1061/TACEAT.0006098]

4. Horne, M. The consolidation of a stratified soil with vertical and horizontal drainage. Int. J. Mech. Sci.; 1964; 6, pp. 187-197. [DOI: https://dx.doi.org/10.1016/0020-7403(64)90015-3]

5. Yoshikuni, H.; Nakanodo, H. Consolidation of soils by vertical drain wells with finite permeability. Soils Found.; 1974; 14, pp. 35-46. [DOI: https://dx.doi.org/10.3208/sandf1972.14.2_35]

6. Yoshikuni, H. Design and Control of Construction in the Vertical Drain Method; Gihodo Publishing: Tokyo, Japan, 1979.

7. Hansbo, S. Consolidation of fine-grained soils by prefabricated drains. Proceedings of the 10th International Conference on Soil Mechanics and Foundation Engineering; Stockholm, Sweden, 15–19 June 1981; Volume 3, pp. 677-682.

8. Onoue, A. Consolidation by vertical drains taking well resistance and smear into consideration. Soils Found.; 1988; 28, pp. 165-174. [DOI: https://dx.doi.org/10.3208/sandf1972.28.4_165]

9. Xie, K.H.; Zeng, G.X. Consolidation Theories for Drain Wells under Equal Strain Condition. Chin. J. Geotech. Eng.; 1989; 11, pp. 3-17.

10. Basak, P. Consolidation by sand drains in initial gradient soils. J. Geotech. Eng. Div.; 1977; 103, pp. 1327-1332. [DOI: https://dx.doi.org/10.1061/AJGEB6.0000523]

11. Lekha, K.; Krishnaswamy, N.; Basak, P. Consolidation of clay by sand drain under time-dependent loading. J. Geotech. Geoenviron. Eng.; 1998; 124, pp. 91-94. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(1998)124:1(91)]

12. Indraratna, B.; Rujikiatkamjorn, C.; Sathananthan, I. Radial consolidation of clay using compressibility indices and varying horizontal permeability. Can. Geotech. J.; 2005; 42, pp. 1330-1341. [DOI: https://dx.doi.org/10.1139/t05-052]

13. Zhou, Q.; Liu, H.; Chen, Z.B. Non-linear radial consolidation of sand-drained ground with varying volume compressibility and horizontal permeability. Rock Soil Mech.; 2007; 28, pp. 855-858.

14. Guo, B.; Han, Y.; Gong, X.N.; Lu, M.M. Nonlinear consolidation behavior of sand foundation with both horizontal and vertical drainage. J. Shenzhen Univ. (Sci. Eng.); 2010; 27, pp. 459-463.

15. Han, W.; Liu, S.; Zhang, D.W. Nonlinear consolidation of subsoil by vacuum preloading based on bi-ogarithmic coordinate compression model. J. Southeast Univ. (Nat. Sci. Ed.); 2013; 43, pp. 967-972.

16. Gao, G.Y.; Nie, C.X.; Zhang, H.Q. Radial Consolidation Solution of Plastic Wick Drain Combined Vacuum Preloading. J. Tongji Univ. (Nat. Sci.); 2017; 45, pp. 1290-1297.

17. Li, F.; Xie, K.; Deng, Y. Analytical solution for consolidation by vertical drains with exponential flow under vacuum preloading. J. Cent. South Univ. (Sci. Technol.); 2015; 46, pp. 1075-1081.

18. Li, C.X.; Xie, K.H.; Lu, M.M. One-dimensional consolidation analysis considering exponential flow law and time-depending loading. Rock Soil Mech.; 2011; 32, pp. 553-559.

19. Lu, M.M.; Wang, S.; Sloan, S.W.; Indraratna, B.; Xie, K.H. Nonlinear radial consolidation of vertical drains under a general time-variable loading. Int. J. Numer. Anal. Methods Geomech.; 2015; 39, pp. 51-62. [DOI: https://dx.doi.org/10.1002/nag.2295]

20. Lu, M.; Wang, S.; Sloan, S.W.; Sheng, D.; Xie, K. Nonlinear consolidation of vertical drains with coupled radial–vertical flow considering well resistance. Geotext. Geomembr.; 2015; 43, pp. 182-189. [DOI: https://dx.doi.org/10.1016/j.geotexmem.2014.12.001]

21. Schiffman, R.L. Consolidation of soil under time-dependent loading and varying permeability. Proceedings of the Proceedings of the Highway Research Board; Washington, DC, UA, 6–10 January 1958; pp. 584-617.

22. Davis, E.; Raymond, G. A non-linear theory of consolidation. Geotechnique; 1965; 15, pp. 161-173. [DOI: https://dx.doi.org/10.1680/geot.1965.15.2.161]

23. Mesri, G.; Rokhsar, A. Theory of consolidation for clays. J. Geotech. Eng. Div.; 1974; 100, pp. 889-904. [DOI: https://dx.doi.org/10.1061/AJGEB6.0000075]

24. Lekha, K.; Krishnaswamy, N.; Basak, P. Consolidation of clays for variable permeability and compressibility. J. Geotech. Geoenviron. Eng.; 2003; 129, pp. 1001-1009. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2003)129:11(1001)]

25. Zhang, X.; Huang, T.; Wu, Y. Soil drainage clogging mechanism under vacuum preloading: A review. Transp. Geotech.; 2024; 45, 101178. [DOI: https://dx.doi.org/10.1016/j.trgeo.2023.101178]

26. Wang, P.; Wu, J.; Ge, X.; Chen, F.; Yang, X. Non-uniform Consolidation of Soil and Influence of Corresponding Clogging Effect During Vacuum Preloading. Int. J. Geosynth. Ground Eng.; 2022; 8, 58. [DOI: https://dx.doi.org/10.1007/s40891-022-00402-1]

27. Xu, B.H.; He, N.; Jiang, Y.B.; Zhou, Y.Z.; Zhan, X.J. Experimental study on the clogging effect of dredged fill surrounding the PVD under vacuum preloading. Geotext. Geomembr.; 2020; 48, pp. 614-624. [DOI: https://dx.doi.org/10.1016/j.geotexmem.2020.03.007]

28. Hvorslev, M.J. Physical components of the shear strength of saturated clay. Proceedings of the Research Conference on Shear Strength of Cohesive Soils; Boulder, CO, USA, 13–15 June 1960; pp. 169-273.

29. Deng, Y.B.; Xie, K.H.; Wang, K. Mesh sensitivity study for the consolidation of composite ground with granular column by FEM. J. Hunan Univ.; 2012; 39, pp. 12-17.

30. Roscoe, K.H.; Burland, J.B. On the Generalized Stress-Strain Behaviour of Wet Clay; Cambridge University Press: Cambridge, UK, 1968.