## 1 Introduction  
The long-term behaviors of rock masses around underground structures have received increasing attention from both academia and industry since these behaviors have significant influence on the structures' stability [1; 2; 3]. Creep behaviours refer to the continuous re-distribution of stress and displacement that occurs when the shear stress exceeds a threshold. For underground structures in weak rock masses or rock masses under high in situ stresses, creep behaviors are particularly common. Moreover, creep behaviors must be investigated for some underground structures that need to be operated for a long period of time, i.e., thousands to tens of thousands years for repositories for radioactive waste disposal [4].  
Significant advances have been achieved towards understanding the creep behaviors of various rock types, such as red sandstone [5], rock salt [6], cataclastic rocks [7], granite [8; 9], green tuff [10], and limestone [11]. At present, the investigation of the creep behaviors of rock can be categorized into experimental studies and numerical modelling. Experimental studies are often performed using creep tests on various rock types and under various conditions, such as different confining pressures and temperatures [12]. Creep tests can be conducted in either single-stage or multi-stage manners [13; 14; 15]. Since single-stage loading is time-consuming, multi-stage loading has been increasingly employed to obtain the complete creep characteristics of rock, including deceleration creep, constant velocity creep and accelerated creep [16]. Despite many investigations on the mechanical properties of argillaceous sandstone (AS), the creep behaviors of AS under different confining pressures have not been well documented.  
Numerical modelling is another way to analyze the creep behaviors of rocks. For example, the interaction between rock masses and a backfill body considering the creep behaviors of rock masses be investigated [3]. Possible explanations were provided for in situ anomalies after the introduction of the creep constitutive model. One important issue in numerical modelling is the selection of the constitutive model and the corresponding parameters [17; 18; 19]. Various constitutive models have been established to represent the creep characteristics of different materials.  
The traditional creep constitutive model is established without proper consideration of the tertiary creep stage [20; 21]. Towards this end, nonlinear viscoplastic creep models have been proposed [22; 23]. For example, [PERSON] proposed a nonlinear elastic viscoplastic damage model based on the modified Mohr-Coulomb criterion and investigated the creep and permeability of clay rock during the construction of a high-level waste repository [24]. [PERSON] investigated the instantaneous deformation characteristics of rock creep by establishing an elastic-plastic model describing the instantaneous behavior [25]. However, the constitutive creep relations for AS are still in the infancy. In other words, more efforts should be directed to the constitutive modelling of the creep behavior of AS.  
Therefore, an intensive study on the creep behaviors of AS under various confining pressures was conducted. Multi-stage creep tests were carried out first on AS specimens with different confining pressures. The creep characteristics of AS were analyzed based on the experimental results. Then, a novel constitutive relation was proposed to describe the creep behaviors of AS, which was embedded into FLAC3D for implementation. Finally, the experimental and numerical results were compared to validate the feasibility of the proposed creep model. A systematic study about experiment, theory model, and numerical calculation of AS' creep characteristics has been presented in this paper. Furthermore, major rock creep models only stay at the stage of theoretical explanation, this paper embeds the constructed creep model into numerical calculation software to verify its correctness.