In every tight formation reservoir, natural fractures play an important role for mass and energy transport and stress distribution. Enhanced Geothermal Systems (EGS) make no exception and stimulation aims at increasing the reservoir permeability to enhance fluid circulation and heat transport. EGS development relies upon the complex task of predicting accurate hydraulic fracture propagation pathway by taking into account reservoir heterogeneities and natural or pre-existing fractures. In this contribution, we employ the variational phase-field method which handles hydraulic fracture initiation, propagation and interaction with natural fractures and is tested under varying conditions of rock mechanical properties and natural fractures distributions. We run bi-dimensional finite element simulations employing the open-source software OpenGeoSys and apply the model to simulate realistic stimulation scenarios, each one built from field data and considering complex natural fracture geometries in the order of a thousand of fractures. Key mechanical properties are derived from laboratory measurements on samples obtained in the field. Simulations results confirm the fundamental role played by natural fractures in stimulation’s predictions, which is essential for developing successful EGS projects.