## 2 Visco-Elastic-Plastic Constitutive Model for Squeezing Rock
### Strain Integration Scheme of the Rheological Model  
The finite element method (FEM) is a powerful numerical tool for simulating the behavior of subsurface structures under different geological and environmental conditions [11; 36]. This study uses the ABAQUS software for numerical simulation, and the  
Figure 3: One-dimensional Komamura-Huang-Bingham model. proposed constitutive model is implemented within the UMAT subroutine. As a result, the visco-elastic-plastic rheological model can simulate the large deformations in the soft rock mass that take place during tunnel excavation. An explicit strain integration scheme was applied to compute the strain increment. For the proposed visco-elastic-plastic constitutive model, the total strain increment (\(\{\Delta\varepsilon\}\)) can be decomposed into the instantaneous elastic strain increment (\(\{\Delta\varepsilon_{\varepsilon}\}\)), the viscoelastic strain increment (\(\{\Delta\varepsilon_{\varepsilon\varepsilon}\}\)), and the viscoplastic strain increment (\(\{\Delta\varepsilon_{\varepsilon\varepsilon}\}\)):  
\[\{\Delta\varepsilon\}=\{\Delta\varepsilon_{\varepsilon}\}+\{\Delta\varepsilon _{\varepsilon\varepsilon}\}+\{\Delta\varepsilon_{\varepsilon\varepsilon}\}. \tag{4}\]  
The elastic strain increment can be derived from [PERSON]'s elasticity theory:  
\[\{\Delta\varepsilon_{\varepsilon}\}=[D]^{-1}\{\Delta\sigma\}, \tag{5}\]  
where \([D]\) represents the elastic stiffness matrix of the rock mass, and \(\{\Delta\sigma\}\) is the stress increment. For the viscoelastic rheology, two Kelvin bodies are adopted in the proposed nonlinear rheological model, as shown in Figure 3. The superposition of the respective strains of the two Kelvin bodies provides the total viscoelastic strain. According to the constitutive law of the Kelvin model, the viscoelastic strain, \(\{\varepsilon_{\varepsilon\varepsilon}\}_{t+\Delta t}\), at a specific time-step can be obtained as  
\[\{\varepsilon_{\varepsilon\varepsilon}\}_{t+\Delta t}=\sum_{i=1}^{2}\biggl{(} \{\varepsilon_{\varepsilon\varepsilon}\}_{t}\cdot e^{-\frac{E_{i}}{\eta_{i} }\Delta t}+\frac{E_{0}}{E_{i}}[D]^{-1}\cdot\{\sigma\}\cdot\biggl{(}1-e^{- \frac{E_{i}}{\eta_{i}}\Delta t}\biggr{)}\biggr{)} \tag{6}\]  
where \(\{\varepsilon_{\varepsilon\varepsilon,i}\}_{t}\) and \(\{\varepsilon_{\varepsilon\varepsilon,i}\}_{t+\Delta t}\) are the viscoelastic strains of the \(i\)-th Kelvin model at two consecutive time steps \(t\) and \(t+\Delta t\), respectively. \(\Delta t\) is the time step increment, and \(\{\sigma\}\) is the stress state exerted on the rock mass. As a result, the increment of the viscoelastic strain at a specific time step is  
\[\{\Delta\varepsilon_{\varepsilon\varepsilon}\}=\{\varepsilon_{\varepsilon \varepsilon}\}_{t+\Delta t}-\{\varepsilon_{\varepsilon\varepsilon}\}_{t} \tag{7}\]  
After the stress state (\(\{\sigma\}\)) of the soft rock mass reaches the plastic yield stress state, the rock mass will experience plastic deformation. The plastic strain develops according to the yield surface of the adopted Drucker-Prager yielding criterion and the flow rule of the nonlinear Bingham body. The yield surface of the DP yielding criterion is defined as  
\[F(\sigma,c)=\sqrt{J_{2}}+\eta I_{1}-\xi c \tag{8}\]  
where \(I_{1}\) is the first stress invariant, \(J_{2}\) is the second deviatory stress invariant, and \(c\) is the cohesion of the rock mass. \(\eta\) and \(\zeta\) are material variables related to the internal frictional angle of the rock mass:  
\[\eta=\frac{\sqrt{3}\sin\varphi}{3\sqrt{\bigl{(}3+\sin^{2}\varphi\bigr{)}}} \tag{9}\]  
\[\zeta=\frac{\sqrt{3}\cdot\cos\varphi}{\sqrt{3}\sqrt{\bigl{(}3+\sin^{2}\varphi \bigr{)}}} \tag{10}\]  
The flow rule specifies the relationship between the plastic strain increment and the stress increment.