## 1 Introduction  
The soft rocks are widely distributed in situ and closely related to various geotechnical facilities, including slope, shallow tunnel, and rock-sacketed pile [1; 2; 3; 4]. Some soft sedimentary rocks contain a few discontinuous structures and are considered continuous media. The mechanical behaviour of soft rocks such as strain-softening, dilatancy, and confining pressure dependency were comprehensively studied within the framework of critical state soil mechanics, where the soft rock was regarded as a heavily overconsolidated soil [5; 6; 7]. It has also been found that temperature plays a significant role in the mechanical behaviour of soft rocks [8; 9; 10; 11]. Specifically, an increase in the temperature will reduce the peak strength of soft rock and lead to a ductile failure pattern with a smaller tangent modulus.  
To investigate the thermo-mechanical behaviour of soft rocks, [PERSON] et al. [10] developed the concept of temperature-induced equivalent stress to account for the thermoelastic volumetric deformation. Following this concept, the effect of temperature \(T\) on the precon-solidated pressure \(p_{cT}\) can be obtained through the isotropic compression test. Increasing the temperature will shrink the normal yield surface along with a so-called temperature-dependent collapse curve, resulting in a reduction of the strength of soft rock. To obtain a smooth transition from the elastic and plastic responses of soft rock, the subloading surface concept [12; 13] is also suggested in developing an elastoplastic model of soft rock. The subloading surface that passes through the current stress point is similar to the normal yield surface, which can be used to capture the irreversible plastic deformation when the current stress point is inside the normal yield surface. With temperature-induced equivalent stress and subloading surface concepts, [PERSON] et al. [10] developed a simple thermo-elastoplastic model for soft rock. [PERSON] et al. [5] extended the model by regarding the soft rock as a structured overconsolidated soil to study the processes of temperature-induced overconsolidated ratio dissipation and structure degradation with great success in modelling the triaxial test data of Ohya stone.  
The associated flow rule was adopted in those models by assuming the strain vectors are normal to the yield surface, which may not be consistent with the experimental observations for some geomaterials [14]. [PERSON] and [PERSON] [15] suggested that the non-associated flow rule is imperative to capture the legitimate behaviour of geomaterials; otherwise, the predicted shear strain will be overestimated. The recently developed fractional-order flow rule can be a promising approach to account for this behaviour without increasing the complexity of the constitutive equations [16; 17; 18]. In general, a so-called fractional stress operator is directly applied to the yield function to determine the plastic strain tensor without introducing an extra plastic potential.  
This paper aims to develop a fractional subloading surface model for soft rock within the framework of critical state soil mechanics. The proposed model contains two yield surfaces representing the normal consolidated and overconsolidated states. The change in the size of those yield surfaces is captured by a new hardening rule, which can consider the effect of temperature on the strain-softening and dilatancy features of soft rocks. A fractional-order flow rule is also included to describe the nonorthogonality between the plastic flow direction and the subloading yield surface to improve model performance. The model is verified by the drained triaxial test data of Ohya rock with different confining pressures and temperatures.