## 2 Formulation and Theoretical Derivation of the Elastoplastic Phase Field Fracture Model
### Crack Surface Energy  
The fracture process in rock involves crack initiation and propagation. Consider a damaged medium occupying a domain \(\Omega\subset R^{n_{\text{dim}}}\) with spatial dimension \(n_{\text{dim}}\in\{1,2,3\}\), bounded by an external boundary \(\partial\Omega\) and containing an internal discontinuous surface \(\Gamma\), as illustrated in Figure 1a. As a non-local regularization technique, the phase field method (PFM) employs a diffuse crack model to transform sharp discontinuities from classical fracture mechanics into a continuous field description [19]. Specifically, the phase field variable \(\phi\) and its gradient approximate the crack topology. A diffuse crack band of finite width replaces the sharp crack surface, with crack evolution governed by the continuous distribution of \(\phi\):\(\phi\) = 0 corresponds to an intact material region, while \(\phi\) = 1 indicates a fully fractured state, as depicted in Figure 1b [20].  
Consistent with [PERSON]'s fracture theory, the fracture problem in PFM is formulated as a total system energy minimization problem, where the total system energy is a functional of the displacement field \(u\) and the phase field \(\phi\), reflecting the competition between elastic energy and crack surface energy [21; 22]. As a gradient-type non-local model, the fracture phase field theory eliminates mesh sensitivity issues. As a regularization approach, the phase field method does not require predefined crack propagation paths and can naturally capture crack initiation, propagation, branching, and merging. Under Dirichlet boundary conditions, the one-dimensional phase field is represented by an exponential function.  
\[\phi=\exp\biggl{(}-\frac{x}{\ell_{0}}\biggr{)} \tag{1}\]  
In this formula, \(\phi\) is the phase field variable, \(\phi\in[0,1]\), where \(\phi\) = 0 denotes an uncracked state and \(\phi\) = 1 signifies complete fracture; \(\ell_{0}\) is the crack regularization length scale, governing the smooth transition of the crack profile.  
Figure 1: Phase field crack characterization method. (**a**) Schematic of phase field crack characterization; (**b**) approximate characterization of crack topology.  
Following [PERSON] et al. [23], the crack surface energy can be expressed as  
\[E_{fra}=\int_{\Omega}G_{c}\gamma(f,\
abla f)\mathrm{d}\Omega \tag{2}\]  
In this formula, \(G_{c}\) is the critical fracture toughness; \(\
abla f\) is the spatial gradient of the phase field variable \(\phi\); and \(\gamma(\phi,\
abla\phi)\) is the regularized crack propagation energy term, defined as  
\[\gamma(\phi,\
abla\phi)=\frac{1}{2\mathcal{E}_{0}}\phi^{2}+\frac{\mathcal{E}_ {0}}{2}\
abla\phi\cdot\
abla\phi \tag{3}\]  
In Equation (3), the first term quantifies the energy required for local crack nucleation, and the second term introduces gradient regularization to ensure smooth crack front transitions. Substituting Equation (3) into Equation (2) yields the crack surface energy:  
\[E_{\text{fra}}=\int_{\Gamma}G_{c}\mathrm{d}\Gamma=\int_{\Omega}G_{c}\gamma( \phi,\
abla\phi)\mathrm{d}\Omega=\int_{\Omega}\frac{G_{c}}{2\mathcal{E}_{0}} \phi^{2}+\frac{\mathcal{E}_{0}}{2}G_{c}|\
abla\phi\,|\,^{2}\mathrm{d}\Omega \tag{4}\]