<p>Constitutive models currently implemented in numerical simulation tools had limited applicability and cannot accurately characterize the complex creep behavior of magnetite. Here, we derive a three-dimensional central-difference formulation from the one-dimensional constitutive equation of the nonstationary viscoelastic-plastic BuCY creep model, based on rock rheology theory. We further enhance the numerical solution framework to perform iterative stress-strain updates within each time step. During the accelerating creep stage of the nonstationary viscoplastic CY component, the plastic behaviour of element Y is governed by the Mohr-Coulomb instability criterion. The BuCY creep model able to reproduce the observed creep characteristics of the magnetite-was integrated into the finite-difference code via secondary development, enabling model compilation and solution algorithm and inverse verification. Numerical simulation and inversion results show close agreement, confirming the success of the implementation. Comparative analysis with conventional creep models further confirms the superior applicability and accuracy of the BuCY model. This research significantly improves the fidelity of numerical simulations and is particularly well suitable for practical engineering problems, such as the analysis of deeply buried mining stopes.</p>