This paper presents a three-dimensional (3D) numerical model of multifracture hot dry rock (HDR) systems using coupled thermo-hydro-mechanical (THM) modelling processes. Several numerical models are implemented in the COMSOL Multiphysics Finite Element (FE) solver to analyse the effect of fracture spacing and number on reservoir productivity. The developed multifracture HDR reservoirs employ vertical fracture systems with deviated wellbores to reflect the state of stress in deep geological formations. The models consider the geomechanical effect of fracture aperture, permeability and stiffness changes on each individual fracture during long-term exploitation of 30 years. The crucial parameters investigated concerning reservoir productivity are production temperature, thermal energy and recovery rate. The results show that both fracture number and spacing affect the productivity of HDR geothermal systems. It is clear from the results that variations of geomechanical parameters of fracture aperture, permeability and stiffness can play a crucial role in enhancing the performance of these systems.