## 2 NGNRM with Non-Stationary Parameter Creep  
The accelerated creep of soft rock plays an important role in the large deformation disasters of tunnels. From the introduction, it is known that the previous studies seldom reflect the accelerated creep in the whole process of soft rock large deformation. In this paper, we used the NGNRM (see Figure 1) representing the accelerated creep characteristics of soft rocks, to investigate the large deformation mechanism of soft rocks. The NGNRM consists of a Hooke body, a Kelvin body, and a viscoplastic body with the viscosity coefficient as a function of stress. In Figure 1, \(E_{0}\) is the elastic modulus of a Hooke body; \(E_{1}\) and \(\eta_{1}\) are the elastic modulus of a Hooke body and the viscosity coefficient of a Newton body in a Kelvin body, respectively; \(\sigma_{\rm f}\) and \(\eta_{2}\) are the yield stress of a St. Venant body and the viscosity coefficient of a dashpot in the non-stationary viscoplastic body, respectively, and \(\eta_{2}\) is as a function of stress \(\sigma\) and time \(t\). When \(\sigma\geq\sigma_{\rm f}\), for the nonlinear Bingham body, the expression of the viscosity coefficient \(\eta_{2}\) is denoted by Equation (1) [28]. \[\left\{ \begin{array}{l}\varepsilon=\frac{\sigma}{E_{0}}+\frac{A(\sigma- \sigma_{\rm f})t}{E_{0}\big{[}T+t\big{(}1-\frac{\sigma-\sigma_{\rm f}}{\sigma_ {\rm s}}A\big{)}\big{]}}\\ \\ \eta_{2}(\sigma)=\frac{E_{0}}{AT}\Big{[}T+t\Big{(}1-\frac{\sigma-\sigma_{\rm f }}{\sigma_{\rm s}}\cdot A\Big{)}\Big{]}^{2}\end{array} \right. \tag{1}\]  
where \(\sigma\) is the current stress; \(\sigma_{\rm s}\) is the long-term strength of a Bingham body; \(A\) is a dimensionless parameter; \(T\) is the time for the surrounding rock to yield; \(t\) is load time. It is noted from Equation (1) that when \(\sigma=\sigma_{\rm f}+\sigma_{\rm s}/A\), nonlinear creep degenerates into linear creep, \(\eta=E_{0}T/A\), and strain increases at the constant rate \(\dot{\varepsilon}=\sigma_{\rm s}/E_{0}T\), that is, steady creep; when \(\sigma>\sigma_{\rm f}+\sigma_{\rm s}/A\), \(\eta_{2}\) gradually decreases with \(t\), reflecting the accelerated creep process of soft rock in a high-stress state. Therefore, when \(\sigma\geq\sigma_{\rm f}\), the one-dimensional creep equation of the nonlinear viscoelastic-plastic model is denoted by Equation (2). \[\varepsilon(t)=\sigma\bigg{\{}\frac{1}{E_{0}}+\frac{1}{E_{1}}\bigg{[}1-\exp \bigg{(}-\frac{E_{1}}{\eta_{1}}t\bigg{)}\bigg{]}\bigg{\}}+\frac{A(\sigma- \sigma_{\rm f})t}{E_{0}\Big{[}T+t\Big{(}1-\frac{\sigma-\sigma_{\rm f}}{\sigma_ {\rm s}}A\Big{)}\Big{]} \tag{2}\]  
In this paper, the initial strain method is used to calculate the viscoelastic-plastic constitutive equations of the NGNRM, and its ABAQUS user material subprogram is developed by FORTRAN, which is embedded in the ABAQUS. The flowchart of the user-defined material mechanical behavior (UMAT) subroutine development for the NGNRM is shown in Figure 2.  
UMAT is a FORTRAN program interface provided by ABAQUS to users to define their own material properties, enabling users to use material models that are not defined in the ABAQUS material library. The UMAT realizes the data exchange with ABAQUS through the interface with the ABAQUS main solver. Figure 1: Rheological model with variable viscoplastic viscosity coefficient. Drucker-Prager yield criterion is adopted in this paper, denoted by Equation (3). \[\left\{ \begin{array}{l}F(\mathbf{\sigma},c)=\sqrt{J_{2}}+\eta I_{1}-\xi c\\ \eta=\frac{\sqrt{3}\sin\varphi}{3\sqrt{(3+\sin^{2}\varphi)}}\\ \bar{\xi}=\frac{\sqrt{3}\cos\varphi}{\sqrt{3}\sqrt{(3+\sin^{2}\varphi)}}\end{array} \right. \tag{3}\]  
where \(c\) and \(\varphi\) are the internal cohesion and internal friction angle, respectively; \(I_{1}\) is the first invariant of stress; \(J_{2}\) is the second invariant of deviatoric stress. The UMAT subroutine of the NGNRM are implemented by following steps:  
_Step 1_ Call the UMAT subroutine.