Enhanced geothermal system is currently the most efficient technology of exploiting geothermal energy from high-temperature rock reservior in deep earth. However, due to the uncertainty of geological structural discontinuities, such as natural or artificial fractures, great challenges are encountered in conducting the precise evaluation of heat extraction within enhanced geothermal system. In this paper, a probability analysis is carried out to incorporate the uncertainty of the fracture network into the analysis of the spatial-time evolution on the seepage and temperature fields, based on a discrete fracture network based modeling scheme developed by the authors. The developed finite element model of the fractured reservior exhibits two-phase structure consisting of rock matrix and fractures, which are practically discretized into triangular elements and zero-thickness elements, respectively. The approach demonstrates a favorable comparison with the results obtained from the thereotical and other numerical solution. Monte Carlo simulation technique is then employed to probabilistically analyze the evolution of coupled thermal–hydraulic behavior, which involves 1000 different realizations considering the uncertainty of fracture structure. The results illustrate that the coefficient of variation of the seepage field decreases from 0 to 0.29 at 1 a to 0.06 to 0.15 at 50 a, averaging at 0.07, and the peak coefficient of variation of the temperature field decreases from 0.61 to 0.52, averaging at 0.27. The uncertainty of the temperature field is varying against time and prominent around the proximity to the advancing cold front, upper and lower boundaries, primarily associated with variations in fractures and constraint of boundary conditions. The existence of highly permeable fractures leads to a pronounced heterogeneity in the seepage and temperature fields within enhanced geothermal system. Furthermore, the importance of taking the uncertainty of the fracture structure into account for geothermal reservoir evaluation is highlighted, to enhance the reliability of numerical modeling for geothermal process prediction.