Carbonate rocks represent the most abundant geothermal aquifers worldwide, excluding active volcanic areas. The most important challenge is to accurately define the conceptual model of a regional hydrothermal system. The focus of modelling is generally on the outflow area of the system and its fluid flow dynamics, while the geological and hydrogeological characterisations of the recharge area and the estimation of the local infiltration rate are rough. Thermal water is usually characterised by long residence time, and classical hydrogeological balance can give misleading results. On the other hand, long-term numerical simulations could represent a useful tool to understand the regional and local fluid flow circulations and to estimate the hydrogeological features of the recharge area. Numerical modelling has to be supported by an accurate characterisation of the regional and local hydraulic and thermal properties of the reservoir. This methodology was applied for the first time in Croatia in the case study of Daruvar hydrothermal system. Its recharge area had been determined decades ago solely by geological mapping, as a wide zone in the mountainous hinterland of natural thermal springs. Mesozoic fractured carbonate rocks represent the geothermal aquifer, and the tectonic juxtaposition of these permeable formations with low permeable Neogene deposits enables the rising of thermal water (38-50 degrees C) up to the surface. 2D coupled flow and heat transport numerical simulations of the hydrothermal system were performed for the first time in this research using measured data on its temperatures, and hydraulic and thermal conductivities. However, their results did not support the existing conceptual model. Conceptual model of the system was therefore revised decreasing the extent of the possible infiltration area and considering the fracturing induced by local faults in the Daruvar subsurface. A new set of simulations was conducted using the proposed model. The obtained results reproduced the temperature of the Daruvar thermal springs and the regional fluid flow. Despite its limitations, the employed modelling approach was useful to perform a first simulation of the hydrothermal system dynamics, and it could be employed in similar systems to evaluate the consistency of different hypotheses on the available conceptual model and their impact on the temperature distributions. In addition, the modelling demonstrated that the recharge area is smaller than previously considered, which has implications on the protection of the whole hydrothermal system. The greatest danger to the system currently comes from the increasing anthropic impact in the recharge area by multiple active dolomite quarries in the region and the related seismic disturbances, which should be limited in the newly defined recharge area.