Geothermal energy is one of the renewable energy sources foreseen in the European Union’s plans for the green energy transition. Its sustainable utilisation mostly depends on the characteristics of the geothermal resource from which it is extracted. Among others, detailed geological modelling and reconstruction of the subsurface is a key factor for estimating the potential of a geothermal resource. In particular, it is crucial to determine the geometry of faults and fractures since the associated highly permeable damage zones represent a preferential pathway for the circulation of thermal fluids and their uprising (BENSE et al., 2013; FAULDS et al., 2013). This research focuses on the modelling and reconstruction of the geological and structural settings in the Daruvar thermal spring area using surface geophysical techniques. Electrical resistivity tomography (ERT) was employed to delineate the geometry of resistivity layers in the subsurface. The ERT results were combined with an integrated approach based on the passive Horizontal to Vertical Spectral Ratio (HVSR) method and the active Multichannel Analysis of Surface Waves (MASW) method to map the thickness of the Quaternary cover. The geophysical data were constrained using the stratigraphic logs of wells in the spring area and its surroundings obtaining a 3D reconstruction of the geological setting. The spatial distribution of resistivity shows relatively low values from 10 to 150 Ωm (Fig. 1). Based on ERT results and the stratigraphic logs of the wells, three resistivity layers/geological units were identified (Fig. 1): (1) the Quaternary alluvial cover with resistivity ranging between 30 and 50 Ωm (layer 1); (2) the Neogene sediments with resistivity values of 10–30 Ωm (layer 2); and (3) the Triassic dolomites that were divided in a compact layer with resistivity ranging from 70 to 150 Ωm (layer 3a) and a fractured layer characterised by resistivity values of 20–30 Ωm (layer 3b). Furthermore, sharp lateral variations in the resistivity distributions were observed. They were generally marked by low resistivity anomalies that were interpreted as the fracture zones along the faults (F1 to F5; Fig. 1). The high secondary porosity of the fault zones and the occurrence of thermal waters decrease the bulk resistivity of the rock mass. The obtained results allowed us to reconstruct the geological setting of the Daruvar thermal spring area. The Quaternary cover has a 5 to 15 m thickness, increasing northward and eastward. Its thickness and geometry were confirmed by seismic investigations. Neogene deposits are generally found below the alluvial cover, except for the central part of the study area where the Carbonate complex, i.e., Triassic dolomites are found. Local scale faults and their fracture zones enhance the upwelling of thermal waters resulting in the occurrence of thermal springs with temperatures up to 48 °C (BOROVIĆ et al., 2019). Two main faults border southward and eastward of the thermal spring area juxtaposing the highly permeable Triassic dolomites with the low permeable Neogene deposits. This lateral contrast fosters the rising of the thermal waters forming a shallow thermal resource in the Daruvar area.