Geo-energy is a comprehensive term used to describe any form of energy that comes from the Earth. This includes hydrocarbons such as gas, oil, and coal, but also geothermal energy (shallow and deep). The focus of this thesis is on Enhanced Geothermal Systems (EGS), which are renewable sources of energy. In EGS, hydraulic stimulation is the technique employed to enhance the permeability e.g. of the granite basements at great depths, which might exceed 5 kilometers. Rocks have natural fractures which affect the behavior of the rock mass. The mechanical goal of hydraulic stimulation is to open these pre-existing fractures. Due to the generally great depths, only indirect information such as in-situ stress, pumping records, and micro-seismicity, are available. Many issues in the field of reservoir geothermal (but also petroleum) production, are attributed to the lack of proper understanding of how these fractures interact with each other. The main goal of this study is to provide an insight into how the pre-existing fractures in rocks interact with each other. Specifically, this research answers the question of how pre-existing flaws of different geometrical parameters contained in rock specimens subjected to various loading configurations, affect the mode and the location of the mechanism of crack initiation. We develop a finite element (FE) linearly elastic MATLAB Partial Differential Equation (PDE) code. This numerical model predicts the mode and the location of the initiation of the local fractures at the inner tips of two interacting pre-existing flaws in Barre granite in the context of EGS. We investigate different combinations of double-flaw geometries (flaw inclination angle β, bridging angle α, and ligament length L) contained in specimens subjected to various loading scenarios (Combination of Vertical Load VL, Lateral Load LL, and Internal Flaw Pressure FP) likely to induce first local tensile and shear fractures, and we determine their corresponding predicted locations around the inner flaw tips. For a given double-flaw pair geometry, we subject the specimen to a certain loading configuration (Combination of VL, LL, and FP). We then retrieve the state of stress (maximum and minimum principal stresses 𝜎1 and 𝜎3, and maximum shear stresses 𝜏𝑚𝑎𝑥) at the Finite Element (FE) corner nodes on the boundary of the inner tips, from which we plot the corresponding Mohr’s circles. We study the interaction of these Mohr’s circles with the Griffith-Coulomb failure envelope. Since the main objective is to predict the first local failure’s mechanism and location, we iteratively vary the magnitude of the loading component until the first Mohr’s circle touches the failure envelope. Depending on where it does, the mechanism of the local fracture is inferred. If the maximum principal stress (𝜎1) of this Mohr’s circle is equal to the tensile strength of the granite (𝜎𝑇=7.5𝑀𝑃𝑎), a local tensile fracture is numerically predicted. If this Mohr’s circle touches the linear shear Coulomb portion of the envelope, it is a predicted local shear failure. The corresponding location of the FE node is then determined at the inner tip of the flaw, and this is where for the given geometry and loading scenarios, the first local fracture’s initiation is numerically expected. We follow this procedure for several double-flaw geometrical configurations. We compare the numerical predictions with the previously obtained experimental results in the MIT CEE Rock Mechanics Laboratory. This allows us to validate our set-up numerical model. 4 This study provides a solid ground to the understanding of the initiation and the interaction of pre-existing fractures in EGS. This thesis comes as an answer to a problem of major practical importance in the field of rock mechanics, and beyond. Typical geotechnical problems, such as rock slope stability and tunnel support, can also be better addressed when clearly understanding the process of crack initiation in rocks. This study also enables the correlation between numerical models and experimental findings, a challenge that has always existed in the research world. This thesis shows how our developed FE code can be used as a tool to predict the location and the mode of fractures in the laboratory, but also in the real world.