Research presented in this paper aimed at establishing a better understanding of natural fracture (NF)/hydraulic fracture (HF) interaction mechanisms and fracture network development in naturally fractured and nonhomogeneous shale formations through numerical modelling using the two-dimensional particle flow code (PFC2D). Hydraulic fracture propagation was firstly modelled in a 30 m x 30 m model representing intact rock by bonded particle method (BPM), which served as a base case in the research. Then a single, deterministic natural fracture was embedded into the same model by a smooth joint contact model (SJM) to investigate different NF/HF interaction mechanisms under a range of different conditions by varying the angle of approach, differential horizontal stress, and the mechanical properties of a fracture within the model. Based on the parametric research findings, number and diversity of natural fractures in the model were increased both deterministically and stochastically, and the results are compared and discussed.