We aim at modeling the main crustal and thermal interfaces of Sicily (Italy), a key area for understanding the geological complexity at the collisional boundary between the African and European plates. To this end, we analyze the gravity and magnetic fields, integrated with information from well logs, geology, heat flow, and seismic data. In order to make the most accurate description of the crustal structure of the area, we modeled with different methodologies the carbonate and crystalline top surfaces, as well as the Moho and the Curie isotherm surface. The reconstruction of the carbonate platform is achieved using a nonlinear 3D method constrained by the available seismic and borehole data. The crystalline top, the Curie, and the Moho are instead estimated by spectral analysis of both gravity and magnetic data. The results show a complex carbonate basement and a deep crystalline crust in central Sicily, with a prominent uplift beneath the Hyblean Plateau. Maps of the Moho and the Curie isotherm surface define a variable thermal and structural setting of Sicily, with very thin crust in the southern and eastern sectors, where high heat flow is found, and deep and cold crust below the Caltanissetta Basin.