## 3 Numerical Model for Multi-Cluster Hydraulic Fracturing in Shale Reservoirs with Natural Fractures
### Modeling Methods for Hydraulic Fracturing  
It is assumed that the fluid flowing in fractures is incompressible. The flow is then split into tangential flow, which occurs along the fracture plane, and normal flow, which occurs perpendicular to the fracture plane, as shown in Figure 7. The tangential flow rate formula for the fractures can be expressed as [55]:  
\[q_{i}=-\frac{d^{3}}{12\mu}\
abla p \tag{3}\]  
where \(q_{i}\) is the flow velocity, \(d\) is the fracture width; \(p\) is the fluid pressure in the fracture, and \(\mu\) is the fluid viscosity. The flow pattern on both the upper and lower surfaces of the fracture can be described as follows [56]:  
\[\left\{ \begin{aligned} q_{t}&=c_{t}(p-p_{t})\\ q_{b}&=c_{b}(p-p_{b})\end{aligned} \right\} \tag{4}\]  
where \(q_{t}\) and \(q_{b}\) are the volume flow rates of the upper and lower surfaces of the fracture, respectively; \(c_{t}\) and \(c_{b}\) represent the filtration coefficient; \(p_{t}\) and \(p_{b}\) are the pore fluid pressures on both the top and bottom sides of the breakage, respectively. Hydraulic fracturing simulation is depicted as shown in Figure 7. Fracturing fluid is injected into the fracture through the wellbore, and under the effect of fluid pressure, the crack is opened. When the fluid pressure at the tip reaches the strength of the formation, the crack will propagate and begin to extend. In hydraulic fracturing simulation calculation, the element damage mode of the CZM follows the traction-separation law, and in its core lies the relationship between the traction force and separation displacement. The process of rock damage and hydraulic fracture propagation is governed by the traction-separation law. The traction force acting on a cohesive element can be expressed as:  
Figure 7: Diagram of fluid flow in a fracture. \[\left\{ \begin{array}{c}t_{n}\\ t_{s}\\ t_{t}\end{array} \right\}=\left[ \begin{array}{ccc}E_{mn}&E_{ns}&E_{nt}\\ E_{sn}&E_{ss}&E_{st}\\ E_{tn}&E_{ts}&E_{tt}\end{array} \right]\left\{ \begin{array}{c}\delta_{n}/T_{0} \\ \delta_{s}/T_{0}\\ \delta_{t}/T_{0}\end{array} \right\} \tag{5}\]  
where \(t_{n}\) is the normal stress exerted on the element, \(t_{s}\) and \(t_{t}\) are the shear stresses exerted on the element in two directions, \(E_{mn}\), \(E_{ns}\), \(E_{nt}\), \(E_{ss}\), \(E_{st}\), and \(E_{tt}\) are the six components of the stiffness matrix of the cohesive elements, \(\delta_{n}\), \(\delta_{s}\), and \(\delta_{t}\) are the normal and shear separation. \(T_{0}\) is the thickness of the initial cohesive element. In numerical simulations, the bedding plane is selected as the \(x\)-\(y\) plane, with the outward normal direction of the bedding plane parallel to the \(z\)-axis. According to this treatment, \(E_{ns}=E_{nt}=E_{sn}=E_{st}=E_{tn}=E_{ts}=0\). The other components of the stiffness matrix \(E_{mt}\), \(E_{ss}\), and \(E_{tt}\) could be obtained from the tensile moduli for the vertical bedding planes and parallel bedding planes as per Tables 4 and 5, respectively. \(T_{0}\) is the thickness of the initial cohesive element. As shown in Figure 8, the total fracture energy is represented by the area of the pink triangle, which can be divided into a yellow part (the shear component) and a green part (the normal component). When the fracture energy reaches the critical energy, the material will be completely failed and unable to support any traction force. As the initial damage criterion, a quadratic nominal stress criterion is widely used in hydraulic fracturing simulations [57]. As shown in Figure 8, \(t_{n}^{0}\) is the critical normal stress for the failure of the element; \(t_{s}^{0}\) and \(t_{t}^{0}\) are the critical tangential stresses of failure in the two directions of the element. The symbol \(<\!\!>\) signifies that the element can withstand tensile stress and does not undergo damage when subjected to compressive stress. Once the damage criterion is reached, the traction force decreases with the increase in the separation displacement according to:  
\[