## 3 Computational Framework
### Work Calculation of the Failure Mechanism  
The above mechanism is utilized to compute the total work in a spatially discretized manner, as the stability issue of the tunnel can be translated into the equilibrium relationship between the internal dissipation and the external work of the entire failed rigid block based on the principle of virtual work. The calculation of work associated with the failure mechanism can be categorized into internal energy dissipation rates governed by effective cohesion and apparent cohesion, as well as external work rates influenced by soil weight and a face support force. To accurately compute the cumulative work rates, the point-to-point generation mechanism facilitates this process. Geometric parameters in space, along with their work expressions for small surface planes and prisms, are discretely formulated. The total work rate is the aggregate of the work performed by the corresponding small elemental units. As shown in Figure 4, the right big cube contains two triangle planes and when two triangular planes are projected onto a symmetrical surface, they can form two cubes. The triangle planes \(S_{\text{i,j}}\) and \(S_{\text{i,j}}{}^{\prime}\) are used to calculate the work rate of effective cohesion and apparent cohesion, while the cube is used to compute the work rate of soil weight. The rectangular plane on the tunnel face is used to compute the work rate of the face support. As illustrated in Figure 4, the large cube on the right contains two triangular planes which each form a cube when these triangular planes are projected onto a symmetrical surface. The triangular planes \(S_{\text{i,j}}\) and \(S_{\text{i,j}}{}^{\prime}\) are employed in determining the work rates for effective cohesion \(W_{c^{\prime}}\) and apparent cohesion \(W_{capp}\), whereas the cubes \(V_{\text{i,j}}\) and \(V_{\text{i,j}}{}^{\prime}\) are utilized for calculating the work rates associated with soil weight \(W_{\gamma}\). The rectangular plane \(\Sigma_{j}\) situated on the tunnel is used to compute the work rate related to face support \(W_{c_{T}}\). According to the principle of virtual work, the balance relationship of these work rates is:  
\[W_{\gamma}+W_{\sigma_{T}}=W_{c^{\prime}}+W_{c_{app}} \tag{26}\]  
in which the formulae of the work terms are as follows:  
\[W_{c^{\prime}}+W_{c_{app}}=\omega\underset{i}{\sum}\underset{j}{\sum}\,(c^{ \prime}+c_{app})(\cos\varphi_{i,j}R_{i,j}S_{i,j}+\cos\varphi^{\prime}{}_{i,j}R ^{\prime}{}_{i,j}S^{\prime}{}_{i,j}) \tag{27}\]  
Figure 4: The tiny surface plane and cube for work calculation in discretization method. where the \(\varphi_{i,j}\) (or \(\varphi^{\prime}_{i,j}\)) can be obtained from the coordinates of the centroid of the small triangle \(S_{i,j}\) (or \(S^{\prime}_{i,j}\)) and the center of rotation. The \(R_{i,j}\) (or \(R^{\prime}_{i,j}\)) represents the distance between these two points. The expression of \(c^{\prime}\) and \(c_{app}\) also vary in space:  
\[c^{\prime}+c_{app}=\left\{ \begin{array}{l}c^{\prime}1+\frac{\tan\varphi^{ \prime}_{1}h_{ml}(y)\gamma_{\pi\eta}}{\left\{1+\left[a_{1}h_{ml}(y)\gamma_{ \pi\eta}\right]^{1/\eta}\right\}^{1/\eta}},\ X_{N}x_{i,j}+Y_{N}(y_{i,j}-dt)+Z_{ N}z_{i,j}\leq 0\\ \\ c^{\prime}2+\frac{\tan\varphi^{\prime}_{2}h_{2m}(y)\gamma_{\pi\eta}}{\left\{1 +\left[a_{2}h_{m2}(y)\gamma_{\pi\eta}\right]^{1/\eta}\right\}^{1/\eta}},\ X_{N}x_{i,j}+Y_{N}(y_{i,j}-dt)+Z_{ N}z_{i,j}>0\end{array} \right. \tag{28}\]  
where \((x_{i,j},y_{i,j},z_{i,j})\) is the Cartesian coordinate of the centroid of the facet \(S_{i,j}\) (or \(S^{\prime}_{i,j}\)). It should also be noted that the values of \(l_{1}\) in \(h_{m2}(y)\) vary with different horizontal positions:  
\[l_{1}=d_{t}-\frac{X_{N}x_{i,j}+Z_{N}z_{i,j}}{Y_{N}}-y_{0} \tag{29}\]  
where \(y_{0}\) is the coordinate of the groundwater level's elevation.