## 3 Solutions to Governing Equations  
In order to simplify the solution process, the following dimensionless parameters are given: \(U=\frac{u}{\sigma_{0}^{\prime}}\), \(Z=\frac{z}{H}\), \(T_{\rm v}=\frac{c_{\rm v0}t}{H^{2}}\), \(c=c_{\rm c}/c_{\rm k}\), \(b=(q_{0}+\sigma_{0}^{\prime})/\sigma_{0}^{\prime}\), \(I=(i\gamma_{\rm w}H)/\sigma_{0}^{\prime}\)and \(I_{1}=(i_{1}\gamma_{\rm w}H)/\sigma_{0}^{\prime}\). Substituting these dimensionless parameters into the governing equations and continuous conditions, the corresponding equations can be given as:  
\[\left\{ \begin{array}{l}(b-U)\frac{\partial}{\partial Z}\left[\frac{(b-U)^{-c}}{ m_{1}^{m_{1}-1}}\Big{(}\frac{\partial U}{\partial Z}\Big{)}^{m-1}\frac{ \partial U}{\partial Z}\right]=\frac{\partial U}{\partial T_{\rm v}}\,I\leq I_{1}\\ (b-U)\frac{\partial}{\partial Z}\left[(b-U)^{-c}\frac{\partial U}{\partial Z} \right]=\frac{\partial U}{\partial T_{\rm v}},I\geq I_{1}\end{array} \right., \tag{15}\]  
\[U(0,T_{\rm v})=(b-1){\rm e}^{-\pi T_{\rm v}}, \tag{16}\]  
\[\left.\frac{\partial U}{\partial Z}\right|_{Z=1}=0, \tag{17}\]  
\[U(Z,0)=b-1. \tag{18}\]  
The soil layer is discretized in space and time, and \(i=0,1,2,\ldots,n\) and \(j=0,1,2,\ldots,N\) represent the space and time notes, respectively. The \(0\)th space note means the top surface of the soil layer, and the \(0\)th time note mean the initial moment. The time and space steps are denoted as \(\tau\) and \(h\), respectively. Setting \(\lambda=\frac{\tau}{h^{2}}\) and using the implicit scheme of a quasilinear diffusion equation, the following difference equation can be obtained:  
\[U_{i}^{j+1}=U_{i}^{j}+\lambda(b-U_{i}^{j})\Big{[}\alpha_{i+1/2}^{j+1}\Big{(}U_ {i+1}^{j+1}-U_{i}^{j+1}\Big{)},-\alpha_{i-1/2}^{j+1}\Big{(}U_{i}^{j+1}-U_{i-1} ^{j+1}\Big{)}\Big{]}, \tag{19}\]  
where  
\[\alpha_{i+1/2}^{j+1}=\left\{\begin{array}{l}\frac{1}{m_{1}^{m_{1}-1}}  
\[\mathbf{V}=\left[ \begin{array}{c}U_{1}\\ U_{2}\\ U_{3}\\ \vdots\\ U_{n-1}\\ U_{n}\end{array} \right], \tag{29}\]  
\[\mathbf{B}=\left[ \begin{array}{c}B_{1}\\ B_{2}\\ B_{3}\\ \vdots\\ B_{n-1}\\ B_{n}\end{array} \right], \tag{30}\]  
\[A_{ii}=1+\lambda(b-U_{i}^{j})\alpha_{i+1/2}^{j+1}+\lambda(b-U_{i}^{j})\alpha_ {i-1/2^{\prime}}^{j+1} \tag{31}\]  
\[A_{i(i-1)}=-\lambda(b-U_{i}^{j})\alpha_{i-1/2^{\prime}}^{j+1} \tag{32}\]  
\[A_{i(i+1)}=-\lambda(b-U_{i}^{j})\alpha_{i+1/2^{\prime}}^{j+1} \tag{33}\]  
\[A_{mn}=1+\lambda(b-U_{n}^{j})\alpha_{n+1/2}^{j+1}+\lambda(b-U_{n}^{j})\alpha_ {n-1/2^{\prime}}^{j+1} \tag{34}\]  
\[A_{n(n-1)}=-\lambda(b-U_{n}^{j})\alpha_{n-1/2}^{j+1}-\lambda(b-U_{n}^{j}) \alpha_{n+1/2^{\prime}}^{j+1} \tag{35}\]  
\[B_{i}=U_{i}, \tag{36}\]  
\[B_{1}=U_{1}^{j-1}+\lambda(b-U_{1}^{j})\alpha_{1-1/2}^{j+1}U_{0}^{j}, \tag{37}\]  
\[\alpha_{n+1/2}^{j-1}=\left(b-\frac{U_{n-1}^{j-1}+U_{n}^{j-1}}{2}\right)^{-c}, \tag{38}\]  
\[\alpha_{n-1/2}^{j-1}=\left(b-\frac{U_{n}^{j-1}+U_{n-1}^{j-1}}{2}\right)^{-c}. \tag{39}\]  
It can be seen that the unknowns of Equation (23) are linear, and the tridiagonal matrix can be solved by the catch-up method with Matlab. In order to speed up the solving process of the difference equation, \(h\) is set as 0.01. When \(T_{\rm v}<0.2\), \(\tau\) is set as 0.000001; when \(T_{\rm v}>0.2\), \(\tau\) is set as 0.0001. Therefore, when \(T_{\rm v}<0.2\), \(\lambda\) is set as 0.01; when \(T_{\rm v}>0.2\), \(\lambda\) is set as 1. With these values, the solution of EPWP can be obtained.