## 2 Numerical Model
### Simulation Method  
Multi-cluster hydraulic fracturing in reservoirs is a typical fluid structure coupling process that involves the interaction between fluids and solids. When considering the interaction between fluids and solids during hydraulic fracturing, the stress balance equation of hydraulic fracturing can be expressed as [31]  
\[\int\limits_{V}(\overline{\sigma}-p_{w}\mathbf{I})\delta\varepsilon\delta dV= \int\limits_{S}\mathbf{t}\cdot\delta\sigma dS+\int\limits_{V}\mathbf{f}\mathbf{v}dV \tag{1}\]  
where \(\overline{\sigma}\) represents the effective stress matrix, MPa; \(p_{w}\) denotes the pore pressure, MPa; \(\delta\varepsilon/s^{-1}\) is the virtual strain rate matrix; \(\mathbf{t}\) is the surface force matrix, \(\mathbf{N}/\text{m}^{2}\); \(\delta_{\text{v}}\) (m/s) and \(\mathbf{f}\) (N/m\({}^{3}\)) are the virtual velocity matrix and the physical force matrix. During the fracturing process, fluid density and matrix porosity influence the behavior of the fluid within the reservoir. The fluid continuity equation captures the conservation of mass during this process, providing a mathematical representation of the fluid's movement and distribution. It enables a more accurate description of the fluid flow dynamics within the reservoir [32]. \[\frac{d}{dt}\left(\int\limits_{V}\left(\frac{\rho_{w}}{\rho_{w}^{0}}ndV\right) \right)+\int\limits_{S}\frac{\rho_{w}}{\rho_{w}^{0}}\mathbf{n}\cdot v_{w}dS=0 \tag{2}\]  
where \(J\) represents the volume change ratio, dimensionless; \(\rho_{w}\) is the fluid density, kg/m\({}^{3}\); \(n_{w}\) is the void ratio, dimensionless; and \(v_{w}\) is the fluid seepage velocity, m/s. It is important to mention that the fracturing process also involves the influence of fluid seepage within the matrix. The seepage behavior of reservoir rocks during fracturing can be calculated using [PERSON]'s law, which determines the velocity and flow rate of fluids through factors such as fluid viscosity, permeability, and pressure gradient [33]:  
\[\mathbf{v}_{w}=-\frac{1}{n_{w}\mathbf{g}\rho_{w}}\mathbf{k}\cdot\left(\frac{\partial p_{w }}{\partial\mathbf{x}}-\rho_{w}\mathbf{g}\right) \tag{3}\]  
where \(\mathbf{k}\) represents a permeability matrix, m/s, and \(\mathbf{g}\) is the gravity acceleration vector, m/s\({}^{2}\). As fluid starts to seep into the reservoir rock mass, the influx of fluid will lead to an elevation in pore pressure, consequently facilitating the propagation of microfractures. It is generally believed that before hydraulic fractures open, the cohesive elements used to simulate hydraulic fractures follow a linear elastic relationship [6],  
\[\mathbf{\sigma}_{coh}= \begin{cases}\sigma_{coh\_n}\\ \sigma_{coh\_s}\\ \sigma_{coh\_t}\end{cases} =\mathbf{K}_{coh}\mathbf{\varepsilon}_{coh}= \begin{bmatrix} K_{coh\_nn}&K_{coh\_ns}&K_{coh\_nt}\\ K_{coh\_ns}&K_{coh\_ss}&K_{coh\_st}\\ K_{coh\_nt}&K_{coh\_st}&K_{coh\_t}\end{bmatrix}\begin{cases}\varepsilon_{coh\_s} \\ \varepsilon_{coh\_s}\\ \varepsilon_{coh\_t}\end{cases} , \tag{4}\]  
where \(\mathbf{\sigma}_{coh}\) denotes the stress vector; \(\mathbf{\sigma}_{coh\_h}\), \(\mathbf{\sigma}_{coh\_s}\), \(\mathbf{\sigma}_{coh\_t}\) are the normal stress, the first tangent stress and the second tangent stress, respectively. Here, \(\mathbf{K}_{coh}\) denotes the stiffness matrix; \(\mathbf{\varepsilon}_{coh}\) is the strain matrix, and \(\mathbf{\varepsilon}_{coh\_n}\), \(\mathbf{\varepsilon}_{coh\_s}\), \(\mathbf{\varepsilon}_{coh\_t}\) represent the normal strain, the first tangential strain, and the second tangential strain, respectively. They are defined as follows [25]:  
\[\varepsilon_{coh\_n}=\frac{d_{n}}{T_{0}},\varepsilon_{coh\_s}=\frac{d_{s}}{T_{ 0}},\varepsilon_{coh\_t}=\frac{d_{t}}{T_{0}} \tag{5}\]where \(d_{n}\), \(d_{s}\), \(d_{l}\) represent the normal displacement, the first displacement, and the second tangential displacement, respectively. Here, \(T_{0}\) is the constitutive thickness.