Geothermal energy is a promising source of clean renewable energy that does not depend on weather conditions and has the potential of providing base load as well as peaking electrical supply. During the last twenty years, application of the geothermal energy for electricity generation, heating and cooling have increased noticeably in many countries, including USA, Iceland, New Zealand, Turkey, and Indonesia. The worldwide installed geothermal power capacity is approximately 16 MWe, with a growth rate estimated between 3.4% to 5.4% per year over the period between 2015 to 2060, according to the World Energy Council [1]. Geothermal energy has the potential to provide 3% of electricity and 5% of heating with respect to the global demand by 2050 [2].
 Enhanced Geothermal System enables producing geothermal energy from Hot Dry Rock (HDR) reservoirs, which are deficient in water and permeability, the two basic components necessary to exploit their geothermal potential (e.g., [3, 4]). EGS implementation requires the drilling of injection and production wells. Through the first one a cold fluid is injected, whereas hot water and/or steam is recovered from the production well. A key feature of this methodology is the network of fractures that connect these two wells. The discontinuities can be either, (naturally) pre-existing fractures, or (artificially) triggered fractures by the thermal shock induced by the contact between the cool injection fluid and the hot natural rock. Water or liquid nitrogen are contemplated as potential injection fluids. This type of project envisages reservoirs at depth with temperature above 180°C. Several field tests in different countries have been conducted to evaluate the feasibility of EGS (e.g., [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]), including the ongoing Utah-FORGE project in the USA (e.g., [13]).
 Numerical modeling of the interaction between hydraulic driven fractures (HF) and natural fractures (NF) is still a challenging problem, since the presence of NFs introduce heterogeneity on the hydromechanical properties of the rock mass and affects the geometry of the HF. Therefore, the development of robust numerical methods able to address this kind of phenomena is crucial for the safe application of stimulation techniques. This work presents a finite element (FE) technique, called mesh fragmentation technique (MFT) [14], capable of tackling this type of problem. The MFT consists in introducing solid high-aspect ratio finite elements in-between the regular (bulk) finite elements to simulate both NFs and the formation and propagation of HFs. This MFT has been successfully used to model drying cracks in soils [15, 16], hydraulic fracturing in rocks [17, 18] and natural fractures deformation due to pressure depletion [19]. The fully coupled hydro-mechanical approach has been implemented in the FE program CODE_BRIGHT. The paper presents the main theoretical components and implementation of the proposed approach. Numerical examples are used to demonstrate the capabilities of the proposed technique. Figure 1 presents the pressure contours related to the modeling of the hydraulic fracturing of two neighbor wells in a naturally fractured rock. A very satisfactory performance of the proposed method is observed in all the analyzed cases.