Determination of hydraulic parameters of non-linear consolidation clay layers by type curve method

The consolidation of clay layers is of great significance for groundwater environmental protection, groundwater storage utilization, and land subsidence. In this study, the governing equation for the excess pore water pressure during the non-linear consolidation process of clay layers under load conditions is obtained based on the one-dimensional non-linear consolidation theory. Analytical solutions are then derived for clay layers with single or double drainage caused by the dissipation of the excess pore water pressure. With these analytical solutions, the groundwater dynamics and deformation of the clay layer are analyzed. Correspondingly, a type curve method is proposed to calculate the hydraulic parameters of the clay layer through laboratory experiments, which verifies the reliability of the analytical solutions. The study results show that the deformation of the clay layer predicted by the non-linear consolidation theory is smaller than that predicted by the linear consolidation theory. The deformation of the clay layer increases with the increase in the thickness of the clay layer, the compressive index, and the overburden load, while it decreases with the increase in the initial void ratio and the initial effective stress. The stable time, at which the consolidation of the clay layer is completed, increases with the increase in the compression index and the thickness of the clay layer, while it decreases with the increase in the initial void ratio, the initial effective stress, and the initial hydraulic conductivity. It does not vary with the load pressure. Conclusively, the deformation prediction based on the non-linear consolidation theory is more accurate and applicable to further load pressures.

Highlights

• Analytical solutions of deformation for clay layer(s) undergoing nonlinear consolidation are derived.

• A type curve method is proposed to determine the hydraulic parameters of clay layers.

• The deformation of the clay layer increases with the increase in the compression index, thickness, and overburden load, and with the decrease in the initial void ratio and initial effective stress.

1 Introduction

With global population growth and economic development, groundwater pollution and geological disasters induced by human activities are becoming increasingly serious, threatening the sustainable development of human society (Gorelick and Zheng, 2015; Jia et al., 2020; Luo et al., 2020; Shu et al., 2022; Zhang et al., 2023). Clayey soils widely cover the surface of basins, plains, seabed, and groundwater aquifer systems, such as the Yangtze Delta in China and the Virginia coastal plain aquifer system in the U.S. (Konikow and Neuzil, 2007; Guo and Li, 2015). Clay and silt are the main components of aquitards, and the hydraulic characteristics of aquitards are closely related to groundwater environmental protection, groundwater storage utilization, and land subsidence (Li et al., 2017; Li et al., 2019; Guo et al., 2022). In addition, clay soil is a significant component of the liner systems installed at hazardous waste disposal sites (e.g., municipal solid waste landfills and nuclear waste disposal sites), thus it has a positive effect on the protection of the water and soil environment around the site, and is an important component of building foundations, dam seepage control systems, mine tailings dams, etc. (Neuzil, 1986; Xue et al., 2005; Zhou et al., 2013; Nicholls et al., 2021; Yan et al., 2021; Zhu et al., 2021). Therefore, it is important to determine the deformation of the clay layer and its hydraulic parameters for groundwater environmental protection and geological disaster prevention.

Many studies have been conducted on the consolidation and drainage of aquitards (Zhuang et al., 2020; Zhuang et al., 2021). Zhou et al. (2013) derived an analytical solution to describe the change of excess pore water pressure in an aquitard based on Terzaghi’s one-dimensional consolidation theory. The drawdown was considered to be a constant in an adjacent aquifer and was verified by laboratory tests. Li et al. (2016) proposed a method for estimating the water release from an aquitard during the entire mining history or a limited period, then applied the method to the aquifer system in the Su-Xi-Chang area, but the groundwater level data were insufficient in the study. In addition, Li et al. (2018) proposed a governing equation of drawdown variation to describe the one-dimensional large-strain consolidation behavior of an aquitard without creep effects and derived an analytical solution for the deformation of a large-strain aquitard under a sudden drawdown in an adjacent aquifer.

In addition, many scholars in soil mechanics and hydrogeology have studied the consolidation and deformation of clay layers (Chen et al., 2021; Liu et al., 2021; Wang et al., 2021). Xie and Leo, (2004) derived a completely explicit analytical solution for one-dimensional large-strain consolidation in thick and thin soil layers and compared it with the classical small-strain theory. They also extensively studied the non-linear consolidation properties of clays. Davis and Raymond, (1965) pointed out that Terzaghi’s one-dimensional linear consolidation theory cannot accurately predict the dissipation rate of pore pressure. Gibson et al. (1981) proposed a consistent theory for finite-strain consolidation of thick and uniform clay layers under load and found that the traditional consolidation theories may seriously underestimate excess pore water pressure in soft soil layers. Considering the complex non-linearity of materials and the geometric characteristics of the soil consolidation process, Wu et al. (2010) derived an analytical solution and obtained the variation law of soil effective stress and void ratio under different conditions through finite element numerical analysis. Ma et al. (2021) established a settlement calculation model for soft soil foundations, developed a method to determine the deformation modulus of the soil before and after damage under load, and proposed a non-linear settlement calculation method considering the structural characteristics of soft clay. The effect of temperature and chemical properties of pore water on the seepage and compression of cohesive soils was also analyzed (Abdullah et al., 1999; Romero et al., 2001; Deng et al., 2014).

Moreover, research on soil permeability is also in full swing. Zhou et al. (2020) and Zhao et al. (2020) proposed a novel double threshold segmentation algorithm to segment cracks, pores, and grains, and pore-scale variables are defined and extracted from these X-ray CT images to study the geometric characteristics of microstructures of porous geomaterials. Zhao and Zhou, (2020a), Zhao and Zhou, (2020b), Zhao and Zhou, (2022); Zhao, et al. (2022); Zhou and Zhao, (2019) successfully predicted the permeability of soil by establishing a model through X-ray CT imaging test, analyzed the internal microstructures of FC and their effects on the compressive strength and permeability, and studied the hydraulic transport properties and fluid flow in single fractures from different projection directions. Zhao et al. (2021) established a digital quantum mechanism-based neural network (DQNN) to study the permeability using digital porosity, coordination number, and pore network size.

The prediction accuracy of consolidation and deformation processes for an aquitard is obviously affected by the uncertainties and errors in the values of hydraulic properties of the clay layer. Based on laboratory testing data, Gregory et al. (2006) proposed several methods to calculate the compressive index and precompression stress from soil compressive test data. An analytical method was proposed to determine the hydrogeological parameters using the data of deformation for a large-strain aquitard (Li et al., 2018). Given that previous analytical solutions were limited by the assumption of a constant head of an overlying unpressurized aquifer, Zhuang et al. (2015) proposed a type curve method for linear drawdown and deformation cases to estimate the hydrogeological parameters of an aquitard in an aquifer system.

Terzaghi’s one-dimensional linear consolidation (Terzaghi, 1943) problem is a classic problem in soil mechanics and is the most generally used theory for the consolidation of saturated soils. However, saturated clay soils have obvious non-linear trends during consolidation (Gibson et al., 1967; Gibson et al., 1981; Xie and Leo, 2004). Under the condition that the drawdown of adjacent aquifers is constant, Luo et al. (2020) derived an analytical solution to calculate the drawdown and water release of an aquitard undergoing non-linear consolidation described by Gibson et al. (1967) and proposed a type curve approach to determine the hydraulic parameters of the aquitard. Moreover, many researchers have studied the non-linear consolidation properties of clay by means of analytical solutions or combined analytical solutions with laboratory experiments under the premise that the reduction of permeability is proportional to the reduction of compressibility ($C_c=C_k$) (Davis and Raymond, 1965; Xie et al., 2002; Wu et al., 2010; Luo et al., 2020). To the authors’ knowledge, no analytical method has been proposed to determine the hydraulic parameters of clay layers undergoing non-linear consolidation under an overburden load condition. Such a new solution would be useful to analyze the non-linear consolidation phenomenon of clay in practical engineering problems, such as soft soil building foundations and landfill clay liners.

Based on the one-dimensional non-linear consolidation theory developed by (Gibson et al., 1967; Gibson et al., 1981), this study first derived the analytical solutions of the excess pore water pressure and the deformation in accordance with Darcy’s law under constant overburden load. Next, the deformation of clay layers with single or double drainage, i.e., singly (SDCL) or doubly draining clay layer (DDCL), was analyzed, and the factors for the deformation were explored by analyzing the sensitivity of multiple parameters. In addition, through a laboratory single and double drainage consolidation test, a type curve approach based on the analytical solutions of non-linear consolidation theory was developed to determine the hydraulic parameters of clay layers and to verify the theoretical research presented in this study.

2 Mathematical model and analytical solutions

2.1 Governing equation for the consolidation of clay layers

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.

FIGURE 1

FIGURE 1. Schematic diagram of 1-D non-linear consolidation of a clay layer under an overlying load.

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.

$\frac{\partial }{\partial a}\left[\frac{k\left(e\right)\left(1+e_0\right)}{{\gamma }_w\left(1+e\right)}\frac{d{\sigma }^{\prime }}{de}\frac{\partial e}{\partial a}\right]+\frac{1}{\left(1+e_0\right)}\frac{\partial e}{\partial t}=0(1)$

where $e$ is the void ratio (dimensionless), $e_0$ is the initial void ratio (dimensionless), $k\left(e\right)$ is the hydraulic conductivity (LT^-1), which is the function of $e$; ${\sigma }^{\prime }$ is the effective stress (ML^-1T^-2), $t$ is the time (T), and ${\gamma }_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) ${\sigma }^{\prime }$ and $k\left(e\right)$ 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}\left({\sigma }^{\prime }/{\sigma }_0^{\prime }\right)(2)$

$e=e_0+C_k{\log }_{10}\left[k\left(e\right)/k\left(e_0\right)\right](3)$

where $C_c$ is the compression index (dimensionless), ${\sigma }_0^{\prime }$ is the initial effective stress (ML^-1T^-2), $C_k$ is the penetration index (dimensionless), and $k\left(e_0\right)$ 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}\left({\sigma }^{\prime }/{\sigma }_0^{\prime }\right)=1+\left(e_0+e_{\mathrm{f}}\right)/2(4)$

where $e_{\mathrm{f}}$ is the final void ratio, which can be calculated by Eq. 2 using the ultimate effective stress ${\sigma }_{\mathrm{f}}^{\prime }$ at the end of clay layer consolidation, and ${\sigma }_{\mathrm{f}}^{\prime }$ is equal to the sum of ${\sigma }_0^{\prime }$ 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,

${\sigma }^{\prime }={\sigma }_0^{\prime }+p-u(5)$

where $p$ is the load due to some overlying structure or layers (ML^-1T^-2), and $u\left(a,t\right)$ 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 (${\sigma }_0^{\prime }$, 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\left[\frac{{\partial }^{2}u}{\partial a^2}+\frac{1}{{\sigma }_0^{\prime }+p-u}{\left(\frac{\partial u}{\partial a}\right)}^2\right]=\frac{\partial u}{\partial t}(6)$

Eq. 6 uses excess pore water pressure $u$ as the primary variable, and

$c_w=\frac{2\ln 10k\left(e_0\right){\sigma }_0^{\prime }{\left(1+e_0\right)}^2}{{\gamma }_w{C}_c\left[1+\left(e_0+e_{\mathrm{f}}\right)/2\right]}(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\left(e_0\right)$, ${\sigma }_0^{\prime }$, $C_c$, $e_0$ and $e_{\mathrm{f}}$.

2.2 Analytical solutions

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\left(a,0\right)=p0<a<l(8)$

$\frac{\partial u}{\partial a}=0\mathrm{o}\mathrm{r}u\left(0,t\right)=0a=0,t>0(9)$

$u\left(l,t\right)=0t>0(10)$

The variations of excess pore water pressure $u\left(a,t\right)$ in SDCL and DDCL can be obtained using the variable transformation and characteristic function method, respectively, and expressed as follows:

$u_1\left(a,t\right)={\sigma }_{\mathrm{f}}^{\prime }\left[1-{\left(\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}\right)}^{\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{\left(-1\right)}^{\mathrm{n}}}{\left(\mathrm{n}-\frac{1}{2}\right)}{\mathrm{e}}^{-\frac{{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2{c}_{w}t}{{2l}^2}}\cos \frac{\left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }a}{l}}\right](11)$

$u_2\left(a,t\right)={\sigma }_{\mathrm{f}}^{\prime }\left[1-{\left(\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}\right)}^{\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{\left.{\left[(-1\right)}^{\mathrm{n}}-1\right]}{\mathrm{n}}{\mathrm{e}}^{-\frac{{\mathrm{n}}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\sin \frac{\mathrm{n}\mathrm{\pi }a}{l}}\right](12)$

where ${\sigma }_{\mathrm{f}}^{\prime }={\sigma }_0^{\prime }+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\left(a,t\right)$ variation for SDCL and DDCL can be obtained by substituting Eqs 9, 10, respectively, with the assumption proposed by Mitchell (1993).

$e_1\left(a,t\right)=e_0-C_c{\log }_{10}\left(\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}\right)(1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{\left(-1\right)}^{\mathrm{n}}}{\left(\mathrm{n}-\frac{1}{2}\right)}{\mathrm{e}}^{-\frac{{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2{c}_{w}t}{{2l}^2}}\cos \frac{\left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }a}{l})(13)$

$e_2\left(a,t\right)=e_0-C_c{\log }_{10}\left(\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}\right)(1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{\left.{\left[(-1\right)}^{\mathrm{n}}-1\right]}{\mathrm{n}}{\mathrm{e}}^{-\frac{{\mathrm{n}}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\sin \frac{\mathrm{n}\mathrm{\pi }a}{l})(14)$

Similarly, the hydraulic conductivity $k\left(a,t\right)$ variation for SDCL and DDCL can be derived, respectively.

$k_1\left(a,t\right)=k\left(e_0\right){\left(\frac{{\sigma }_0^{\prime }}{{\sigma }_{\mathrm{f}}^{\prime }}\right)}^{\left[1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{\left(-1\right)}^{\mathrm{n}}}{\left(\mathrm{n}-\frac{1}{2}\right)}{\mathrm{e}}^{-\frac{{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\cos \frac{\left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }a}{l}\right]}(15)$

$k_2\left(a,t\right)=k\left(e_0\right){\left(\frac{{\sigma }_0^{\prime }}{{\sigma }_{\mathrm{f}}^{\prime }}\right)}^{\left[1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{\left.{\left[(-1\right)}^{\mathrm{n}}-1\right]}{\mathrm{n}}{\mathrm{e}}^{-\frac{{\mathrm{n}}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\sin \frac{\mathrm{n}\mathrm{\pi }a}{l}\right]}(16)$

Based on Darcy’s law of seepage, the dimensionless flow velocity $\overline{q}\left(\overline{a},\overline{t}\right)$ for SDCL and DDCL can be obtained, respectively.

${\overline{q}}_1\left(\overline{a},\overline{t}\right)=2{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\left(-1\right)}^{\mathrm{n}+1}{\mathrm{e}}^{-{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2\overline{t}}\sin \left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }\overline{a}(17)$

${\overline{q}}_2\left(\overline{a},\overline{t}\right)=2{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\left[(-1\right)}^{\mathrm{n}}-1]{\mathrm{e}}^{-{\mathrm{n}}^2{\mathrm{\pi }}^2\overline{t}}\cos \mathrm{n}\mathrm{\pi }\overline{a}(18)$

where $q\left(a,t\right)=\overline{q}\left(\overline{a},\overline{t}\right)\frac{\ln 10{\sigma }_0^{\prime }k\left(e_0\right)}{{\gamma }_{w}l}{\log }_{10}\left(\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}\right)$ is the flow velocity at position $a$ and time $t$ in SDCL and DDCL (L³T⁻¹), $\overline{a}=a/l$ is the dimensionless form of $a$, and $\overline{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 $\overline{q}\left(\overline{t}\right)=\overline{q}\left(1,\overline{t}\right)-\overline{q}\left(0,\overline{t}\right)$, respectively, where $\overline{q}\left(0,\overline{t}\right)$ and $\overline{q}\left(1,\overline{t}\right)$ represent the flow velocity at the upper and lower clay layer interfaces, respectively.

${\overline{q}}_1\left(\overline{t}\right)=2{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\mathrm{e}}^{-{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2\overline{t}}(19)$

${\overline{q}}_2\left(\overline{t}\right)=4{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }{\left[1-\left(-1\right)\right.}^{\mathrm{n}}]{\mathrm{e}}^{-{\mathrm{n}}^2{\mathrm{\pi }}^2\overline{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.

${\overline{Q}}_1\left(\overline{a},\overline{t}\right)=\frac{2}{{\mathrm{\pi }}^2}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{1}{{\left(\mathrm{n}-\frac{1}{2}\right)}^2}\left({\mathrm{e}}^{-{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2\overline{t}}-1\right)\sin \left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }\overline{a}(21)$

${\overline{Q}}_2\left(\overline{a},\overline{t}\right)=\frac{2}{{\mathrm{\pi }}^2}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{\left[1-{\left(-1\right)}^{\mathrm{n}}\right]}{{\mathrm{n}}^2}\left({\mathrm{e}}^{-{\mathrm{n}}^2{\mathrm{\pi }}^2\overline{t}}-1\right)\cos \mathrm{n}\mathrm{\pi }\overline{a}(22)$

Owing to the above assumption about the incompressibility of soil particles, the accumulated released water per unit area of SDCL and DDCL $\overline{Q}\left(\overline{t}\right)=\overline{Q}\left(1,\overline{t}\right)-\overline{Q}\left(0,\overline{t}\right)$, which is equal to the deformation of the clay layer.

${\overline{Q}}_1\left(\overline{t}\right)=\frac{2}{{\mathrm{\pi }}^2}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{1}{{\left(\mathrm{n}-\frac{1}{2}\right)}^2}\left({1-\mathrm{e}}^{-{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2\overline{t}}\right)(23)$

${\overline{Q}}_2\left(\overline{t}\right)=\frac{4}{{\mathrm{\pi }}^2}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{1-\left(-1\right)}^{\mathrm{n}}}{{\mathrm{n}}^2}\left({1-\mathrm{e}}^{-{\mathrm{n}}^2{\mathrm{\pi }}^2\overline{t}}\right)(24)$

where $Q\left(t\right)=\overline{Q}\left(\overline{t}\right)\frac{C_{c}l\left[f\left(e_0\right)+f\left(e_{\mathrm{f}}\right)\right]}{2{\left(1+e_0\right)}^2}{\log }_{10}\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}$ 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.

2.3 Dimensionless analysis

Figure 2A shows that the evolution of $\overline{q}$ with $\overline{t}$. $\overline{q}\left(0,\overline{t}\right)$ is always equal to 0, which leads to $\overline{q}\left(1,\overline{t}\right)=\overline{q}\left(\overline{t}\right)$ 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 $\overline{q}\left(\overline{t}\right)$ decreases significantly over time. The value of $\overline{q}\left(\overline{t}\right)$ of DDCL is greater than that of SDCL. 2) The value of $\overline{q}\left(\overline{t}\right)$ gradually decreases and tends to 0 in the second stage of the consolidation process, while $u$ gradually dissipates. When $\overline{q}\left(\overline{t}\right)$ approaches 0, the consolidation of the clay layer is completed.

FIGURE 2

FIGURE 2. (A) $\overline{\mathit{q}}$ versus $\overline{\mathit{t}}$ curve of SDCL and DDCL. (B) $\overline{\mathit{Q}}$ versus $\overline{\mathit{t}}$ curve of SDCL and DDCL.

Figure 2B shows the deformation in SDCL and DDCL. Obviously, the $\overline{Q}\mathrm{v}\mathrm{s}\overline{t}$ curve also has the same variation characteristics as the $\overline{q}\mathrm{v}\mathrm{s}\overline{t}$ curve. The $\overline{Q}\left(\overline{t}\right)$ of SDCL is equal to $\overline{Q}\left(1,\overline{t}\right)$, and the $\overline{Q}\left(\overline{t}\right)$ of DDCL is equal to the sum of $\overline{Q}\left(0,\overline{t}\right)$ and $\overline{Q}\left(1,\overline{t}\right)$. 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.

3 Results and discussion

In order to evaluate the above analytical solutions, the values of $u$, $e$, $k,$ and $Q$ of SDCL and DDCL were calculated by case analysis. Then, the values of these parameters based on the one-dimensional linear consolidation theory and the analytical solutions proposed in this study (based on the one-dimensional non-linear consolidation theory) were compared. According to previous studies (Wu et al., 2010; Luo et al., 2020), the main parameters of the clay layer used in the case study are as follows: $C_c$ = $C_k$ = 0.0532, ${\sigma }_0^{\prime }$ = 51 kPa, $e_0$ = 1.101, $k\left(e_0\right)$ = 2.10 × 10⁻⁵ m/day, $p$ = 100 kPa, and $l$ = 10 m. In the one-dimensional linear consolidation theory, the parameters of the clay layer were considered to be constant, such as $e_0$ and $k\left(e_0\right)$, and the coefficient of volume compressibility ($m_v$) of the clay layer was obtained by $m_v=\frac{0.434C_c}{\left(1+e_0\right){\sigma }_0^{\prime }}$.

The values of $u$ at different times ($t$ = 1, 3, 5, 10 years) and positions ($a$ = 1, 3, 5, 7, 9 m) were calculated by the one-dimensional linear consolidation theory, and it was calculated by Eqs 11, 12 for SDCL and DDCL, respectively. The results are shown in Figure 3. Obviously, when $t$ = 0, $u$ is equal to the overlying load $p$. From Figure 3A and Figure 3C, the value of $u$ in SDCL decreases with the increase in $a$ and $t$. From Figure 3B and Figure 3D, the value of $u$ in DDCL decreases with the increase in $t$. Meanwhile, the curve of $u$ is symmetric around $a$ =5 m, which means that the clay layer in the middle is the greatest. It was also found that the value of $u$ predicted by Eqs 11, 12 was greater than that by the linear consolidation theory, which indicated that the dissipation of $u$ predicted by the non-linear consolidation theory was lower. The stable time $t_s=0.75l^2/c_w$ predicted by the two theories for the completion of the consolidation of the clay layer is the same. From Figure 3, $t_s$ is about 21 years for DDCL, which is close to a quarter of that of SDCL. This is because the length of the drainage path in SDCL is $l$, which is approximately twice that of DDCL.

FIGURE 3

FIGURE 3. The $\mathit{u}$ versus $\mathit{a}$ curves at different times $\mathit{t}$ in (A) SDCL and (B) DDCL; the $\mathit{u}$ versus $\mathit{t}$ curves at different positions $\mathit{a}$ in (C) SDCL and (D) DDCL.

The values of $e$ at different times and positions in SDCL and DDCL were calculated, and the results are shown in Figure 4. The results show that in SDCL, $e$ decreases with the increase in $a$, while in DDCL, it is symmetric around $a$ = 5 m. In both SDCL and DDCL, $e$ eventually decreases to a stable value while the consolidation of the clay layer completes. The values of $e_{\mathrm{f}}$ predicted by linear and non-linear consolidation theories are 1.076 and 1.056, respectively, for the whole section of the clay layer. Moreover, the value of $e$ predicted by non-linear consolidation theory is greater than that predicted by linear consolidation theory. The shapes of the graphs in Figure 4 are similar to those in Figure 3, that is, the variation of $e$ predicted by non-linear consolidation theory is similar to that of $u$ in the consolidation process of the clay layer.

FIGURE 4

FIGURE 4. The $\mathit{u}$ versus $\mathit{e}$ curves at different times $\mathit{t}$ in (A) SDCL and (B) DDCL; the $\mathit{e}$ versus $\mathit{t}$ curves at different positions $\mathit{a}$ in (C) SDCL and (D) DDCL.

The spatial and temporal distributions of $k\left(e\right)$ in SDCL and DDCL were calculated by Eqs 15, 16, respectively. During the consolidation process of the clay layer, in SDCL, $k\left(e\right)$ decreases with the increase in $a$ and $t$ (Figure 5A; Figure 5C); while in DDCL, $k\left(e\right)$ has a maximum value in the middle of the clay layer, and the upper and lower parts are symmetrical around $a$ = 5 m (Figure 5B; Figure 5D). For both conditions, $k\left(e\right)$ of the clay layer decreases by 66.2% from the initial value of 7.7 × 10⁻³ m/year to 2.6 × 10⁻³ m/year when the consolidation of the clay layer completes. That is, the variation of $k\left(e\right)$ corresponds to the changes in the void ratio in the consolidation process of the clay layer.

FIGURE 5

FIGURE 5. The $\mathit{k}$ versus $\mathit{a}$ curves at different times $\mathit{t}$ in (A) SDCL and (B) DDCL; the $\mathit{k}$ versus $\mathit{t}$ curves at different positions $\mathit{a}$ in (C) SDCL and (D) DDCL.

The values of $Q$ in SDCL and DDCL were calculated based on linear and non-linear consolidation theories, and the results are shown in Figure 6. During the consolidation process of the clay layer, the deformation gradually increases to a stable value using both theories. In both SDCL and DDCL, the deformation predicted by each theory ends up being equal. Among them, the deformation predicted by Eqs 23, 24 is 0.116 m, which is significantly smaller than the deformation predicted by the one-dimensional linear consolidation theory (0.216 m). The value of $t_s$ predicted by both theories are the same, 21 years for SDCL and 84 years for DDCL.

FIGURE 6

FIGURE 6. Predicted $\mathit{Q}$ values in SDCL and DDCL.

4 Parametric sensitivity analysis

Under overlying load conditions, the deformation [$Q$] of a clay layer depends on the hydraulic [$k\left(e_0\right)$] and soil mechanical parameters [$C_c$, $e_0$, ${\sigma }_0^{\prime }$, $l$, and $p$] of the clay layer. Since the $Q$ values were the same for SDCL and DDCL, we analyzed the parametric sensitivity of $Q$ for DDCL by varying $C_c$, $e_0$, ${\sigma }_0^{\prime }$, $l$, $k\left(e_0\right)$, and $p$ values taken in the above case study one by one through non-linear and linear consolidation theories.

Figure 7 shows that the value of $Q$ predicted by Eq. 24 is smaller than that predicted by the one-dimensional linear consolidation theory. The difference in the final stable deformation calculated by linear and non-linear consolidation theories increases with increasing $C_c$ (Figure 7A), $l$ (Figure 7D), and $p$ (Figure 7F) and decreasing $e_0$ (Figure 7B) and ${\sigma }_0^{\prime }$ (Figure 7C), and it does not depend on the value of $k\left(e_0\right)$ (Figure 7E). The value of $C_c$ has a great effect on the $Q$ value of the clay layer (Figure 7A). The value of $Q$ predicted by Eq. 24 increases non-linearly with the increase in $C_c$, while that predicted by the one-dimensional linear consolidation theory increases linearly with the increase in $C_c$. When $C_c$ = 0.1, 0.2, 0.3, 0.4, and 0.5, the corresponding $Q$ values predicted by Eq. 24 are 0.213, 0.421, 0.623, 0.823, and 1.016 m, respectively; while the corresponding $Q$ values predicted by the one-dimensional linear consolidation theory are 0.405, 0.811, 1.216, 1.621, and 2.027 m, respectively. The predicted value of $t_s$ increases with the increase in $C_c$. When $C_c$ = 0.1, 0.2, 0.3, 0.4, and 0.5, the predicted values of $t_s$ are 24.7, 54.6, 87.5, 140.4, and 142.1 years, respectively.

FIGURE 7

FIGURE 7. The $\mathit{Q}$ of DDCL for different values of (A) ${\mathit{C}}_{\mathit{c}}$, (B) ${\mathit{e}}_{\mathbf{0}}$, (C) ${\mathit{\sigma }}_{\mathbf{0}}^{\prime }$, (D) $\mathit{l}$, (E) $\mathit{k}\left({\mathit{e}}_{\mathbf{0}}\right)$, and (F) $\mathit{p}$.

The impact of $e_0$ on $Q$ is shown in Figure 7B. From the figure, the $Q$ values of the clay layer predicted by Eq. 24 and the linear consolidation theory decrease non-linearly with the increase in $e_0$. When $e_0$ is 0.5, 1, 1.5, and 2, the calculated $Q$ values by Eq. 24 are 0.159, 0.120, 0.096, and 0.080 m, respectively, and the $Q$ values from the one-dimensional linear consolidation theory are 0.302, 0.226, 0.181, and 0.151 m, respectively. The values of $t_s$ increase with the increase in $e_0$. When $e_0$ is 0.5, 1, 1.5, and 2, the obtained $t_s$ values are 33.1, 15.8, 14.4, and 13.2 years in DDCL, respectively.

The effect of ${\sigma }_0^{\prime }$ on $Q$ is shown in Figure 7C. When ${\sigma }_0^{\prime }$ is 10, 20, 50, and 100 kPa, the calculated $Q$ values by Eq. 24 are 0.250, 0.187, 0.115, and 0.073 m, respectively, and the calculated $Q$ values by the one-dimensional linear consolidation theory are 1.10, 0.55, 0.22, and 0.11 m, respectively. The value of $Q$ of the clay layer predicted by Eq.24 decreases non-linearly with the increase in ${\sigma }_0^{\prime }$, and that predicted by the one-dimensional linear consolidation theory decreases linearly with the increase in ${\sigma }_0^{\prime }$. The predicted $t_s$ decreases with the increase of ${\sigma }_0^{\prime }$. When ${\sigma }_0^{\prime }$ is 10, 20, 50, and 100 kPa, the obtained $t_s$ values are 92.8, 33.1, 13.7, and 6.5 years, respectively.

The impact of $l$ on $Q$ is shown in Figure 7D. The $Q$ values predicted by both Eq. 24 and the linear theory increase linearly with the increase in $l$. When $l$ is 5, 10, 15, and 20 m, the $Q$ values predicted by Eq. 24 are 0.057, 0.114, 0.171, and 0.228 m, respectively, and that calculated by the one-dimensional linear consolidation theory are 0.108, 0.216, 0.323, and 0.431 m, respectively. The predicted value of $t_s$ increases with the increase in $l$. When $l$ is 5, 10, 15, and 20 m, the $t_s$ values are 4.2, 22.8, 33.0, and 67.6 years, respectively.

The impact of $k\left(e_0\right)$ on $Q$ is shown in Figure 7E. The $Q$ values obtained by Eq. 24 and the linear consolidation theory for all cases are 0.114 and 0.216 m, respectively. The values of $t_s$ decreases linearly with the increase in $k\left(e_0\right)$. When $k\left(e_0\right)$ values are 3.15 × 10⁻³, 6.31× 10⁻³, 1.58 × 10⁻², and 3.15 × 10⁻² m/year, the $t_s$ values are 41.9, 20.8, 8.4, and 4.1 years, respectively.

The impact of $p$ on $Q$ is shown in Figure 7F. When $p$ is 5, 10, 15, and 20 m, the $Q$ values predicted by Eq. 24 are 0.072, 0.114, 0.144, and 0.167 m, respectively, and those by the one-dimensional linear consolidation theory are 0.108, 0.216, 0.323, and 0.431 m, respectively. The predicted value of $Q$ by Eq.24 increases non-linearly with the increase in $p$, and that by the one-dimensional linear consolidation theory increase linearly with the increase in $p$. The value of $t_s$ is not affected by $p$, and has the constant value of about 21 years by both Eq. 24 and the linear consolidation theory.

5 Laboratory experimental studies

In order to validate the analytical solutions proposed in this study, two laboratory experimental studies were performed for SDCL and DDCL. Besides, a type curve approach was developed to determine the clay layer parameters based on the experimental results.

5.1 Determination of clay layer parameters through type curve analysis

According to Eqs 23, 24, the deformation of the clay layer can be expressed as:

$Q\left(t\right)=\overline{Q}\left(\overline{t}\right)\frac{C_{c}l\left[1+\left(e_0+e_{\mathrm{f}}\right)/2\right]}{2{\left(1+e_0\right)}^2}{\log }_{10}\frac{{\sigma }_{\mathrm{f}}^{\prime }}{{\sigma }_0^{\prime }}(25)$

where $\overline{Q}\left(\overline{t}\right)$ is the dimensionless deformation of the clay layer.

By changing Eq. 25 and dimensionless time ($\overline{t}$ = $c_{w}t/l^2$) into logarithmic forms, the following equations can be obtained:

${\log }_{10}\overline{Q}\left(\overline{t}\right)={\log }_{10}Q\left(t\right)+{\log }_{10}\left[\frac{2{\left(1+e_0\right)}^2}{C_{c}l\left[1+\left(e_0+e_{\mathrm{f}}\right)/2\right]{\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)}\right](26)$

${\log }_{10}\overline{t}={\log }_{10}t+{\log }_{10}\frac{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 $\overline{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; $\frac{2{\left(1+e_0\right)}^2}{C_{c}l\left[1+\left(e_0+e_{\mathrm{f}}\right)/2\right]{\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)}$ 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$; $\overline{Q}$, $t$,; $\overline{t}$) into Eqs 28–30, the following parameters can be determined.

$c_w=\frac{\overline{t}}{t}·l^2(28)$

$k\left(e_0\right)=\frac{Q}{\overline{Q}}·\frac{{\gamma }_w{c}_w}{\ln 10{\sigma }_0^{\prime }l{\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)}(29)$

$C_c=\frac{2\ln 10k\left(e_0\right){\sigma }_0^{\prime }{\left(1+e_0\right)}^2}{{\gamma }_w{c}_w\left[1+\left(e_0+e_{\mathrm{f}}\right)/2\right]}(30)$

5.2 Results of the SDCL experiment

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 ${\sigma }_0^{\prime }$ = 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 $\overline{Q}\left(\overline{t}\right)$ 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 $\overline{Q}\left(\overline{t}\right)$ with the data curve of $Q\left(t\right)$. The coordinates of the matching point were $Q=$ 6.220 mm; $\overline{Q}=$ 0.980; $t=$ 48,303 s,; $\overline{t}=$ 9.105 in both systems. Substituting these values into Eqs 28–30, the values of $c_w$; $C_c$,; $k\left(e_0\right)$ 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.

FIGURE 8

FIGURE 8. Determination of parameters of SDCL by the type curve approach.

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\left(t\right)$ 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\left(t\right)$ 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\left(t\right)$ from the experiment is 6.23×10⁻³ m. During the load stage of 100–200 kPa, the calculated $Q\left(t\right)$ 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.

FIGURE 9

FIGURE 9. Comparison of predicted and measured deformation of SDCL.

The $Q\left(t\right)$ 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.

5.3 Results of the doubly draining clay layer experiment

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$, ${\sigma }_0^{\prime }$, 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.

FIGURE 10

FIGURE 10. DDCL experiment equipment (A) electronic balance, (B) compacter, (C) consolidation container, (D) soil sample to be saturated, (E) rock vacuum saturated water test device, and (F) single-lever arm oedometer.

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, $\overline{Q}$ = 0.980, $t$ = 3,600 s, and $\overline{t}$ = 11.218 in both systems. According to the analytical solution proposed in this study, $c_w$ = 1.03×10⁻⁶ m²/s, $k\left(e_0\right)$ = 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\left(e_0\right)$ 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 11

FIGURE 11. Determination of parameters of DDCL by the type curve approach.

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\left(t\right)$ 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\left(t\right)$ 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\left(t\right)$ value predicted by the analytical solutions in this study agrees well with the experimental results, while the $Q\left(t\right)$ 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\left(t\right)$ 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).

FIGURE 12

FIGURE 12. Comparison of predicted and measured deformation of DDCL.

6 Conclusion

In this study, analytical solutions (ignoring creep effects) for SDCL and DDCL undergoing non-linear consolidation, excess pore water pressure, and deformation under overlying load conditions were derived. Combined with the case analysis based on the analytical solution presented in this paper, it is found that under the condition of an overlying load, the deformation in SDCL is equal to that in DDCL. There is a lag in $t_s$ for SDCL compared to DDCL. The drainage path of the consolidated drainage of SDCL is twice as large as that of DDCL, resulting in a $t_s$ lag that is about four times that of DDCL. In addition, there is a difference in $u$ under the two theories, i.e., the $u$ value by the non-linear consolidation theory is larger than that by the linear consolidation theory. The deformation predicted by the non-linear consolidation theory is smaller than that by the linear consolidation theory. Larger $C_c$, $l$, and $p$ and smaller $e_0$ and ${\sigma }_0^{\prime }$ lead to a larger difference in the amount of deformation under the two theories. The $t_s$ predicted by the two theories are the same. Through the sensitivity analysis of the parameters, it is obvious that the deformation in SDCL and DDCL is controlled by $C_c$, $e_0$, ${\sigma }_0^{\prime }$, $l$, and $p$. However, the variation of $k\left(e_0\right)$ does not affect the deformation in SDCL and DDCL, and $t_s$ does not depend on the overlying load of the clay layer. Besides, a type curve approach is proposed in this study based on the non-linear consolidation theory, and it agrees well with the measured data from the test. Compared with the linear consolidation theory and the analytical solution of Xie and Leo, (2004), the analytical solution based on the non-linear consolidation theory in this study is more accurate and suitable for predicting the consolidation and deformation of the clay layer under further load conditions.

7 Open research

Parameter values in Sections 2, 3, and 4 are available through Wu et al. (2010) (https://doi.org/10.16285/j.rsm.2010.01.045) and Li et al. (2017) (https://doi.org/10.1002/hyp.11373). The test data in Section 5.2 are available through Zhuo et al. (2022) (https://doi.org/10.1016/j.jhydrol.2022.127994). The test data in Section 5.3 are available through the repository of Dryad (https://datadryad.org/stash/share/frvva1Jd1ETIo-VSViK4-ChOTT6hW9i9EjP2Foj1ExA) and the supporting information files.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

RW is responsible for manuscript preparation and data collection; ZL is responsible for the preparation, data correction and approval of manuscripts, as well as providing fund support; MX, QZ, WI, and HL are responsible for the correction and approval of manuscripts.

Funding

This study was financed by the Open Fund of State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology) (SKLGP2021K023), Sichuan Provincial Science and Technology Department Youth Fund (2023NSFSC0802), the National Natural Science Foundation of China (41702253) and the China Scholarship Council for working at the University of Waterloo as a visiting scholar. WA. I appreciates the support of the Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (NSERC), which made this collaboration possible.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1131128/full#supplementary-material

References

Abdullah, W. S., Alshibli, K. A., and Al-Zou'bi, M. S. (1999). Influence of pore water chemistry on the swelling behavior of compacted clays. Appl. Clay Sci. 15 (5-6), 447–462. doi:10.1016/S0169-1317(99)00034-4

CrossRef Full Text | Google Scholar

Chen, P. S., Li, J. C., Huang, M. H., and Li, D. J. (2021). Consolidation of viscoelastic soil with vertical drains for continuous drainage boundary conditions incorporating a fractional derivative model. Front. Mater. 8, 670150. doi:10.3389/fmats.2021.670150

CrossRef Full Text | Google Scholar

Davis, E. H., and Raymond, G. P. (1965). A non-linear theory of consolidation. Geotechnique 15 (2), 161–173. doi:10.1680/geot.1965.15.2.161

CrossRef Full Text | Google Scholar

Deng, Y. F., Yue, X. B., Cui, Y. J., Shao, G. H., Liu, S. Y., and Zhang, D. W. (2014). Effect of pore water chemistry on the hydro-mechanical behaviour of Lianyungang soft marine clay. Appl. Clay Sci. 95, 167–175. doi:10.1016/j.clay.2014.04.007

CrossRef Full Text | Google Scholar

Gibson, R. E., England, G. L., and Hussey, M. J. L. (1967). The theory of one-dimensional consolidation of saturated clays. Geotechnique 17 (3), 261–273. doi:10.1680/geot.1967.17.3.261

CrossRef Full Text | Google Scholar

Gibson, R. E., Schiffman, R. L., and Cargill, K. W. (1981). The theory of one-dimensional consolidation of saturated clays. II. Finite nonlinear consolidation of thick homogeneous layers. Geotechnique 18, 280–293. doi:10.1139/t81-030

CrossRef Full Text | Google Scholar

Gorelick, S. M., and Zheng, C. (2015). Global change and the groundwater management challenge. Water Resour. Res. 51 (5), 3031–3051. doi:10.1002/2014wr016825

CrossRef Full Text | Google Scholar

Gregory, A. S., Whalley, W. R., Watts, C. W., Bird, N. R. A., Hallett, P. D., and Whitmore, A. P. (2006). Calculation of the compression index and precompression stress from soil compression test data. Soil Tillage Res. 89 (1), 45–57. doi:10.1016/j.still.2005.06.012

CrossRef Full Text | Google Scholar

Guo, H. P., Hao, A. B., Li, W. P., Zang, X. S., Wang, Y. L., Zhu, J. Y., et al. (2022). Land subsidence and its affecting factors in cangzhou, north China plain. Front. Environ. Sci. 10. doi:10.3389/fenvs.2022.1053362

CrossRef Full Text | Google Scholar

Guo, Q. N., and Li, H. L. (2015). Terrestrial-originated submarine groundwater discharge through deep multilayered aquifer systems beneath the seafloor. Hydrol. Process. 29 (4), 295–309. doi:10.1002/hyp.10163

CrossRef Full Text | Google Scholar

Jia, X. Y., Hou, D. Y., Wang, L. W., O'Commor, D., and Luo, J. (2020). The development of groundwater research in the past 40 years: A burgeoning trend in groundwater depletion and sustainable management. J. Hydrology 587, 125006. doi:10.1016/j.jhydrol.2020.125006

CrossRef Full Text | Google Scholar

Konikow, L. F., and Neuzil, C. E. (2007). A method to estimate groundwater depletion from confining layers. Water Resour. Res. 43 (7), W07417. doi:10.1029/2006wr005597

CrossRef Full Text | Google Scholar

Li, Z. F., Zhou, Z. F., Chen, Z., Liu, G. Q., and Zhou, C. Y. (2016). An analytical method to estimate groundwater depletion from a confining layer. Nat. Hazards 85 (2), 887–901. doi:10.1007/s11069-016-2609-1

CrossRef Full Text | Google Scholar

Li, Z. F., Zhou, Z. F., Dai, Y. F., and Dai, B. B. (2019). Contaminant transport in a largely-deformed aquitard affected by delayed drainage. J. Contam. Hydrology 221, 118–126. doi:10.1016/j.jconhyd.2019.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z. F., Zhou, Z. F., Dai, Y. F., Huang, Y., Guo, Q. N., and Zhou, C. Y. (2017). Contaminant transport in a small deformation aquitard affected by the delayed drainage phenomenon. Hydrol. Process. 31 (24), 4429–4439. doi:10.1002/hyp.11373

CrossRef Full Text | Google Scholar

Li, Z. F., Zhou, Z. F., Li, M. Y., Zhang, B. R., and Dai, B. B. (2018). Delayed drainage of a largely deformed aquitard due to abrupt water head decline in adjacent aquifer. Geofluids 2018, 1–12. doi:10.1155/2018/2326491

CrossRef Full Text | Google Scholar

Liu, C., Shi, B., Gu, K., Zhang, T. S., Tang, C. S., Wang, Y., et al. (2021). Negative pore water pressure in aquitard enhances land subsidence: Field, laboratory, and numerical evidence. Water Resour. Res. 58 (1), e2021WR030085. doi:10.1029/2021wr030085

CrossRef Full Text | Google Scholar

Luo, H., Li, Z. F., Luo, M. L., Zhang, Q., and Illman, W. A. (2020). An analytical method to calculate groundwater released from an aquitard undergoing nonlinear consolidation. Water Resour. Res. 56 (9), e2020WR027320. doi:10.1029/2020wr027320

CrossRef Full Text | Google Scholar

Ma, B. H., Li, Z., Cai, K., Hu, Z. Y., Zhao, M. H., He, C. B., et al. (2021). An improved nonlinear settlement calculation method for soft clay considering structural characteristics. Geofluids 2021, 1–7. doi:10.1155/2021/8837889

CrossRef Full Text | Google Scholar

Mitchell, J. K. (1993). Fundamentals of soil behavior. New York: Wiley.

Google Scholar

Neuzil, C. E. (1986). Groundwater flow in low-permeability environments. Water Resour. Res. 22 (8), 1163–1195. doi:10.1029/wr022i008p01163

CrossRef Full Text | Google Scholar

Nicholls, R. J., Beaven, R. P., Stringfellow, A., Monfort, D., Le Cozannet, G., Wahl, T., et al. (2021). Coastal landfills and rising sea levels: A challenge for the 21st century. Front. Mar. Sci. 8. doi:10.3389/fmars.2021.710342

CrossRef Full Text | Google Scholar

Romero, E., Gens, A., and Lloret, A. (2001). Temperature effects on the hydraulic behaviour of an unsaturated clay. Geotechnical Geol. Eng. 19, 311–332. doi:10.1023/A:1013133809333

CrossRef Full Text | Google Scholar

Shu, B., Chen, Y., Amani-Beni, M., and Zhang, R. Z. (2022). Spatial distribution and influencing factors of mountainous geological disasters in southwest China: A fine-scale multi-type assessment. Front. Environ. Sci. 10, 1049333. doi:10.3389/fenvs.2022.1049333

CrossRef Full Text | Google Scholar

Terzaghi, K. (1943). Theoretical soil mechanics. New York: Wiley. doi:10.1002/9780470172766

CrossRef Full Text | Google Scholar

Toufigh, M. M., and Ouria, A. (2009). Consolidation of inelastic clays under rectangular cyclic loading. Soil Dyn. Earthq. Eng. 29 (2), 356–363. doi:10.1016/j.soildyn.2008.03.006

CrossRef Full Text | Google Scholar

Wang, F. Y., Huang, H. W., Yin, Z. Y., and Huang, Q. (2021). Probabilistic characteristics analysis for the time-dependent deformation of clay soils due to spatial variability. Eur. J. Environ. Civ. Eng. 26, 6096–6114. doi:10.1080/19648189.2021.1933604

CrossRef Full Text | Google Scholar

Wu, J., Xie, X. Y., and Zhu, X. R. (2010). Study of properties of 1-D complex nonlinear consolidation of saturated soils (in Chinese). Rock Soil Mech. 31 (1), 81–86. doi:10.16285/j.rsm.2010.01.045

CrossRef Full Text | Google Scholar

Xie, K. H., and Leo, C. J. (2004). Analytical solutions of one-dimensional large strain consolidation of saturated and homogeneous clays. Comput. Geotechnics 31 (4), 301–314. doi:10.1016/j.compgeo.2004.02.006

CrossRef Full Text | Google Scholar

Xie, K. H., Xie, X. Y., and Jiang, W. (2002). A study on one-dimensional nonlinear consolidation of double-layered soil. Comput. Geotechnics 29 (2), 151–168. doi:10.1016/s0266-352x(01)00017-9

CrossRef Full Text | Google Scholar

Xue, Y. Q., Zhang, Y., Ye, S. J., Wu, J. C., and Li, Q. F. (2005). Land subsidence in China. Environ. Geol. 48 (6), 713–720. doi:10.1007/s00254-005-0010-6

CrossRef Full Text | Google Scholar

Yan, H. X., Xie, H. J., Wu, J. W., Ding, H., Qiu, Z. H., and Sun, Z. L. (2021). Analytical model for transient coupled consolidation and contaminant transport in landfill liner system. Comput. Geotechnics 138, 104345. doi:10.1016/j.compgeo.2021.104345

CrossRef Full Text | Google Scholar

Zhang, Y. H., Huo, X. F., and Luo, Y. (2023). Prediction of groundwater pollution diffusion path based on multi-source data fusion. Front. Environ. Sci. 10, 1116309. doi:10.3389/fenvs.2022.1116309

CrossRef Full Text | Google Scholar

Zhao, Z., Zhou, X. P., and Chen, J. W. (2022). Pore-scale hydraulic properties of virtual sandstone microstructures: Spatial variations and voxel scale effects. Archives Civ. Mech. Eng. 23 (1), 22. doi:10.1007/s43452-022-00566-7

CrossRef Full Text | Google Scholar

Zhao, Z., and Zhou, X. P. (2022). Digital analysis for pore-scale compressive strength and permeability of foamed cement with realistic microstructures by X-ray-μCT imaging. Constr. Build. Mater. 346, 128456. doi:10.1016/j.conbuildmat.2022.128456

CrossRef Full Text | Google Scholar

Zhao, Z., and Zhou, X. P. (2020a). Digital voxel-based fracture extraction: Insights to characterization of single fracture flow and anisotropy permeability. J. Nat. Gas Sci. Eng. 84, 103635. doi:10.1016/j.jngse.2020.103635

CrossRef Full Text | Google Scholar

Zhao, Z., and Zhou, X. P. (2020b). Pore-scale effect on the hydrate variation and flow behaviors in microstructures using X-ray CT imaging. J. Hydrology 584, 124678. doi:10.1016/j.jhydrol.2020.124678

CrossRef Full Text | Google Scholar

Zhao, Z., Zhou, X. P., and Qian, Q. H. (2021). Dqnn: Pore-scale variables-based digital permeability assessment of carbonates using quantum mechanism-based machine-learning. Sci. China Technol. Sci. 65 (2), 458–469. doi:10.1007/s11431-021-1906-1

CrossRef Full Text | Google Scholar

Zhao, Z., Zhou, X. P., and Qian, Q. H. (2020). Fracture characterization and permeability prediction by pore scale variables extracted from X-ray CT images of porous geomaterials. Sci. China Technol. Sci. 63 (5), 755–767. doi:10.1007/s11431-019-1449-4

CrossRef Full Text | Google Scholar

Zhou, X. P., and Zhao, Z. (2019). Digital evaluation of nanoscale-pore shale fractal dimension with microstructural insights into shale permeability. J. Nat. Gas Sci. Eng. 75, 103137. doi:10.1016/j.jngse.2019.103137

CrossRef Full Text | Google Scholar

Zhou, X. P., Zhao, Z., and Li, Z. (2020). Cracking behaviors and hydraulic properties evaluation based on fractural microstructure models in geomaterials. Int. J. Rock Mech. Min. Sci. 130, 104304. doi:10.1016/j.ijrmms.2020.104304

CrossRef Full Text | Google Scholar

Zhou, Z. F., Guo, Q. N., and Dou, Z. (2013). Delayed drainage of aquitard in response to sudden change in groundwater level in adjacent confined aquifer: Analytical and experimental studies. Chin. Sci. Bull. 58 (25), 3060–3069. doi:10.1007/s11434-013-5730-5

CrossRef Full Text | Google Scholar

Zhu, W., Yang, X., He, J., Wang, X., Lu, R., and Zhang, Z. (2021). Investigation and systematic risk assessment in a typical contaminated site of hazardous waste treatment and disposal. Front. Public Health 9, 764788. doi:10.3389/fpubh.2021.764788

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuang, C., Yan, L., Zhou, Z. F., Wang, J. G., and Dou, Z. (2021). Estimation of aquitard hydraulic conductivity and skeletal specific storage considering non-Darcy flow. Water Sci. Eng. 14 (4), 269–276. doi:10.1016/j.wse.2021.09.003

CrossRef Full Text | Google Scholar

Zhuang, C., Zhou, Z. F., Illman, W. A., Dou, Z., and Wang, J. G. (2020). Parameter estimation of an overconsolidated aquitard subjected to periodic hydraulic head variations within adjacent aquifers. J. Hydrology 583, 124555. doi:10.1016/j.jhydrol.2020.124555

CrossRef Full Text | Google Scholar

Zhuang, C., Zhou, Z. F., Zhan, H. B., and Wang, G. Y. (2015). A new type curve method for estimating aquitard hydraulic parameters in a multi-layered aquifer system. J. Hydrology 527, 212–220. doi:10.1016/j.jhydrol.2015.04.062

CrossRef Full Text | Google Scholar

Zhuo, Y., Wang, J. G., Zheng, S. Y., and Li, W. J. (2022). Experimental study of one-dimensional nonlinear consolidation of an aquitard under multistage loading. J. Hydrology 611, 127994. doi:10.1016/j.jhydrol.2022.127994

CrossRef Full Text | Google Scholar

Appendix A Derivation of Eq.11

Using the substitution method:

$w={\log }_{10}\left({\sigma }^{\prime }/{\sigma }_0^{\prime }\right)={\log }_{10}\left(\left({\sigma }_0^{\prime }+p-u\right)/{\sigma }_0^{\prime }\right)$(A-1)

Substituting Eq. A-1 into –Eq.10, can simplify the original equation to $w$ as the independent variable:

$c_w\frac{{\partial }^{2}w}{\partial a^2}=\frac{\partial w}{\partial t}$(A-2)

$w\left(a,0\right)=00<a<l$(A-3)

$\frac{\partial w}{\partial a}=0a=0,t>0$(A-4)

$w\left(l,t\right)={\log }_{10}\left(\left({\sigma }_0^{\prime }+p\right)/{\sigma }_0^{\prime }\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)t>0$(A-5)

By applying the boundary subordination processing to Eq. A-5:

$w\left(a,t\right)=p\left(a,t\right)+q\left(a,t\right)$(A-6)

$q\left(a,t\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)$(A-7)

The governing equations and boundary conditions for $p\left(a,t\right)$ as the independent variable can be obtained as:

$c_w\frac{{\partial }^{2}p}{\partial a^2}=\frac{\partial p}{\partial t}$(A-8)

$p\left(a,0\right)=-{\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)0<a<l$(A-9)

$\frac{\partial p}{\partial a}=0a=0,t>0$(A-10)

$p\left(l,t\right)=0t>0$(A-11)

Using the separation of variables method, the solution to Eq. A-8–Eq. A-11 can be derived as:

$p\left(a,t\right)=\frac{2{\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)}{\pi }{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{\left(-1\right)}^{\mathrm{n}}}{\left(\mathrm{n}-\frac{1}{2}\right)}{\mathrm{e}}^{-\frac{{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\cos \frac{\left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }a}{l}$(A-12)

Combining Eq. A-7 with Eq. A-12, the solution to Eq. A-2 can be given,

$w\left(a,t\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)\left(1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{{\left(-1\right)}^{\mathrm{n}}}{\left(\mathrm{n}-\frac{1}{2}\right)}{\mathrm{e}}^{-\frac{{\left(\mathrm{n}-\frac{1}{2}\right)}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\cos \frac{\left(\mathrm{n}-\frac{1}{2}\right)\mathrm{\pi }a}{l}\right)$(A-13)

Finally, the solution to Eq. 11 is obtained by combining Eq. A-1 and Eq. A-13.

Appendix B Derivation of Eq. 12

As the variable passing method in Appendix A, the independent variable $u$ is replaced by $w$.

$c_w\frac{{\partial }^{2}w}{\partial a^2}=\frac{\partial w}{\partial t}$(B-1)

$w\left(a,0\right)=00<a<l$(B-2)

$w\left(0,t\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)t>0$(B-3)

$w\left(l,t\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)t>0$(B-4)

The corresponding solution procedure is the same as given in Appendix A, and the solution to Eq. B-1–Eq. B-4 can be obtained as:

$w\left(a,t\right)={\log }_{10}\left({\sigma }_{\mathrm{f}}^{\prime }/{\sigma }_0^{\prime }\right)(1+\frac{2}{\mathrm{\pi }}{\displaystyle \sum }_{\mathrm{n}=1}^{\infty }\frac{\left.{\left[(-1\right)}^{\mathrm{n}}-1\right]}{\mathrm{n}}{\mathrm{e}}^{-\frac{{\mathrm{n}}^2{\mathrm{\pi }}^2{c}_{w}t}{l^2}}\sin \frac{\mathrm{n}\mathrm{\pi }a}{l})$(B-5)

Finally, the solution to Eq. 12 is obtained by combining Eq. A-1 and Eq. B-5.