## 2 Numerical Model
### Simulation Method
With continued fluid injection, hydraulic fractures will gradually emerge within the reservoir rock mass, a phenomenon that can be simulated by modeling the damage and failure of interface elements. When the interface element is damaged, the anti-deformation ability of the interface element will significantly decrease, leading to the rapid opening of hydraulic fractures. The presence or absence of interface elements can be determined using the following formula [34]:  
\[\left\{\frac{\langle\sigma_{n}\rangle}{\sigma_{n}^{0}}\right\}^{2}+\left\{ \frac{\sigma_{s}}{\sigma_{s}^{0}}\right\}^{2}+\left\{\frac{\sigma_{t}}{\sigma_ {t}^{0}}\right\}^{2}=\lambda \tag{6}\]  
where \(\sigma_{n}\) is the normal stress; \(\sigma_{s},\sigma_{t}\) represent the tangential stress, MPa; \(\sigma_{s}^{0},\sigma_{t}^{0}\) are the threshold stress, MPa; \(\lambda\) indicates that cohesive elements resist tension stress but not compression stress: \(1\leq\lambda\leq 1.05\). After the hydraulic fractures are initiated, the fractures will not continue to open indefinitely due to factors such as interfacial forces. The deformation behavior after reaching the peak can be described by damage factors [28; 35],  
\[E=(1-d)\times E^{0} \tag{7}\]  
\[d=\frac{\delta_{m}^{f}\big{(}\delta_{m}^{max}-\delta_{m}^{0}\big{)}}{\delta_{ m}^{max}\big{(}\delta_{m}^{f}-\delta_{m}^{0}\big{)}}, \tag{8}\]  
where \(E_{0}\), \(E\) are the initial elastic modulus and the elastic modulus after damage, Pa, respectively. Here, \(d\) is a damage factor, dimensionless. \(\delta_{m}^{max}\), \(\delta_{m}^{f}\), \(\delta_{m}^{0}\) signify the distinct displacements in the context of element behaviour. Specifically, \(\delta_{m}^{max}\) refers to the maximum displacement experienced. \(\delta_{m}^{f}\) and \(\delta_{m}^{0}\) denote the displacement. When performing multi-cluster fracturing, the injected fluid will flow into the formation through the wellbore, resulting in flow distribution within the wellbore. Meanwhile, due to the long wellbore, there is also a corresponding pressure drop. \[Q=\sum_{i=1}^{N}Q_{i} \tag{9}\]  
\[p_{0}=p_{pf,i}+p_{cf,i}+p_{wf,i} \tag{10}\]  
When fluid flows into the wellbore, flow resistance is easily generated due to the roughness of the wellbore's surface, and this behavior can be described by [PERSON]'s equation [36]:  
\[ \begin{array}{l}\
abla p-\rho g\Delta Z=(C_{L}+K_{i})\frac{\rho v^{2}}{2}, \\ C_{L}=\frac{fL}{D_{h}},\\ f=8\Big{[}\big{(}\frac{8}{\mathrm{Re}}\big{)}^{12}+\frac{1}{(A+B)^{13}} \Big{]}^{1/12},\\ A=\Big{[}-2.457\mathrm{Im}\Big{(}\big{(}\frac{7}{\mathrm{Re}}\big{)}^{0.9}+0.2 7\Big{(}\frac{K_{i}}{D_{h}}\Big{)}\Big{)}\Big{]}^{16},\\ B=\big{(}\frac{37350}{\mathrm{Re}}\big{)}^{16}\\ \end{array}  \tag{11}\]  
\[p_{f}=0.087249\times\frac{\rho}{n^{2}D_{p}^{4}C^{2}}q_{i}^{2} \tag{12}\]  
where \(\Delta P\) is the pressure difference at the node of the cluster, \(\Delta Z\) is the elevation difference of the node, \(v\) is the fluid velocity in the wellbore, \(\rho\) is the fluid density, \(G\) is the acceleration of gravity, \(C_{L}\) is the loss coefficient, \(f\) is the friction force on the wellbore, \(L\) is the wellbore length, \(K_{i}\) is the loss term in a fixed direction, \(K_{s}\) is the roughness of the wellbore, \(D_{h}\) is the wellbore diameter, Re is the Reynolds number, \(n\) is the number of perforations, \(D_{p}\) is the perforation diameter, and \(C\) is a coefficient, generally 0.56-0.9 [6]. After the fracturing fluid is diverted from the wellbore, it will enter the fracturing fracture, and the fluid pressure drop inside can be described as  
\[q=\frac{t^{3}}{12u}\
abla p, \tag{13}\]  
where \(q\) is the volume flow vector per tangential element length; \(t\) is the thickness, m; \(\mu\) is the viscous coefficient; \(p\) is the fluid pressure, MPa. It is worth mentioning that there is a certain phenomenon of crossflow between the fluid used for fracturing and the existing fluid in the reservoir.