The schematic diagram of the one-dimensional non-linear consolidation of a clay layer under an overlying load condition is shown in Figure 1, in which the thickness of the clay layer is l . Using the Lagrangian coordinate system, a (L) is positive vertically downward, and the origin is located at the top of the clay layer.

Gibson et al. (1967); Gibson et al. (1981) gave a governing equation for the non-linear consolidation of large-deformed soil layers (one-dimensional non-linear consolidation theory) governed by the void ratio e , while ignoring the self-weight effect of the soil. ∂ ∂ a k e 1 + e 0 γ w 1 + e d σ ′ d e ∂ e ∂ a + 1 1 + e 0 ∂ e ∂ t = 0 (1) where e is the void ratio (dimensionless), e 0 is the initial void ratio (dimensionless), k e is the hydraulic conductivity (LT^-1), which is the function of e ; σ ′ is the effective stress (ML^-1T^-2), t is the time (T), and γ w is the unit weight of water (ML^-2T^-2).

The one-dimensional non-linear consolidation theory is assumed here (Luo et al., 2020): 1) the water flow in the clay layer obeys Darcy’s law while one-dimensional consolidation of the saturated clay layer occurs; 2) the compression of the clay layer equals the cumulative water release from the clay layer per unit area, because the compressibility of soil particles and pore water is significantly less than that of soil skeleton, and the soil particles and pore water are assumed to be incompressible; 3) σ ′ and k e of the clay layer corresponding to the e are in accordance with the compression and permeability models (Eqs 2, 3), respectively (Mitchell, 1993): e = e 0 − C c log 10 σ ′ / σ 0 ′ (2) e = e 0 + C k log 10 k e / k e 0 (3) where C c is the compression index (dimensionless), σ 0 ′ is the initial effective stress (ML^-1T^-2), C k is the penetration index (dimensionless), and k e 0 is the initial hydraulic conductivity; 4) the clay layer undergoes normal consolidation, and the non-linear variations of the soil compressibility during the consolidation period are assumed to obey Eq. 4 (Wu et al., 2010), and the creep effect is not considered here. 1 + e 0 − C c log 10 σ ′ / σ 0 ′ = 1 + e 0 + e f / 2 (4) where e f is the final void ratio, which can be calculated by Eq. 2 using the ultimate effective stress σ f ′ at the end of clay layer consolidation, and σ f ′ is equal to the sum of σ 0 ′ and additional effective stress caused by the overlying load above the clay layer.

Furthermore, the principle of effective stress is observed in the clay layer, σ ′ = σ 0 ′ + p − u (5) where p is the load due to some overlying structure or layers (ML^-1T^-2), and u a , t is the excess pore water pressure at position a in the clay layer at time t (ML^-1T^-2).

It is assumed that the volumetric compressibility and hydraulic conductivity of the clay layer decrease proportionally (i.e., C k = C c ) during consolidation and the corresponding consolidation process (Davis and Raymond, 1965; Xie et al., 2002; Wu et al., 2010; Luo et al., 2020). The initial effective stress ( σ 0 ′ , geostatic stress) is distributed uniformly through the entire thickness and has little effect on a thin clay layer, and the calculation error of the water released from a thick clay layer will be greater (Gibson et al., 1967; Gibson et al., 1981). Substituting Eqs 2–5 into Eq. 1, the governing equation for the dissipation of excess pore water pressure through a clay layer undergoing one-dimensional non-linear consolidation can be obtained, namely, c w ∂ 2 u ∂ a 2 + 1 σ 0 ′ + p − u ∂ u ∂ a 2 = ∂ u ∂ t (6)

Eq. 6 uses excess pore water pressure u as the primary variable, and c w = 2 ln 10 k e 0 σ 0 ′ 1 + e 0 2 γ w C c 1 + e 0 + e f / 2 (7) which is closely related to the coefficient of consolidation ( c v ) in Terzaghi’s one-dimensional linear consolidation theory (L²T^-1). In addition, c w relies on the values of k e 0 , σ 0 ′ , C c , e 0 and e f .

Under the conditions of SDCL and DDCL, u in the clay layer can be obtained through the governing equation (Eq. 6) with the following initial and boundary conditions (Eqs 8–10). Here, the overlying load is assumed to be constant. u a , 0 = p   0 < a < l (8) ∂ u ∂ a = 0   o r   u 0 , t = 0   a = 0 , t > 0 (9) u l , t = 0   t > 0 (10)

The variations of excess pore water pressure u a , t in SDCL and DDCL can be obtained using the variable transformation and characteristic function method, respectively, and expressed as follows: u 1 a , t = σ f ′ 1 − σ f ′ σ 0 ′ 2 π ∑ n = 1 ∞ − 1 n n − 1 2 e − n − 1 2 2 π 2 c w t 2 l 2 cos n − 1 2 π a l (11) u 2 a , t = σ f ′ 1 − σ f ′ σ 0 ′ 2 π ∑ n = 1 ∞ ( − 1 n − 1 n e − n 2 π 2 c w t l 2 sin n π a l (12) where σ f ′ = σ 0 ′ + p , the subscripts 1 and 2 of the symbols correspond to the case of single-sided drainage and double-sided drainage, respectively, and the detailed derivations are provided in Appendix A, B, respectively.

The void ratio e a , t variation for SDCL and DDCL can be obtained by substituting Eqs 9, 10, respectively, with the assumption proposed by Mitchell (1993). e 1 a , t = e 0 − C c log 10 σ f ′ σ 0 ′ ( 1 + 2 π ∑ n = 1 ∞ − 1 n n − 1 2 e − n − 1 2 2 π 2 c w t 2 l 2 cos n − 1 2 π a l ) (13) e 2 a , t = e 0 − C c log 10 σ f ′ σ 0 ′ ( 1 + 2 π ∑ n = 1 ∞ ( − 1 n − 1 n e − n 2 π 2 c w t l 2 sin n π a l ) (14)

Similarly, the hydraulic conductivity k a , t variation for SDCL and DDCL can be derived, respectively. k 1 a , t = k e 0 σ 0 ′ σ f ′ 1 + 2 π ∑ n = 1 ∞ − 1 n n − 1 2 e − n − 1 2 2 π 2 c w t l 2 cos n − 1 2 π a l (15) k 2 a , t = k e 0 σ 0 ′ σ f ′ 1 + 2 π ∑ n = 1 ∞ ( − 1 n − 1 n e − n 2 π 2 c w t l 2 sin n π a l (16)

Based on Darcy’s law of seepage, the dimensionless flow velocity q ¯ a ¯ , t ¯ for SDCL and DDCL can be obtained, respectively. q ¯ 1 a ¯ , t ¯ = 2 ∑ n = 1 ∞ − 1 n + 1 e − n − 1 2 2 π 2 t ¯ sin n − 1 2 π a ¯ (17) q ¯ 2 a ¯ , t ¯ = 2 ∑ n = 1 ∞ ( − 1 n − 1 ] e − n 2 π 2 t ¯ cos n π a ¯ (18) where q a , t = q ¯ a ¯ , t ¯ ln 10 σ 0 ′ k e 0 γ w l log 10 σ f ′ σ 0 ′ is the flow velocity at position a and time t in SDCL and DDCL (L³T⁻¹), a ¯ = a / l is the dimensionless form of a , and t ¯ = c w t / l 2 is the dimensionless time.

According to the water balance principle, the dimensionless rate of water release for the SDCL and DDCL per unit horizontal area q ¯ t ¯ = q ¯ 1 , t ¯ − q ¯ 0 , t ¯ , respectively, where q ¯ 0 , t ¯ and q ¯ 1 , t ¯ represent the flow velocity at the upper and lower clay layer interfaces, respectively. q ¯ 1 t ¯ = 2 ∑ n = 1 ∞ e − n − 1 2 2 π 2 t ¯ (19) q ¯ 2 t ¯ = 4 ∑ n = 1 ∞ 1 − − 1 n ] e − n 2 π 2 t ¯ (20)

By integrating Eqs 17, 18 over time, the cumulative water fluxes per unit horizontal area of SDCL and DDCL can be obtained, respectively. Q ¯ 1 a ¯ , t ¯ = 2 π 2 ∑ n = 1 ∞ 1 n − 1 2 2 e − n − 1 2 2 π 2 t ¯ − 1 sin n − 1 2 π a ¯ (21) Q ¯ 2 a ¯ , t ¯ = 2 π 2 ∑ n = 1 ∞ 1 − − 1 n n 2 e − n 2 π 2 t ¯ − 1 cos n π a ¯ (22)

Owing to the above assumption about the incompressibility of soil particles, the accumulated released water per unit area of SDCL and DDCL Q ¯ t ¯ = Q ¯ 1 , t ¯ − Q ¯ 0 , t ¯ , which is equal to the deformation of the clay layer. Q ¯ 1 t ¯ = 2 π 2 ∑ n = 1 ∞ 1 n − 1 2 2 1 − e − n − 1 2 2 π 2 t ¯ (23) Q ¯ 2 t ¯ = 4 π 2 ∑ n = 1 ∞ 1 − − 1 n n 2 1 − e − n 2 π 2 t ¯ (24) where Q t = Q ¯ t ¯ C c l f e 0 + f e f 2 1 + e 0 2 log 10 σ f ′ σ 0 ′ is the deformation of SDCL and DDCL (L). The analytical solutions proposed above are compared with the solution based on one-dimensional linear consolidation theory, where the former considers the non-linear variation of compressibility and permeability during the consolidation of the clay layer (Mitchell, 1993). Therefore, the proposed solution is more practical and accurate in estimating deformation and consolidation problems than the conventional linear theory, especially for newly formed soft sedimentary soils.

Figure 2A shows that the evolution of q ¯ with t ¯ . q ¯ 0 , t ¯ is always equal to 0, which leads to q ¯ 1 , t ¯ = q ¯ t ¯ for SDCL; while the flow values at the upper and lower clay layer interfaces for DDCL are equal, but in the opposite direction. It is obvious that both SDCL and DDCL have the same variation law: 1) In the early stage of the consolidation process, due to the abrupt overlying load at the initial time, u in the clay layer instantaneously reaches the maximum value, and q ¯ t ¯ decreases significantly over time. The value of q ¯ t ¯ of DDCL is greater than that of SDCL. 2) The value of q ¯ t ¯ gradually decreases and tends to 0 in the second stage of the consolidation process, while u gradually dissipates. When q ¯ t ¯ approaches 0, the consolidation of the clay layer is completed.

Figure 2B shows the deformation in SDCL and DDCL. Obviously, the Q ¯   v s   t ¯ curve also has the same variation characteristics as the q ¯   v s   t ¯ curve. The Q ¯ t ¯ of SDCL is equal to Q ¯ 1 , t ¯ , and the Q ¯ t ¯ of DDCL is equal to the sum of Q ¯ 0 , t ¯ and Q ¯ 1 , t ¯ . The deformation of the clay layer increases significantly in the early stage of the consolidation process, and then gradually approaches a stable value. At the end of consolidation, the dimensionless deformation of DDCL is the same as that of SDCL. The time to complete consolidation for SDCL is four times longer than that for DDCL.

According to Eqs 23, 24, the deformation of the clay layer can be expressed as: Q t = Q ¯ t ¯ C c l 1 + e 0 + e f / 2 2 1 + e 0 2 log 10 σ f ′ σ 0 ′ (25) where Q ¯ t ¯ is the dimensionless deformation of the clay layer.

By changing Eq. 25 and dimensionless time ( t ¯ = c w t / l 2 ) into logarithmic forms, the following equations can be obtained: log 10 Q ¯ t ¯ = log 10 Q t + log 10 2 1 + e 0 2 C c l 1 + e 0 + e f / 2 log 10 σ f ′ / σ 0 ′ (26) log 10 t ¯ = log 10 t + log 10 c w l 2 (27)

Eq.(26) and (27) are both decomposed into modes, where the corresponding variable and constant terms are added. Therefore, in the double logarithmic coordinate system, the shapes of the curves presented by the deformation Q and the dimensionless deformation Q ¯ are analogous. However, there is an offset in the position of the two curves in the coordinate system ( c w / l 2 in the horizontal direction; 2 1 + e 0 2 C c l 1 + e 0 + e f / 2 log 10 σ f ′ / σ 0 ′ in the vertical direction), but the order of magnitude of the X and Y-axes corresponds to the same. Therefore, the two curves were translated to enable the data on the two curves to overlap as much as possible, and one overlapping point was selected as the matching point. Substituting the coordinates of the matching point ( Q ; Q ¯ , t ,; t ¯ ) into Eqs 28–30, the following parameters can be determined. c w = t ¯ t · l 2 (28) k e 0 = Q Q ¯ · γ w c w ln 10 σ 0 ′ l log 10 σ f ′ / σ 0 ′ (29) C c = 2 ln 10 k e 0 σ 0 ′ 1 + e 0 2 γ w c w 1 + e 0 + e f / 2 (30)

The single-sided drainage experiment carried out by Zhuo et al. (2022) was designed to test the one-dimensional consolidation of the aquitard under load. The thickness of the clay layer was 20 cm, the initial void ratio was 1.28, the load stage was 0–50–100–200 kPa, and the deformations of the clay layer corresponding to the three stages were 6.9, 6.2, and 6.6 mm, respectively.

In this study, the second-level load condition (50–100 kPa) was taken as an example to determine the hydraulic parameters of SDCL with σ 0 ′ = 50 kPa; l = 19.31 cm at this load stage. The deformation of SDCL and the dimensionless deformation curve with time are plotted in Figure 8. The type curve of Q ¯ t ¯ was then superposed onto this figure, by keeping the axes of the two graphs parallel to each other. The matching point refers to the intersection point of the type curve Q ¯ t ¯ with the data curve of Q t . The coordinates of the matching point were Q = 6.220 mm; Q ¯ = 0.980; t = 48,303 s,; t ¯ = 9.105 in both systems. Substituting these values into Eqs 28–30, the values of c w ; C c ,; k e 0 were calculated to be 7.03×10⁻⁶ m²/s, 0.25, and 6.89×10⁻⁸ m/s, respectively. Using the one-dimensional linear consolidation theory; c v was the same as c w , volume compressibility ( m v ) and hydraulic conductivity were calculated to be 0.66 MPa^-1 and 4.62×10⁻⁸ m/s (Toufigh and Ouria, 2009), respectively. By the analytical solutions of Xie and Leo, (2004), the volume compressibility and hydraulic conductivity were calculated to be 0.66 MPa^-1 and 4.60×10⁻⁸ m/s, respectively. The calculated m v and k values from the analytical solutions of Xie and Leo, (2004) are close to those obtained by linear consolidation theory, and the k value predicted by Eq. 29 is larger than that determined by the other two methods.

Using the hydraulic parameters determined by the proposed type curve approach, the predicted and measured deformations of SDCL were compared under two load stages (50–100–200 kPa), as shown in Figure 9. From Figure 9, when the load is in the range of 50–100 kPa, the Q t values calculated by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions by Xie and Leo, (2004) are consistent. Under this load stage, the calculated Q t values by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions of Xie and Leo, (2004) are 6.27×10⁻³ m, 6.35×10⁻³ m, and 6.22×10⁻³ m, respectively. In contrast, the measured value of Q t from the experiment is 6.23×10⁻³ m. During the load stage of 100–200 kPa, the calculated Q t by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions of Xie and Leo, (2004) are 6.50×10⁻³ m, 1.27×10⁻² m, and 1.15×10⁻² m, respectively, while the measured value obtained in the laboratory test is 6.57×10⁻³ m.

The Q t value predicted by the analytical solutions in this study agrees well with the experimental results, while those predicted by the one-dimensional linear consolidation theory and the analytical solutions of Xie and Leo, (2004) are significantly larger than the experimental results. Among them, the value calculated by the analytical solutions of Xie and Leo, (2004) is smaller than that by the one-dimensional linear consolidation theory.

The consolidation deformation characteristics of DDCL under overlying load conditions were studied by laboratory tests. The steps of the test were as follows: (a) The inner walls of the consolidation vessel were cleaned. The cohesive soil was air-dried, crushed, screened, milled, and dried. After cooling to room temperature, the soil sample was weighed with an electronic balance and put into a consolidation container (20 mm high). The soil sample was compacted in three layers with the corresponding heights of 6.67mm, 13.33mm, and 20mm, respectively (Figures 10A–C). (b) The sample (Figure 10D) was put into the rock vacuum saturated water test device (DP SJA, institute for Beijing Yaou Depeng Technology Co., Ltd., Beijing, China), the saturated water pressure was set to 20 kPa, the saturated pressure time was set to 300 min, the saturated water time was set to 240 min, and the immersion time was set to 2,440 min. After the soil saturation test was complete, the container of soil samples was removed and the soil samples were taken out (Figure 10E). (c) The double-sided drainage consolidation test was performed on the saturated soil samples. A single-lever arm oedometer (Institute for Nanjing Ningxi Soil Instrument Co., Ltd., Nanjing, China), which met the requirements of the standard consolidation test, was used in the test (Figure 10F). The load in the consolidation vessel was gradually increased through the following pressure sequence (e.g., 50, 100, and 200 kPa). The dial indicator reading was recorded during the test until the deflection stabilized. The initial thickness of the test soil sample was 2 cm, and the initial void ratio was 0.63. In order to avoid the influence of the stress state in the early stage caused by the difference in the previous sample preparation, the first load stage (0–50 kPa) was used as the pre-load stage in this test. Therefore, the characteristics of the clay layer, l , σ 0 ′ , and e 0 at the end of the consolidation under the load stage of 0–50 kPa are 1.82 cm, 50 kPa, and 0.54, respectively.

The deformation of DDCL and the dimensionless deformation curve with time are plotted in Figure 11. The coordinates of the matching point for DDCL are Q = 0.632 mm, Q ¯ = 0.980, t = 3,600 s, and t ¯ = 11.218 in both systems. According to the analytical solution proposed in this study, c w = 1.03×10⁻⁶ m²/s, k e 0 = 1.06×10⁻⁸ m/s, and C c = 0.18. In contrast, according to the one-dimensional linear consolidation theory, c v was the same as c w , while m v and k e 0 were calculated to be 0.71 MPa^-1 and 7.32×10⁻⁹ m/s, respectively. By the analytical solutions of Xie and Leo, (2004), the volume compressibility and hydraulic conductivity were calculated to be 0.71 MPa^-1 and 7.3×10⁻⁹ m/s, respectively. The values of m v and k calculated by the analytical solutions of Xie and Leo, (2004) are close to those obtained by the one-dimensional linear consolidation theory, and the k predicted by the analytical solutions in this study has a greater value than those predicted by the other two methods.

Figure 12 compares the deformations of DDCL under two load stages (50–100–200 kPa) by using the hydraulic parameters determined by the type curve approach. Similar to the performance of the SDCL, when the load stage is 50–100 kPa, the calculated Q t values by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions of Xie and Leo, (2004) agree with the experimental results. Under this load stage, the calculated Q t by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions of Xie and Leo, (2004) are 6.03 × 10⁻³ m, 6.45 × 10⁻³ m, and 6.32 × 10⁻³ m, respectively, while the measured value of the experiment is 6.10 × 10⁻³ m. During the load stage of 100–200 kPa, the Q t value predicted by the analytical solutions in this study agrees well with the experimental results, while the Q t values predicted by the one-dimensional linear consolidation theory and the analytical solutions of Xie and Leo, (2004) are significantly larger than that predicted by the analytical solutions in this study. Under the load stage of 100–200 kPa, the calculated Q t values by the analytical solutions in this study, the one-dimensional linear consolidation theory, and the analytical solutions of Xie and Leo, (2004) are 6.23 × 10⁻³ m, 1.29 × 10⁻² m, and 1.20 × 10⁻² m, respectively, while the measured values obtained in the laboratory test is 6.25 × 10⁻⁴ m. Therefore, based on these experimental results, the analytical solutions in this study predict the deformation of the clay layer more accurately than the one-dimensional linear consolidation theory and the analytical solutions of Xie and Leo, (2004).