Transition to cleaner energy while sustaining the heat demand necessity has always been an emphasized topic in the recent years. With the slowly declining utilization and production of natural gas in Netherlands, geothermal energy shows a promising future in providing sustainable power for heating and various applications. For majority of economically feasible geothermal project, either natural fractures or hydraulic fractures are present. These fractures served as a preferential fluid pathway between geothermal well couplets, which determines crucial parameters such as breakthrough time and project life. In this study, the main focus is the impact of stochastic distribution of fracture aperture with or without pressure dependence and gravity effect on geothermal energy production. The simulation conditions as well as reservoir physical model are collected from analogous enhanced geothermal systems (EGS)that are classified as deep geothermal energy projects. The test cases include base case, case with no pressure or stress dependence, case with half well doublet distance, and case with gravity that has well positioned in the dip or strike direction in relative to the fault orientation. The ensemble simulation of 100 realization is run for each scenario considered. The relevant parameters analyzed are pressure difference between well couplet, production temperature, NPV, and most importantly, net cumulative energy production. It is found out that the geothermal energy production greatly relies on the fracture aperture distribution with or without pressure dependence. For base case, the range of energy production is from 1.50E13 to 3.75E13 KJ with temperature variation of 35 K between the minimum and maximum value. Without pressure dependence, only pressure difference between injector and producer slightly changes, and production temperature as well net energy pro-duction remains nearly identical. For the case with half of original well-spacing, the spread of energy production distribution reduces from 1.50E13 to 3.25E13 KJ. The diminishing energy production results from less total enthalpy in the flow path with shorter distance, yet the outcome of energy production still demonstrates a considerable range. By taking account of gravity effect and populate the reservoir with non-uniform pressure and temperature distribution with corresponding gradient, the net energy production drops to approximately two times less than the base case. However, the ratio of maximum energy production to the minimum over100 realization is similar to that of the base case. Last but not least, two types of well orientation, along the strike and along the dip of fracture plane, are considered to observe the impact of gravity-assisted flow. It is concluded that with the current model setting, the low temperature gradient of injector in the dip configuration case overweighs the effect of gravity assistance, resulting in a smaller net energy production than the strike case. This project concludes that the fracture aperture heterogeneity is an important factor in predicting the outcome of deep geothermal production. The gravity is another crucial component to portray accurate flow behavior between the wells, but it does not result in additional increase in the range of energy production outcome in comparison to the base come with the simulation conditions applied.