## 4 Solution Derivation  
The solution to Equation (21), based on eigenfunction expansion, can be expressed as  
\[u(z,t)=\sum_{m=1}^{\infty}T_{m}(t)\sin\bigg{(}\frac{M}{H}z\bigg{)} \tag{23}\]  
where \(M=\frac{1}{2}(2m-1)\pi,m=1,2,3\cdots\), and \(T_{m}(t)\) is the generalized Fourier coefficient that varies with time \(t\). Equation (23) satisfies the boundary conditions given in Equation (10). By substituting Equation (23) into Equation (21), a family of ordinary integro-differential equations can be obtained as  
\[\bigg{(}\frac{c_{\
u}M^{2}}{H^{2}}+\frac{\alpha_{0}}{\lambda_{0}}\bigg{)}T_{m }(t)+\bigg{(}2-\alpha_{0}+\frac{1-\alpha_{1}}{2-\alpha_{1}}\frac{1}{\kappa} \bigg{)}\frac{dT_{m}(t)}{dt}+\frac{\alpha_{1}}{\left(2-\alpha_{1}\right)^{2} \kappa\lambda_{1}}\int_{0}^{t}\frac{dT_{m}(\tau)}{d\tau}e^{-\frac{\alpha_{1}} {2-\alpha_{1}}\cdot\frac{t-\tau}{\lambda_{1}}}d\tau+C_{m}f_{q}(t)=0 \tag{24}\]  
where  
\[C_{m}=\frac{2}{M}\bigg{[}1+(-1)^{m+1}\frac{\zeta}{M}\bigg{]} \tag{25}\]  
Considering the initial condition given in Equation (11), the Laplace transform of Equation (24) yields  
\[\bigg{[}\frac{c_{\
u}M^{2}}{H^{2}}+\frac{\alpha_{0}}{\lambda_{0}}+(2-\alpha_{ 0})s+\frac{(1-\alpha_{1})\lambda_{1}s+\alpha_{1}}{(2-\alpha_{1})\lambda_{1}s+ \alpha_{1}}\frac{s}{\kappa}\bigg{]}\overline{T}_{m}(s)-\bigg{[}\frac{\alpha_ {0}}{\lambda_{0}}+(2-\alpha_{0})s+\frac{(1-\alpha_{1})\lambda_{1}s+\alpha_{1} }{(2-\alpha_{1})\lambda_{1}s+\alpha_{1}}\frac{s}{\kappa}\bigg{]}C_{m}Q(s)=0 \tag{26}\]  
where \(\overline{T}_{m}(s)=L[T_{m}(t)]\) and \(Q(s)=L[q(t)]\). Solving for \(\overline{T}_{m}(s)\) in Equation (26) results in  
\[\overline{T}_{m}(s)=\Bigg{\{}1-\frac{c_{\
u}M^{2}}{H^{2}}\bigg{[}\frac{c_{\
u }M^{2}}{H^{2}}+\frac{\alpha_{0}}{\lambda_{0}}+(2-\alpha_{0})s+\frac{(1-\alpha _{1})\lambda_{1}s+\alpha_{1}}{(2-\alpha_{1})\lambda_{1}s+\alpha_{1}}\frac{s}{ \kappa}\bigg{]}^{-1}\Bigg{\}}\cdot C_{m}Q(s) \tag{27}\]  
The inverse Laplace transform into Equation (27) leads to the general solution of Equation (24). The result is  
\[T_{m}(t)=C_{m}\Bigg{[}q(T_{\
u})-\int_{0}^{T_{\
u}}q(T_{\tau})\sum_{j=1}^{2}A_ {j}e^{x_{j}(T_{\
u}-T_{\tau})}dT_{\tau}\Bigg{]} \tag{28}\]  
where  
\[x_{1,2}=-\frac{\left(\kappa M^{2}T_{\lambda_{1}}+\alpha_{0}\
otin)(2-\alpha_{ 1})+[(2-\alpha_{0})\kappa+1]\lambda_{1}}{2[(2-\alpha_{0})(2-\alpha_{1})\kappa +(1-\alpha_{1})]T_{\lambda_{1}}}\pm\frac{\sqrt{\left\{(\kappa M^{2}T_{ \lambda_{1}}+\alpha_{0}\
otin)(2-\alpha_{1})-[(2-\alpha_{0})\kappa+1]\alpha_{1 }\right\}^{2}+4\alpha_{1}\left(\kappa M^{2}T_{\lambda_{1}}+\alpha_{0}\
otin \right)}}{2[(2-\alpha_{0})(2-\alpha_{1})\kappa+(1-\alpha_{1})]T_{\lambda_{1}}} \tag{29}\]  
\[\left\{ \begin{array}{c}A_{1}=\frac{\kappa M^{2}}{(x_{1}-x_{2})T_{\lambda_{1}}} \frac{(2-\alpha_{1})x_{1}T_{\lambda_{1}}+\alpha_{1}}{(2-\alpha_{0})(2-\alpha_ {1})\kappa+(1-\alpha_{1})}\\ A_{2}=\frac{\kappa M^{2}}{(x_{2}-\lambda_{1})T_{\lambda_{1}}}\frac{(2-\alpha_{1 })x_{2}T_{\lambda_{1}}+\alpha_{1}}{(2-\alpha_{0})(2-\alpha_{1})\kappa+(1- \alpha_{1})}\end{array} \right. \tag{30}\]  
\[\left\{ \begin{array}{c}T_{\
u}=c_{\
u}t/H^{2}\\ T_{\tau}=c_{\
u}\tau/H^{2}\\ T_{\lambda_{1}}=c_{\
u}\lambda_{1}/H^{2}\end{array} \right. \tag{31}\]  
\[\
otin=\eta_{1}/\eta_{0} \tag{32}\]Substituting Equation (28) into Equation (23) produces the general solution of excess pore-water pressure in convolution form. The result is  
\[u(z,t)=q_{0}\sum_{m=1}^{\infty}\left[\phi(T_{\rm v})-\int_{0}^{T_{\rm v}}\phi(T_{ \rm r})\sum_{j=1}^{2}A_{j}e^{x_{j}(T_{\rm v}-T_{\rm r})}dT_{\rm r}\right]C_{m} \sin\left(\frac{M}{H}z\right) \tag{33}\]  
where \(\phi(t)\) is a dimensionless loading function given by  
\[0\leq\phi(t)=\frac{q(t)}{q_{0}}\leq 1 \tag{34}\]  
and \(q_{0}=\max[q(t)]\) is the maximum of time-dependent loading \(q(t)\).