## 2 Governing Equations
### Temperature-Dependent Subloading Surface  
This study regards soft rock as a heavily overconsolidated clay and investigates the effect of temperature on its strain-softening and dilatancy features within the framework of the critical state mechanics. The concept of temperature-induced equivalent stress developed by [PERSON] et al. [10] is adopted here and can be calculated by  
\[p_{T}=p\text{exp}\left[\frac{3\alpha_{t}(T-T_{0})(1+e_{0})}{\kappa}\right] \tag{7}\]  
where \(p_{T}\) and \(p\) are the current and reference mean effective stress at temperatures \(T\) and \(T_{0}\), respectively. The value of \(T_{0}\) can be arbitrarily selected and is chosen to be \(15\,^{\circ}\text{C}\). \(\alpha_{t}\) is a negative linear thermal expansion coefficient (e.g., compaction is assumed to be positive). \(e\) and \(e_{0}\) are the current and reference void ratios related to the deformation of the soft rock, respectively. \(\kappa\) is the slope of the swelling line in the \(e-p\) plane. Figure 2 shows the Matsuoka-Nakai criteria in the principal stress space with different temperatures using the approach developed by [PERSON] and [PERSON] [20]. It can be found that increasing the temperature will decrease the size of the criteria, leading to a smaller elastic zone. Furthermore, rock samples at higher temperatures will be yielded earlier for a single yield surface elastoplastic model. The phenomenon of temperature-induced strength reduction can be observed at the peak strength state. However, it is difficult to capture  
Figure 1: The schematic of the fractional-order flow rule. the accurate residual strength of rocks because the samples usually fail before entering the critical state. In that case, the effect of temperature on the residual strength of soft rocks remains unclear in the literature and is ignored in the current study for simplicity. To describe the effect of temperature on the mechanical behaviour of soft rocks, the concept of subleading surface [12; 13] is also required, which can capture the plastic deformation at the early loading stage due to the change in temperature. As illustrated in Figure 3, in the principal stress space, the subleading surface always passes through the current stress point. The location corresponding to the normal yield surface is determined by a similarity ratio \(R\) which is defined as the ratio between the current mean effective stress \(p\) and the reference mean effective stress \(\overline{p}\). Similar to the modified Cam-clay (MCC) model [21], the normal yield surface is represented as  
\[\overline{f}=\overline{q}^{2}+M_{cs}^{2}\overline{p}(\overline{p}-\overline{p} _{\mathcal{CT}})=0 \tag{8}\]  
where \(M_{cs}\) is the slope of the critical state line in the \(p-q\) plane. \(\overline{p}_{\mathcal{CT}}\) is the temperature-dependent preconsolidated pressure that controls the size of the normal yield surface and is calculated by Equation (7). Figure 3: The schematic of the temperature-dependent subleading surface concept. Figure 2: Matsuoka–Nakai criteria in the principal stress space with different temperatures. The similarity ratio \(R\) has a similar physical meaning to the overconsolidated ratio which can be further illustrated as  
\[R=\frac{p}{\overline{p}_{0i}}\frac{M_{cs}^{2}}{M_{cs}^{2}+\eta^{2}}\cdot\exp \left[-\frac{3\alpha_{t}(T-T_{0})(1+e_{0})}{\kappa}\right] \tag{9}\]  
where \(\eta\) is the stress ratio defined as \(\eta=q/p\). Similarly, the subloading surface is defined as  
\[f=q^{2}+M_{cs}^{2}p(p-p_{cT})=0 \tag{10}\]  
where \(p_{cT}\) is the intersection of the subloading surface with the \(p\)-axis. The change in the size of the subloading surface is controlled by \(p_{cT}\), which is chosen as the hardening parameter with the following evolution law:  
\[dp_{cT}=\frac{1+e_{0}}{\Lambda-\kappa}p_{cT}\bigg{(}1-\frac{\ln R}{RS}\bigg{)} de_{v}^{p} \tag{11}\]  
where \(S\) is a so-called softening factor that ensures the current stress point passes through the critical state line to produce dilative deformation.