## 2 Methodology
### The Numerical Calculation Platform and the Implementation of the Creep Model  
This study employs the secondary development of the Nishihara creep model in Abaqus software, version 6.14. As Abaqus's native creep models inadequately capture the visco-plastic creep characteristics of fractured rock masses under varying depth-dependent environmental conditions, this study necessitates secondary development by consequently integrating the Nishihara model.  
In Abaqus software, the finite element governing equations for nonlinear problems can be simplified as  
\[\sigma(\varepsilon)=0 \tag{2}\]  
where \(\sigma\) is the stress vector, and \(\varepsilon\) is the displacement/strain vector. Assuming Abaqus obtains an approximate solution \(\varepsilon_{\mathrm{i}}\) from the input stress vector \(\sigma_{\mathrm{i}}\) during the \(i\)-th iteration, with a residual error \(c_{i+1}\) relative to the exact solution, this relationship is expressed as  
\[\begin{split}\sigma(\varepsilon_{i})\
eq 0\\ \sigma(\varepsilon_{i}+c_{i+1})=0\end{split} \tag{3}\]  
Figure 1: A schematic diagram of the Nishihara model.  
Expanding Equation (3) via Taylor series and neglecting higher-order terms yields  
\[Jc_{i+1}=-\sigma_{i} \tag{4}\]  
where \(J=\partial\Delta\sigma/\partial\Delta\varepsilon\) is the Jacobian stiffness matrix, and \(\varepsilon_{i+1}=\varepsilon_{i}+c_{i+1}\) is the strain increment for the next iteration.  
Thus, the core requirement for the secondary development of user-defined constitutive models of the rock mechanics lies in providing an explicit formulation of the Jacobian stiffness matrix that describes the material's constitutive relationship. When the desired material model is unavailable in Abaqus's native library, the user material subroutine (UMAT) interface allows for the integration of custom constitutive equations into the software for numerical implementation.  
According to existing successful cases upon the numerical implementation of the Nishihara creep model in the UMAT subroutine in Abaqus in previous studies [19; 20; 21; 22], the secondary development workflow adopted in this study is illustrated in Figure 2. Using the Fortran language in the Microsoft Visual Studio 2010 compiler, we programmed the Nishihara creep model subroutine.  
Step 1: Formulating the initial Jacobian matrix DDSDDE \(\sigma_{\mathbf{t}}\) based on the incremental form of the constitutive equations;  
Step 2: Calculating the stress increment \(\Delta\sigma_{\mathbf{t+1}}\) from the strain increment \(\Delta\varepsilon_{\mathbf{t}}\).  
Step 3 Updating the stress components \(\sigma_{\mathbf{t+1}}\) at the \(\mathbf{t+1}\) timestep and assigning them to the new Jacobian matrix DDSDDE.  
Note: The strain matrix STRAN and the strain increment matrix DSTRAN in Abaqus employ the engineering shear strain, which is twice the magnitude of the mathematical shear strain.  
Using this workflow, the Nishihara creep model was implemented into Abaqus. Subsequent numerical studies focused on the long-term mechanical response of the surrounding rock in critical sections of the Jinping diversion tunnels at the engineering scale. For the seepage equations, Abaqus's default porous medium flow was adopted within an equivalent continuous-medium framework, with the effects of seepage primarily incorporated via boundary conditions into subsequent simulations. The UMAT implementation proved stable throughout the long-term creep analyses, exhibiting no significant convergence difficulties, timestep limitations, or element distortion.