Efficient heat exploitation strategies from geothermal systems demand for accurate and efficient simulation of coupled flow-heat equations on large-scale heterogeneous fractured formations. While the accuracy depends on honouring high-resolution discrete fractures and rock heterogeneities, specially avoiding excessive upscaled quantities, the efficiency can be maintained if scalable model-reduction computational frameworks are developed. Addressing both aspects, this work presents a multiscale formulation for geothermal reservoirs. To this end, the nonlinear time-dependent (transient) multiscale coarse-scale system is obtained, for both pressure and temperature unknowns, based on elliptic locally solved basis functions. These basis functions account for fine-scale heterogeneity and discrete fractures, leading to accurate and efficient simulation strategies. The flow-heat coupling is treated in a sequential implicit loop, where in each stage, the multiscale stage is complemented by an ILU(0) smoother stage to guarantee convergence to any desired accuracy. Numerical results are presented in 2D to systematically analyze the multiscale approximate solutions compared with the fine scale ones for many challenging cases, including the outcrop-based geological fractured field. These results show that the developed multiscale formulation casts a promising framework for the real-field enhanced geothermal formations.