## 1 Introduction  
The successful and efficient exploitation of shale gas heavily relies on hydraulic fracturing to achieve commercial production [1], which is the most effective way to recover the gas from low permeability naturally fractured formations [2]. Creating a complex fracture network system is the most important goal of hydraulic fracturing in ultra-low permeability reservoirs. However, geological materials are discontinuous and heterogeneous due to the tectonic action of strata. Discontinuous joints and natural fractures that are randomly distributed in the reservoirs have a significant influence on the fracture network geometries. Observations of core samples and image logging indicates that discontinuities in shale, such as natural fractures and bedding interfaces, exist in most shale reservoirs and act as weak planes to affect the propagation of hydraulic fracture [3]. Field fracturing tests and microseismic monitoring reveal that hydro-induced fractures in hydraulic fracturing can activate natural fractures and weak planes in the reservoir, thus generating a large number of branch fractures, and significantly improving the permeability of the reservoir [4; 5; 6]. In the past decades, many studies have been conducted to investigate the effects of discontinuities (e.g., natural fractures and joints) on perturbing the propagation of hydraulic fractures, and the results reveal that natural fractures can observably improve the ultimate morphology of the fracture network [7; 8; 9]. For instance, [PERSON] and [PERSON] [10] used microseismic monitoring technology to image hydraulic fracturing fractures in shale, revealing a positive correlation between gas production and fracture network complexity. [PERSON] and [PERSON] [11] established an orthogonal intersection criterion of hydraulic fractures and natural fractures to predict whether the hydrofractures in brittle rocks would propagate through natural fractures or terminate after they intersected vertically. Later, [PERSON] [12] supplemented the orthogonal intersection criterion and found that a hydraulic fracture was more likely to extend along the interface of a natural fracture when their intersection angle was less than 90\({}^{\circ}\), which has been verified by experiments. [PERSON] [13] reported that there might be bypass, separation of weak planes, fracture deflection, and mixed expansion of the above three modes at the intersection of hydraulic fractures and natural fractures. They concluded that the presence of these natural fractures can have beneficial or detrimental effects on the complex geometry of hydrofractures. On the one hand, the activation of natural fractures may expand the flow paths of oil and gas and improve fracture complexity, making stimulation treatments more effective [14]. On the other hand, reactivated natural fractures may result in additional leak-off and reduced fluid pressure in fractures, hence a lower treatment fluid efficiency [15]. Despite the existing research results, most experiments or numerical model have been devoted to investigating the interaction of a single natural fracture or several hydraulic fractures but ignore the natural fractures in the form of a network. A thorough understanding of the fracture propagation behavior in a natural fracture network is still lacking. In recent years, some studies have focused on large-scale fracture network simulation. New hydraulic fracturing models have been developed and are still evolving. So far, the main methods used to simulate complex fracture networks are the Displacement Discontinuous Method (DDM) [16; 17; 18; 19], Boundary Element Method (BEM) [20], Finite Element Method (FEM) [21; 22], Extended Finite Element Method (XFEM) [23; 24], Discrete Fracture Network (DFN) [25; 26; 27], and the Discrete Element Method (DEM) [28; 29].