With the uncertain fractures involved in the water-rock interaction in the enhanced geothermal system (EGS) taken into account, numerical simulation of heat extraction is challenging to implement, especially in the situation where the EGS is conducted in a complex fracture network. In this study, a discrete fracture network (DFN) modeling procedure based on finite element method (FEM) is proposed to simulate the uncertainty of heat extraction performance in the fractured hot dry rock (HDR), in which the rock matrix and fractures are discretized into triangular elements and zero-thickness elements, respectively. To verify the effectiveness of the established numerical model, the simulation results are compared with the analytical solution and a good agreement is reached. Further, Monte Carlo simulation (MCS) method is adopted to probabilistically study the production performance of EGS considering the randomness of the fracture network. According to the MCS results, heat production performance, i.e., outlet average temperature, EGS lifetime, output thermal power, heat extraction ratio, and levelized cost of energy (LCOE), are statistically analyzed. The study results indicate that the heat extraction performance can be more effectively assessed by the probabilistic analysis, and the uncertainties involved in the heat extraction process cannot be ignored. •Fracture uncertainty should be taken into account in water-rock interaction in EGS.•Single deterministic simulation might be insufficiently accurate.•MCS method is adopted to probabilistic study the production performance.•Probabilistic analysis provides guidance for geothermal exploitation in EGS.•Performance assessment parameters data have statistical distribution feature.