## 1 Introduction
Each of these methods has advantages and disadvantages: the DDM has some difficulties in the deformation of heterogeneous and anisotropic materials; the DEM has a high computing cost in the computing domain; the FEM is an effective and widely used numerical technique, but it is time consuming due to grid repartition. The Cohesive Zone Method (CZM), an efficient tool for fracture initiation and propagation modeling that has been developed in recent years, can help to solve these problems [30; 31; 32; 33]. In the CZM, cohesive elements with pore pressure DOF and displacement DOF are used to describe hydrofractures with consideration of fluid flow in the fracture and leak-off to the rock [34; 35], as shown in Figure 1. The fracture process is characterized by the Traction-Separation Law, which mainly depends on the material properties of the cohesive zone. The CZM successfully circumvents the significant perturbance that the singularity at the fracture tip causes in the solution convergence of the equations. In addition, the CZM imposes restrictions on the hydraulic fractures' propagation paths and only allows fractures to extend along the edges of cohesive elements. Hence, the FEM in conjunction with the cohesive layer can explicitly model the simulation of the intersection between hydraulic fractures and natural fractures without grid remeshing, which reduces the computational expense [36]. However, this propagation restriction is also the main shortcoming of CZM. It requires that the direction of the propagation of the fracture must be known or actively updated during propagation. In other words, the intersection at a natural fracture is predefined and not predicted. Three types of interactions are observed between a hydraulic and a single natural fracture [15], in which we can preset the expansion path of fractures, but this cannot be achieved in the natural fracture network due to the much more complicated fracture geometries. A key consideration of a complex fracture model is the propagation of a hydraulic fracture in a connective natural fracture network. In this case, fracture propagation may occur along the reactivated natural fracture paths with complex geometries, including deflections and bifurcations. In this way, the intersection mechanism between the fracture and the single natural fracture is not so important, but the influence of the fracture network should be considered. In field tests or experimental tests, the development process of induced fractures in rocks is actually unknown, making it hard to understand the role that natural fractures play in induced fracture propagation. In this paper, we investigated the step-by-step propagation pattern of hydraulic fractures in a natural fracture network and realized the non-prefabricated propagation of hydraulic fractures by using global cohesive elements. The effects of multiple factors such as in situ horizontal stress, rock tensile strength, and natural fracture cementation strength were discussed. The results presented in this paper could provide a visible propagation process of a hydraulic fracture and improve the basic understanding of the fracture development in a natural network. It helps the prediction of a complex fracture network in a real situation.