## 2 Theoretical Background
### Energy Decomposition and Driving Force of Phase-Field Evolution  
While the classical phase-field formulation (Equation (14)) does not distinguish between material failure modes, geological materials exhibit complex fracture behaviors under multi-physics coupling, including tensile, compressive, tensile-shear, and compressive-shear failure. To address this limitation, researchers have developed enhanced phase-field frameworks. Notably, [PERSON] et al. [41; 42] pioneered an energy decomposition approach that partitions the elastic strain energy density into volumetric \(\psi_{\epsilon_{\epsilon}}^{sph}\) and deviatoric \(\psi_{\epsilon_{\epsilon}}^{dev}\) components:  
\[\epsilon_{\epsilon}=\epsilon_{sph}+\epsilon_{dev},\ \epsilon_{sph}=1/3 \text{tr}(\epsilon_{\epsilon})\mathbf{I},\ \epsilon_{dev}=\epsilon_{\epsilon}-\epsilon_{sph} \tag{18}\]  
\[\psi_{\epsilon_{\epsilon}}^{sph}=\frac{K}{2}\text{tr}^{2}(\epsilon_{\epsilon} ),\ \psi_{\epsilon_{\epsilon}}^{dev}=\mu\text{tr}\Big{(}\epsilon_{dev}^{2} \Big{)} \tag{19}\]  
where \(K\) is the bulk modulus, \(K=(3\lambda+2\mu)/3\).  
Complementarily, [PERSON] et al. [21] proposed a spectral decomposition approach that partitions the elastic strain energy density into the tensile and compressive part \(\psi_{\epsilon_{\epsilon}}^{\pm}(\epsilon_{\epsilon})\) according to the spectral decomposition method.  
\[\epsilon_{\epsilon}^{\pm}=\sum\limits_{i=1}^{3}<\epsilon_{i}>_{\pm}\mathbf{n}_{i} \otimes\mathbf{n}_{i} \tag{20}\]\[\psi^{\pm}_{\epsilon_{e}}(\epsilon_{e})=\frac{\lambda}{2}<\text{tr}(\epsilon_{e})>^{ 2}_{\pm}+\mu\text{tr}(<\epsilon_{e}>^{2}_{\pm}) \tag{21}\]  
where \(\epsilon_{i}\) denotes the principal strain, and the operator \(<x>_{\pm}\) is defined as \(<x>_{\pm}=(x+|x|)/2\).  
Thermal fracturing in rocks under elevated temperatures predominantly manifests as tensile cracking. Although the classical PFM cannot capture complex mixed-mode failures, it demonstrates proven effectiveness in simulating tensile-dominated fracture. Consequently, this study adopts the classic PFM, specifically restricting material stiffness degradation to the tensile strain energy contribution.  
\[\psi_{\epsilon_{e}}(\epsilon_{e})=\Big{[}(1-d)^{2}+\delta_{0}\Big{]}\psi^{+}_{ \epsilon_{e}}(\epsilon_{e})+\psi^{-}_{\epsilon_{e}}(\epsilon_{e}) \tag{22}\]  
Therefore, Equation (14) can be rewritten as:  
\[\frac{G_{c}}{l_{0}}\Big{(}d-l_{0}^{2}\
abla^{2}d\Big{)}-2(1-d)H_{t}=0 \tag{23}\]  
where \(H_{t}=\max\limits_{s\in[0,d]}\left[\psi^{+}_{\epsilon_{e}}(\epsilon_{e})-\psi _{c}\right]_{+}\), \(\psi_{c}=\frac{f_{t}^{2}}{2E}\) represents the stored energy threshold for damage evolution, and \(f_{t}\) is the uniaxial tensile strength.  
The historical strain field is dynamically updated through COMSOL Multiphysics' state variables, eliminating the need for auxiliary ODE/DAE interfaces traditionally required to track fracture history.  
Additionally, the stress tensor is rewritten as:  
\[\mathbf{\sigma}(\epsilon_{e},d)=\frac{\partial\psi_{\epsilon_{e}}(\epsilon_{e},d) }{\partial\epsilon_{e}}=g(d)\frac{\partial\psi^{+}_{\epsilon_{e}}(\epsilon_{ e})}{\partial\epsilon_{e}}+\frac{\partial\psi^{-}_{\epsilon_{e}}(\epsilon_{e})}{ \partial\epsilon_{e}} \tag{24}\]  
The stiffness matrix of the material is:  
\[\mathbb{C}(\epsilon_{e},d)=\frac{\partial\sigma(\epsilon_{e},d)}{\partial \epsilon_{e}}=g(d)\frac{\partial^{2}\psi^{+}_{\epsilon_{e}}(\epsilon_{e})}{ \partial\epsilon^{2}_{e}}+\frac{\partial^{2}\psi^{-}_{\epsilon_{e}}(\epsilon _{e})}{\partial\epsilon^{2}_{e}} \tag{25}\]  
The stiffness matrix in Equation (21) exhibits pronounced nonlinearity due to the degradation function's dependence on the evolving phase-field variable. To overcome the limitations, the Hybrid phase field solution strategy is adopted [43; 44].