Effects of fracture surface roughness on the coupled hydro-thermal processes in rocks with single fracture have been systematically studied, whereas the variability of fracture aperture in a three-dimensional(3D) discrete fracture network(DFN) is generally neglected. The present work conducted to investigate the influences of stochastic fracture surface roughness in the 3D macroscopic water flow and heat transfer in fractured rocks. Two different Joint Roughness Coefficient(JRC) distribution functions were fitted on the basis of the published JRC data at Bakhtiary dam site and Oskarshamn/Forsmark mountain, and two different fracture hydro-mechanical aperture relationship models were used to calculate the fracture hydraulic aperture. The role of JRC distribution varieties and hydro-mechanical aperture relationship diversities in flow rate and temperature distribution in fractured rocks were numerically analyzed by using COMSOL Multiphysics finite element software. Numerical results indicate that: Lognormal distribution and normal distribution are adopted to appropriately describe the JRC distribution at Bakhtiary dam site and Oskarshamn/Forsmark mountain, respectively; Fractured rocks with JRC lognormal distribution are more permeable than those with JRC normal distribution; Li & Jiang model provides a greater fracture hydraulic aperture than Barton model resulting in faster moving cold front, the water outflow, hydraulic pressure and temperature in Barton model are more sensitive to JRC distribution.