The hydraulic fracturing propagation is strongly influenced by the existence of natural fractures. This is a very important factor in hydraulic fracturing operations in unconventional reservoirs. Various studies have been done to consider the effect of different parameters such as stress anisotropy, toughness, angle of approach and fluid properties on interaction mechanisms including crossing, arresting and opening. Analytical solutions can only be used for simple fracture geometries and are not usually able to provide good predictions due to many simplified assumptions. Laboratory tests are also conducted under certain constraints like sample size and conditions that are different from the real field conditions. Numerical simulations, including continuum and dis-continuum based models have been used extensively to simulate hydraulic fracture propagation and its interaction with natural interfaces. However, calibration of simulated models with real field data is necessary to ensure the accuracy of the results. A calibrated numerical simulation can be used to model complex geometries. In this study, a Lattice numerical simulator, which is the advanced version of Particle flow Code (PFC) based on the granular particle physics, was used for numerical simulation of lab scale hydraulic fracturing. The scaling laws were also used to increase the dimensions of the simulated samples to allow increasing the rate of fluid injection and reducing its viscosity, hence reduce the simulation time. The interaction of hydraulic fractures and orthogonal fractures with angles of approach of 90°, 60° as well as non-orthogonal fracture planes with different filling materials ranging from strong to very weak were studied. The results showed good agreement with lab observations. In general the larger the angle of approach and stronger the filling material, the higher the likelihood of the crossing mode. Also, networks of regular natural fractures with two fracture sets were simulated. The results showed that the combination of different parameters define the preferred fracture propagation (PFP) which is not easy to predict using analytical solutions. In this situation and more complex real field cases, the use of numerical simulations are necessary to predict the propagation of hydraulic fracture and interaction modes.