Hydraulic fracturing is currently the primary method for stimulating low-permeability geothermal reservoirs and creating Enhanced (or Engineered) Geothermal Systems (EGS) with improved permeability and heat production efficiency. Complex natural fracture systems usually exist in the formations to be stimulated and it is therefore critical to understand the interactions between existing fractures and newly created fractures before optimal stimulation strategies can be developed. Our study aims to improve the understanding of EGS stimulation-response relationships by developing and applying computer-based models that can effectively reflect the key mechanisms governing interactions between complex existing fracture networks and newly created hydraulic fractures. In this paper, we first briefly describe the key modules of our methodology, namely a geomechanics solver, a discrete fracture flow solver, a rock joint response model, an adaptive remeshing module, and most importantly their effective coupling. After verifying the numerical model against classical closed-form solutions, we investigate responses of reservoirs with different preexisting natural fractures to a variety of stimulation strategies. The factors investigated include: the in situ stress states (orientation of the principal stresses and the degree of stress anisotropy), pumping pressure, and stimulation sequences of multiple wells.