Fractures occur at varying scales and orientations in the subsurface. The role of fractures as conduits for fluid flow in a reservoir must be well constrained for planning and development of a geothermal system. This study examines the uncertainty associated with fractures in a rock matrix using a discrete fracture network modelling approach; fractures are considered as discrete elements embedded in a rock matrix. Many fractures are below the resolution of geophysical investigation methods and are therefore typically described by a statistical distribution of permeability, length and orientation. A Monte Carlo approach is used whereby multiple reservoir simulations are conducted using fracture network parameters drawn randomly from a pre-defined distribution. The simulation results show that fracture permeability is the major factor which influences the production. An increase in fracture network permeability from 8:33x10-^12 m2 to 8:33x10-^10 m2 showed an increase of uncertainty in reservoir production by up to a factor of six when considering the standard deviation. The results show that this uncertainty increases by up to 200 % with the life of the reservoir for the high permeability case. Stochastic analysis of the azimuth uncertainty revealed an opposite trend where a decrease in uncertainty of 34 % was observed over the 30 year simulated lifetime of the reservoir. The orientation of fractures in the reservoir was found to play an important role in the temperature field in the reservoir with increases in production temperature ranging from 0:5 K to 3 K resulting from the proximity of fractures to the production well. The results also indicate that discrete fracture networks which possess similar orientations of fractures produce statistically similar production profiles. The findings of this study show that a stochastic approach can be used for efficient planning of geothermal doublet systems based on expected trends in production.