Accurate simulation of fluid flow and heat transfer in geothermal reservoirs is a crucial necessity for optimising energy extraction strategies. However, natural formations (including geothermal ones) extend large length scales (in the order of km), while their properties (e.g., heat and flow conductivity) can change in small (fine) scales (e.g., cm or even below). As such, accurate simulations of filed-scale models are too expensive to be handled by the state-of-the-art commercial simulators. As a matter of fact, these geological models are upscaled excessively in order to reduce the computational costs. Excessive upscaling leads to loss of accuracy and details of the heterogeneous properties, which can result in non-optimum production estimation and operation strategies. As a remedy, in this work, we propose a dynamic multilevel method (ADM) which captures small-scale heterogeneity (i.e., accurate) while preserving the computational efficiency (thus applicable to field-scale models). This development is achieved by combining two major concepts: (1) multiscale basis functions for accurate coarse-scale treatment of heat and flow conductive properties at their original fine-scale, and (2) adaptive mesh refinement strategy to minimise the requirement for employing the fine-scale grid, i.e., when and where needed. These two developments combined in one framework allows for both accurate and efficient simulation of coupled flow-heat equations in subsurface geothermal reservoirs. The fine-scale grid is employed only at the cold water front, where most of the nonlinear (and grid resolution sensitive) interaction is taking place. The rest of the domain is solved at the coarser scales, depending on the gradients (slope of change) of the unknowns. Note that no upscaling is needed due to the employment of multiscale basis functions. Moreover, these basis functions are calculated only once at the beginning and reused for the rest of the time-dependent simulations. Our method allows for non-thermal equilibrium between rock and injected fluid, so to allow for full flexibility and possible added accuracy. Through several test cases, the accuracy of the proposed ADM is investigated by measuring its error compared with the fine-scale fully resolved simulation. Its accuracy, on the other hand, is measured through calculating the average number of active grid cells. Our results, for both homogeneous and heterogeneous models, show that the proposed method employs a fraction of the fine-scale grids to deliver accurate solutions. Therefore, it provides a promising framework for field-scale simulation of geothermal reservoirs.