## 1 Introduction  
The role of shale oil and gas resources in the world's energy structure is becoming increasingly prominent [1]. However, the production of shale oil and gas faces some challenges due to the low permeability of shale reservoirs, which sets it apart from conventional reservoirs. One of the main methods to boost shale gas output is hydraulic fracturing, which has been proven to be highly effective in increasing the output of shale gas [2]. Some of the most common simulation methods for fracture propagation in a shale reservoir include the extended finite element method (XFEM) [3; 4], discrete element method (DEM) [5], displacement discontinuity method (DDM) [6; 7], and finite element method (FEM) [8]. The extension of hydraulic fractures is often simulated using the cohesive zone method (CZM), which relies on the principles of viscoelastic-plastic damage theory and is widely applied in finite element method simulations. [PERSON] et al. [9] studied the hydraulic fracturing procedure using the CZM and established two hydraulic fracturing models. This includes the Khristianovic-Geersma-De Klerk (KGD) problem as well as the challenges related to the interaction between hydraulic and natural fractures. A meshing method has been proposed for the fracture intersection problem. [PERSON] et al. [10] developed a hydraulic fracturing model that includes both the deformation of a solid rock and the fluid flow within fractures, with a focus on 2D aspects. [PERSON] et al. [11] built a 2D numerical model of hydraulic fracturing and used it to investigate the impact of model size on the propagation of fractures. [PERSON] et al. [12] used CZM elements to simulate the change in the hydraulic fracture propagation under different reservoir porosity conditions and fracturing conditions. [PERSON] et al. [13] studied the propagation behavior of 3D hydraulic fractures and found that the ground stress conditions affect the height of hydraulic fractures. To understand the interference phenomena of multiple clusters of hydraulic fractures, [PERSON] et al. [14] and [PERSON] et al. [15] studied the deflection phenomenon of multiple fracture propagation due to interference. [PERSON] et al. [16] proposed a 3D model of fractures, taking into account factors such as fluid flow within the fractures [17; 18]. [PERSON] et al. [19], [PERSON] et al. [20], and [PERSON] et al. [21] simulated the spatially extended perturbation behavior of a fracturing network in 3D hydraulic fracturing with multi-cluster fractures and different shot-hole cluster spacings. After an extensive review of existing fracture propagation simulation methods, the cohesive element method has been widely applied in simulating the fracture behavior of materials, including in shale reservoir fracture propagation simulations, etc. Hydraulic fracture propagation and growth depend on the reservoir's tensile and compressive strength properties. In terms of research on the compressive performance of shale, [PERSON] et al. [22] and [PERSON] et al. [23] conducted uniaxial compression tests on stratified Longmaxi shale in Sichuan, finding significant anisotropy in the uniaxial compression parameters of shale, with the bedding plane direction having a considerable impact on the stress-strain curve of shale. The uniaxial compressive strength of shale increases with the increase in the bedding dip angle. [PERSON] et al. [24] conducted triaxial compression tests on Longmaxi Formation shale in Sichuan, discovering that the bedding plane direction significantly influences the elastic modulus of shale. The elastic modulus of shale vertical to the bedding plane is larger than that parallel to it. Different core angles vertical to the bedding direction exhibit strong anisotropy, but the anisotropy is not obvious in the direction parallel to the bedding. They proposed that shale can be approximated as a transversely isotropic medium.