**2. Improved Double-Layer Soil Consolidation Theory Considering Continuous Drainage Boundary Conditions**

In addition to the load, the same basic assumptions as in the one-dimensional consolidation theory by Terzaghi were made. Equations (1)–(5) are the basic assumptions of one-dimensional consolidation theory by Terzaghi:


roadbed.


(2) Soil particles and water are incompressible.

(5) In the osmotic consolidation, the permeability coefficient and the compression coefficient of the soil are constants. (5) In the osmotic consolidation, the permeability coefficient and the compression coefficient of the soil are constants.

*J. Mar. Sci. Eng.* **2019**, *7*, x FOR PEER REVIEW 4 of 18


Based on the above assumptions, the double-layer soil consolidation equations considering continuous drainage boundary conditions were proposed. Figure 1a presents a simplified diagram of the theoretical calculation of an improved double-layered soil foundation consolidation considering continuous drainage boundary conditions. We take the ground table as the coordinate origin. In the figure, *q*(*t*) is an arbitrary loading function. *h*1, *E*s1, and *k*<sup>1</sup> are the parameters of topsoil. *h*2, *E*s2, and *k*<sup>2</sup> are the parameters of subsoil. The load simplification is applied in two stages. *q*<sup>1</sup> and *q*<sup>2</sup> are the primary and secondary load increments, respectively, as shown in Figure 1b. Based on the above assumptions, the double-layer soil consolidation equations considering continuous drainage boundary conditions were proposed. Figure 1a presents a simplified diagram of the theoretical calculation of an improved double-layered soil foundation consolidation considering continuous drainage boundary conditions. We take the ground table as the coordinate origin. In the figure, *q*(*t*) is an arbitrary loading function. *h*1, *E*s1, and *k*1 are the parameters of topsoil. *h*2, *E*s2, and *k*2 are the parameters of subsoil. The load simplification is applied in two stages. *q*1 and *q*<sup>2</sup> are the primary and secondary load increments, respectively, as shown in Figure 1b.

**Figure 1.** Schematic diagram: (**a**) Double-layer soil foundation calculation model, and (**b**) the curve of loading. **Figure 1.** Schematic diagram: (**a**) Double-layer soil foundation calculation model, and (**b**) the curve of loading.

Therefore, one-dimensional consolidation differential equations for a double-layer soil foundation can be obtained: Therefore, one-dimensional consolidation differential equations for a double-layer soil foundation can be obtained:

$$c\_{v1}\frac{\partial^2 u\_1(z,t)}{\partial z^2} = \frac{\partial u\_1(z,t)}{\partial t} - \frac{\partial q(t)}{\partial t}(0 \le z \le h\_1) \tag{6}$$

$$\begin{aligned} c\_{v2} \frac{\partial^2 u\_2(z, t)}{\partial z^2} &= \frac{\partial u\_2(z, t)}{\partial t} - \frac{\partial q(t)}{\partial t} (h\_1 \le z \le h\_1 + h\_2) \\ \dots &\dots & \dots \end{aligned} \tag{7}$$

2 1 12 2 *z tt* ∂ ∂∂ (7) where *u*1 and *u*2 are the first and second layers of excess pore water pressure, respectively; and *v*<sup>1</sup> *c* where *u*<sup>1</sup> and *u*<sup>2</sup> are the first and second layers of excess pore water pressure, respectively; and *cv*<sup>1</sup> and *cv*<sup>2</sup> are the first and second layers of consolidation coefficients, respectively. The latter are calculated by:

$$c\_{v1} = \frac{k\_1 E\_{s1}}{\gamma\_w},\ c\_{v2} = \frac{k\_2 E\_{s2}}{\gamma\_w}.$$

1 1 1 *s v k E c* γ <sup>=</sup> , 2 2 2 *s v k E c* γ <sup>=</sup> . The pore pressures in Equations (6) and (7) can be expressed by the effective stress and converted into the following equations:

$$\frac{\partial \sigma\_1'(z,t)}{\partial t} = c\_{\upsilon 1} \frac{\partial^2 \sigma\_1'(z,t)}{\partial z^2} (0 \le z \le h\_1) \tag{8}$$

$$\frac{\partial \sigma\_2'(z,t)}{\partial t} = c\_{v2} \frac{\partial^2 \sigma\_2'(z,t)}{\partial z^2} (h\_1 \le z \le h\_1 + h\_2) \tag{9}$$

The initial condition is:

*w*

calculated by:

*w*

$$
\sigma'\_{1,2}(z,0) = 0 \tag{10}
$$

The boundary conditions are:

$$
\sigma\_1'(0, t) = q(t) - q(t) \exp(-bt) \tag{11}
$$

$$
\sigma\_2'(h\_1 + h\_2, t) = q(t) - q(t) \exp(-ct) \tag{12}
$$

These boundary conditions are continuous drainage boundary conditions [35,36], and *b* and *c* are parameters related to soil drainage properties. They are interface parameters that reflect the drainage properties of the top and bottom drainage surfaces of the soil layer [41]. As shown in Table 1, the permeability of the drainage boundary can be controlled by adjusting parameters *b* and *c*.

**Table 1.** Description of parameters *b* (top surface) and *c* (bottom surface).


The continuous condition between layers is:

$$
\sigma\_1'|\_{z=h\_1} = \sigma\_2'|\_{z=h\_1} \tag{13}
$$

The continuous flow condition is:

$$\frac{k\_1}{\gamma\_w} \left. \frac{\partial \sigma\_1'}{\partial z} \right|\_{z=h\_1} = \frac{k\_2}{\gamma\_w} \left. \frac{\partial \sigma\_2'}{\partial z} \right|\_{z=h\_1} \tag{14}$$

Equations (8)–(14) constitute the improved double-layer soil consolidation theory equation considering continuous drainage boundary conditions. A Laplace transform is performed on Equations (8)–(14):

$$c\overline{\sigma\_1'}(z,s) - \overline{\sigma\_1'}(z,0) = c\_{v1}\frac{\partial^2 \overline{\sigma\_1'}(z,s)}{\partial z^2} (0 \le z \le h\_1) \tag{15}$$

$$c\overline{\sigma\_2'}(z,s) - \overline{\sigma\_2'}(z,0) = c\_{v2} \frac{\partial^2 \overline{\sigma\_2'}(z,s)}{\partial z^2} (h\_1 \le z \le h\_1 + h\_2) \tag{16}$$

$$\begin{aligned} \overline{\sigma'\_{1,2}}(z,0) = 0 \end{aligned} \tag{17}$$

$$
\overline{\sigma'\_1}(0, \mathbf{s}) = \overline{q}(\mathbf{s}) - \overline{q}(\mathbf{s} + \mathbf{b}) \tag{18}
$$

$$
\overline{\sigma'\_2}(h\_1 + h\_2, s) = \overline{q}(s) - \overline{q}(s + c) \tag{19}
$$

$$\left. \overline{\sigma'\_{1}} \right|\_{z=h\_{1}} = \left. \overline{\sigma'\_{2}} \right|\_{z=h\_{1}} \tag{20}$$

$$\left. \frac{k\_1}{\gamma\_w} \frac{\overline{\partial \sigma\_1'}}{\overline{\partial z}} \right|\_{z=h\_1} = \left. \frac{k\_2}{\gamma\_w} \frac{\overline{\partial \sigma\_2'}}{\overline{\partial z}} \right|\_{z=h\_1} \tag{21}$$

The general solutions of Equations (8) and (9) obtained by Equations (15)–(17) are as follows in the Laplace transform domain:

$$\overline{\sigma\_1'}(z,s) = A\_{11} \exp(r\_1 z) + A\_{12} \exp(-r\_1 z)(0 \le z \le h\_1) \tag{22}$$

$$\overline{\sigma\_2'}(z,s) = A\_{21} \exp(r\_2 z) + A\_{22} \exp(-r\_2 z)(h\_1 \le z \le h\_1 + h\_2) \tag{23}$$

Here:

*r* 2 1 = *<sup>s</sup> cv*<sup>1</sup> , *r* 2 <sup>2</sup> <sup>=</sup> *<sup>s</sup> cv*<sup>2</sup> . Bands Equations (18)–(21) substituted into Equations (22) and (23) yields the following:

$$\begin{cases} A\_{11} + A\_{12} = \overline{q}(s) - \overline{q}(s+b) \\ A\_{21}e^{r\_2(h\_1+h\_2)} + A\_{22}e^{-r\_2(h\_1+h\_2)} = \overline{q}(s) - \overline{q}(s+c) \\ A\_{11}e^{r\_1h\_1} + A\_{12}e^{-r\_1h\_1} = A\_{21}e^{r\_2h\_1} + A\_{22}e^{-r\_2h\_1} \\ \frac{k\_1}{\gamma\_w} \Big( A\_{11}r\_1e^{r\_1h\_1} + A\_{12}(-r\_1)e^{-r\_1h\_1} \Big) = \frac{k\_2}{\gamma\_w} \Big( A\_{21}r\_2e^{r\_2h\_1} + A\_{22}(-r\_2)e^{-r\_2h\_1} \Big) \end{cases} \tag{24}$$

where:

$$\begin{array}{l} \overline{q}(s) - \overline{q}(s+b) = \int\_{o}^{+\infty} q(t) \left(1 - e^{-bt}\right) e^{-st} dt, \\\overline{q}(s) - \overline{q}(s+c) = \int\_{o}^{+\infty} q(t) \left(1 - e^{-ct}\right) e^{-st} dt. \end{array}$$

*q*(*t*) can be expressed as follows (Figure 1b):

$$q(t) = \begin{cases} \frac{t}{t\_1} q\_1 & (0 \le t \le t\_1) \\ q\_1 & (t\_1 \le t \le t\_2) \\ q\_1 + \frac{t - t\_2}{t\_3 - t\_2} q\_2 & (t\_2 \le t \le t\_3) \\ q\_1 + q\_2 & (t\_3 \le t) \end{cases} \tag{25}$$

According to Equations (24) and (25), *q*(*s*) − *q*(*s* + *b*), *q*(*s*) − *q*(*s* + *c*), *A*11, *A*12, *A*21, *A*<sup>22</sup> can be obtained, see the Appendix A for details.

By putting *A*11, *A*12, *A*21, and *A*<sup>22</sup> into Equations (22) and (23) and then performing Laplace inverse transformation, the solutions of consolidation Equations (6) and (7) can be obtained. However, for the complex Laplace solution, it is difficult to carry out the Laplace inverse transform. Then it needs to be solved by the numerical solution of Laplace inverse transform. According to the prior literature, the Stehfest algorithm has better stability and has the advantage of fewer computational parameters [50]. In this case, the Stehfest algorithm is used to write the corresponding program for numerical inversion of Equations (22) and (23). The Stehfest inversion equation is as follows:

$$f(T) = \frac{\ln 2}{T} \sum\_{i=1}^{N} V\_i \overline{f}(\frac{\ln 2}{T} i) \tag{26}$$

where *f*(*s*) is the Laplace function of *f*(*t*), *f*(*s*) = *L*[ *f*(*t*)] = R ∞ 0 *f*(*t*)*e* <sup>−</sup>*stdt*, and *V<sup>i</sup>* = (−1) *N*/2+*i Min*( P *i*,*N*/2) *k*=[ *<sup>i</sup>*+<sup>1</sup> 2 ] *k N*/2 (2*k*)! (*N*/2−*k*)!*k*!(*k*−1)!(*i*−*k*)!(2*k*−*i*)! , where *N* must be a positive even number. Stehfest [46,47]

recommended taking *N* as between 4 and 32. Through repeated verification and reference to the relevant literature [50], we found that 8 was the best choice for *N*.

The Stehfest algorithm can be used to invert the numerical solution of the improved double-layer soil consolidation equation considering the continuous drainage boundary conditions. Then the total consolidation settlement of the double-layer soil can be calculated as follows:

$$S\_{l} = \int\_{0}^{\mathfrak{h}\_{1}} \frac{\sigma\_{1}'(z, t)}{E\_{s1}} dz + \int\_{\mathfrak{h}\_{1}}^{\mathfrak{h}\_{1} + \mathfrak{h}\_{2}} \frac{\sigma\_{2}'(z, t)}{E\_{s2}} dz \tag{27}$$

The average consolidation degree of the double-layer soil is:

$$\mathcal{U} = \frac{S\_t}{S\_{\infty}} = \frac{S\_t}{\left(\int\_0^{l\_1} \frac{p(t)}{E\_{s1}} dz + \int\_{l\_1}^{l\_1 + l\_2} \frac{p(t)}{E\_{s2}} dz\right)}\tag{28}$$

accuracy.

Degree of consolidation (%)
