Advances in the modeling of the Iberian thermal lithosphere and perspectives on deep geothermal studies

Renewable energy sources are key to achieve the transition toward clean energy system. Among them, the geothermal energy has a production whose effectiveness requires sufficient understanding of the temperature distribution and fluid circulation at depth, as well as of the lithological and petrophysical properties of the crust. The focus of this paper is twofold: first, we summarize the main advances in the development of new methodologies and numerical codes to characterize the properties of the thermal lithosphere in terms of its, temperature, density and composition; second, based on the compilation of available thermal modelling results, we present the depth of the thermal Lithosphere–Asthenosphere Boundary (LAB) of the Iberian Peninsula and the temperature distribution at crustal depths of 5, 10, and 20 km, in addition to at Moho level. At 5 km depth, the temperature is above 110 °C with local anomalies (> 130 °C) located in the Iberian Massif and Cenozoic volcanic provinces. A similar pattern is observed at 10 and 20 km depth, where temperatures are above 190 °C and 350 °C, respectively. At 20 km depth, anomalies above > 500 °C, delineate the SE and NE Cenozoic volcanic provinces. At Moho depths, temperature ranges from 450 to 800 °C with hot regions mainly located along the Iberian Massif and the SE and NE volcanic provinces. The compiled results do not show any lithospheric anomaly that could give rise to high temperatures at shallow depths, but they do show an acceptable exploitation potential at intermediate depths. With regard to the direct use of district and greenhouse heating and for industrial processes, the potential is great throughout the Peninsula, the main challenges being the availability of groundwater and drilling costs.


Introduction
The successful mitigation of climate change and preservation of the environment rely on the transition toward clean and renewable energy sources.Geothermal energy can help achieve the sought-after decarbonization of energy system, thanks to its ability to decarbonize heating and its energy storage potential (Zhao and You 2020).Near surface geothermal energy is a widely used resource; however, in many cases, low temperatures (< 50 °C) restrict its use.On the other hand, in stable continental settings, high-temperature geothermal energy suitable for electricity production is often a deep difficultto-locate resource.Finding deep geothermal resources requires a good understanding of the regional temperature distribution and fluid circulation at depth, together with knowledge of the lithological and petrophysical properties of the lithospheric crust and mantle, hence the importance of modeling the thermal structure of the lithosphere.
The study of the lithospheric structure of the Iberian Peninsula (Fig. 1) and, in particular, its thermal properties, has attracted the attention of numerous authors who, from different disciplines and using different methodologies, have tried to interpret the geodynamic complexity of the area.
The large amount of seismic data acquired during the last decades on the Iberian Peninsula and its surrounding margins has provided us with relatively detailed knowledge of the crust.This information has improved thermal modeling of the Iberian Peninsula, since they have helped define the structure of the crust and the lithospheric mantle.
The focus of this work is twofold: the first objective is to summarize the main advances in the development of new methodologies and numerical codes that have characterized the properties of the thermal lithosphere in terms of its temperature, density, and composition.The second objective is to report the most noteworthy results from the thermal modeling of the Iberian lithosphere in terms of its temperature distribution at crustal depths and implications of its medium-high geothermal potential.
The study also raises open issues, such as: (i) how to best explain the differences observed between the thermal and the seismic LAB; (ii) obtaining a more even distribution of surface heat flux (SHF) measurements covering the entire Iberian Peninsula, particularly along the majority of the Iberian Massif, Iberian Chain, Ebro and Tagus foreland basins, and in those areas with potential geothermal sources; (iii) the need for laboratory measurements of thermal rock properties to better constrain thermal modeling parameters and prospection at the survey site; (iv) to develop detailed crustal thermal models, with higher lateral resolution and with thermal parameters that account for the lithologies of the modeled crustal rocks, and (v) to advance geodynamic and thermal modelling to better account for thermal transient effects, as it occurs in the southern areas of the Iberian Peninsula, and to a lesser degree at its northern boundary.
In this study, we have gathered available surface heat flow measurements (Fig. 3a) and thermal modeling results from: (i) 2D profiles from Fernandez et al. (2004) ;Carballo et al. (2015a and2015b); Palomeras et al. (2011); Pedreira et al. (2015) and Jiménez-Munt et al. (2019), (ii) 1D inversion of geoid and elevation plus 3D gravity modelling of Torne et al. (2015), and (iii) the 1D non-linear Bayesian inversion of Fullea et al. (2021) (Fig. 3b).In addition, we also discuss results of the thermal crustal structure deduced from magnetic data presented in Andrés et al. (2018).All thermal models presented in this work are based on three basic assumptions: a planar approximation (Cartesian coordinates), thermal isostasy and thermal steady state.The steady-state thermal assumption with no advective transport poses some limitations in young tectono-thermal regions, as in case for the northern and southern boundaries of Iberia (Fig. 1), where main tectonic activity ended recently-northern boundary-or it is still active-southern boundary- (Vergés and Fernàndez 2006).A transient thermal model could improve the resulting density and temperature distribution; however, the vertically averaged density would have to remain similar to reproduce the observed elevation.Quantifying the transient effects requires time evolving numerical modelling of geodynamic processes for which the timing and even the occurrence are poorly constrained.
In addition, the 2D lithospheric profiles and 1D Bayesian inversion model integrate geophysical and petrological data to determine the thermo-chemical structure of the crust and upper mantle down to 400 km (Afonso et al. 2008 andKumar et al. 2021).The crustal and upper mantle structure is constrained by simultaneously fitting surface heat flux (SHF), elevation, Bouguer anomaly, and geoid height.Since all of these observables depend on the thermophysical properties of the materials used, which in turn depend on temperature, pressure and composition, and furthermore, they have different sensitivities to thermal and/or density anomalies at different depths, the approach allows for improved control of density, temperature, and composition variations at different depths.
Moreover, available seismic data and tomographic models are integrated to constrain the geometry and densities of the crustal layers, the Moho depth and compare V p and V s distribution at mantle lithospheric levels.For a more detailed discussion on the method and the model parameters used, we refer the reader to Afonso et al. (2008), Kumar et al. (2020), andFullea et al. (2021).

Tectonic setting
The Iberian Peninsula, which once comprised the Iberian microplate, is located at the western end of the boundary of the Eurasian and African plates that extends from the Himalayas at its easternmost end, to the triple junction at the Azores, in the Central Atlantic (Fig. 1a).In the study area, the plate boundary is diffuse and has remained active since the Mesozoic, passing from a generalized extensional regime during the latest Paleozoic-Triassic to a compressional one that began in the Late Cretaceous because of the NW drift of the African plate relative to the Eurasian plate (e.g., Macchiavelli et al. 2017).
The lithospheric structure of the Iberian Peninsula and, specifically, its thermal structure, results from the superposition of three large orogenic cycles, and, in particular, of the two most recent, the Variscan Cycle, which spans from 480 to 250 Ma and the subsequent Alpine Cycle, which started in the Permo-Triassic and continues to the present day.These orogenic cycles have left their imprint on the current configuration of the Iberian Peninsula, where two large domains can be distinguished: the Variscan domain, mostly located in the west, and the Alpine domain, in the east, bordered by the Atlantic Ocean and the Western Mediterranean Neogene basins, respectively (Fig. 1b).
The Variscan terrains (pre-Permian rocks) form the basement of the Iberian Peninsula that crops out in the Iberian Massif and also in small areas of the Axial Zone of the Pyrenees, Betics, Iberian System, Catalan Coastal Range and on the island of Menorca.The Iberian Massif formed at the end of the Paleozoic by the collision of Laurasia and Gondwana.The European segment of the Variscan Cordillera, which is characterized by folded and thrusted structures and large-scale metamorphic processes, has been built by the accretion and collision of continental fragments bordering Gondwana during the subduction of the Rheic and PaleoTethys oceans, before the collision between the two large continental masses with Laurasia to the north and Gondwana to the south (Matte 2001;Martínez-Catalán, 2007).
The Iberian Massif, together with the Armorican Massif in France, form the Iberian-Armorican Arc, the main outcrop of the Variscan Cordillera in Western Europe.Depending on the position within the arc, six zones can be distinguished: the Cantabrian and South Portuguese zones, located to the north and south, respectively, show typical characteristics of external zones (abundant synorogenic sediments), while the West Asturian-Leonese, Central Iberian, and Ossa-Morena zones show typical internal zone characteristics, including very intense deformation, magmatism and metamorphism.In addition, two large suture zones can be distinguished: on the limits of the Ossa-Morena and South Portuguese zones and in the Galicia-Tras-Os-Montes Zone, which is an allochthonous terrain with ophiolites, formed as product of the closure of the Rheic Ocean (Arenas et al. 2004;Pérez-Estaún and Bea 2004;Martínez-Catalán, 2007, 2011;Azor et al. 2008, 2019and Simancas et al. 2013) (Fig. 1b).The Iberian Massif has remained practically stable over the last 300 Myr (Gibbons and Moreno 2002).
The continental part of the Alpine domain includes the Cantabrian-Pyrenees mountain belts and their associated foreland basins in the north, the Iberian Ranges and the Catalan Coastal Range in the northeast, the Central System in the central region, and in the south at the current plate boundary, the Gulf of Cádiz-Gibraltar Arc zone, with the Betic-Rif system and the Alboran back-arc basin (e.g., Vergés et al. 2019, among others).
The Alpine Cycle began with an initial extensional period in the Permian and Triassic that shaped the Iberian microplate.This first phase of extension is followed by a second one that spanned from the Jurassic to the middle Late Cretaceous (201-83.5 Ma), which gave rise to the opening of the Central Atlantic and the North Atlantic, and to the propagation of a rift system within and bordering the Iberian Peninsula.The rifting ended with the formation of oceanic crust in the Ligurian-Tethys, the opening of the Bay of Biscay (Vergés et al. 2019) and the displacement of Iberia toward the east (Ziegler 1999;Salas and Casas 1993).In the middle of the Late Cretaceous (~ 83.5 Ma), the increase in the spreading rate of the South Atlantic caused Africa to begin a counterclockwise rotation with respect to Eurasia at the same time that it began to move northward, initiating its long-lasting convergence against Europe with the concomitant closure of the Ligurian-Tethys Ocean.The Alpine orogeny continues to the present day, with the inversion of most of the extensional Mesozoic basins, giving rise to the current mountain ranges and the opening of the Mediterranean Neogene basins (Fig. 1b) (Quesada and Oliveira 2019).
In terms of volcanism, we can broadly distinguish Paleozoic and Cenozoic volcanism, with some manifestations in the Oligocene and Quaternary.The first one is generally located in the Iberian Massif and the second one in the eastern half of the Iberian Peninsula and in the Neogene basins of the Western Mediterranean.According to López- Ruiz et al. (2002), three large provinces can be distinguished on the Iberian Peninsula, which are from south to north: the South-East Volcanic Province (SEVP) encompassing the area of Cabo de Gata-Mazarrón-Cartagena, with a very heterogeneous and complex orogenic-type volcanism, which ranges from calc-alkaline, calc-alkaline with high potassium content, or even to shoshonite and ultrapotassic rocks; the Calatrava Volcanic Province (CVP); and the North-East Volcanic Province (NEVP) including the Empordà-Selva-La Garrotxa zone, both with intraplate alkaline volcanism (anorogenic).Offshore, in the Gulf of Valencia, two main volcanic periods stand out.The calc-alkaline orogenic volcanism of Oligocene-Miocene (24-18.6Ma) age and the alkaline type that extends from the Tortonian to the beginning of the Holocene (Marti et al. 1992;Melchiorre et al. 2017a;Marti and Bolós, 2019).

Advances in modeling the thermal lithosphere of Iberia
The exploration for shallow and deep geothermal resources of the Iberian Peninsula was initiated by the Geological and Mining Institute of Spain (IGME) in the 1970s, with the development of the General Inventory of Geothermal Manifestations, in which the first geothermal potential results were presented for the entire region.Subsequently, Albert-Beltran (1979a, b) presented an initial attempt to correlate geothermal anomalies and crustal thicknesses.Later on, Fernandez and Banda (1989) published the map of geothermal gradients and heat flux for the NE peninsular region first and subsequently for the entire Iberian Peninsula (Banda et al. 1991;Fernandez et al. 1998).
Integrated thermal modelling begins with the work of Fernàndez et al. (1990) and Zeyen and Fernandez (1994) in which, for the first time a methodology is presented that, using finite elements, allows 2D modeling of the thermal structure and density distribution in the lithosphere through joint modeling of elevation data, gravimetry, and heat flow (CAGES code).
The initial CAGES algorithm developed by Fernàndez et al. (1990) was improved by incorporating the geoid anomaly, which allowed for better control of the density distribution at the lithospheric mantle scale.The algorithm was used to study the lithosphere in different geodynamic scenarios, as is the case of the Atlantic margin (Torne et al. 1995), the Cantabrian margin (Ayarza et al. 2004), the transition from the Iberian Massif to the oceanic crust of the Gulf of Cádiz (Fernàndez et al. 2004), the NW margin of Morocco and the Atlas mountains (Zeyen et al. 2005;Jiménez-Munt et al. 2010), the Western Mediterranean (Roca et al.2004), or the Variscan terrain of the SW of the Iberian Peninsula (Palomeras et al. 2011), among other studies.Subsequently, geochemical and petrological data of the mantle were incorporated (Afonso et al. 2008) and the possibility of introducing thermal, seismic, or compositional anomalies, or a combination of them into the sublithospheric mantle (Kumar et al. 2020), has improved the potential of these algorithms (LitMod2D code) considerably.For examples applied to the study of the Iberian Peninsula, we refer to the works of Carballo et al. (2015a, b), Pedreira et al. (2015), Jiménez-Munt et al. (2019), and Kumar et al. (2021).
In parallel, 3D forward thermal lithospheric modeling algorithms, GEO3Dmod (Fullea et al. 2007), and the 1D inversion of elevation and geoid height were developed to obtain the thermal structure of the lithosphere (Fullea et al. 2007).Subsequently, the LitMod3D code was developed, which integrates the joint direct modeling of geophysical and petrological data of the lithosphere and sublithospheric mantle in 3D (Fullea et al. 2009) and the LitMod4I code, which is dedicated to the non-linear Bayesian geophysical-petrological inversion of surface waves, surface heat flux, elevation, and geoid anomalies (Afonso et al. 2013a, b).The works of Fullea et al. (2007Fullea et al. ( , 2010) ) and Torne et al. (2015) are examples of the application of these algorithms to the Iberian Peninsula.In the study of Fullea et al. (2007) and Torne et al. (2015), geoid and elevation inversion are integrated with 3D gravimetric modeling and inversion, while the work of Fullea et al. (2010) applies the LitMod3D algorithm to study the structure of the lithosphere-asthenosphere system in the Atlantic-Mediterranean transition zone for the first time.More recently, Andrés et al. (2018) propose a thermal structure of Iberia and its margins from magnetic data, assuming the 580 °C isotherm as the Curie temperature or Curie limit.Fullea et al. (2021), based on geophysical integrated Bayesian inversion and petrological studies address the origin of the topography of the Iberian Peninsula using the thermal and compositional structure of the lithosphere.
To these results we must add the studies that, through different approaches, deduce the thermal structure of the lithosphere of the central zone of Iberia (Tejero and Ruiz, 2002), and southern margin (Torne et al. 2000;Soto et al. 2008), among others.

The thermal structure of the lithosphere in the Iberian Peninsula
In this section, we summarize the main characteristics of the crust and lithospheric mantle that can directly or indirectly affect the distribution of temperatures at depth and the thermal regime that currently prevails in Iberia.To this aim, we have compiled the available results from 2D profiles and 3D regional thermal models located in Fig. 3b.

Thickness and thermal structure of the Iberian crust
In the Moho map of Fig. 4, we observe that there is a clear correlation between crustal thickness and elevation for both the continental and oceanic domains most importantly with regard to medium and long wavelengths (from tens of kilometers) (Figs. 3  and 4).This correlation is an expression of the principle of isostasy, hence the relief of the Alpine chains correspond to crustal thicknesses greater than 40 km, while locally they can reach values greater than 55 km (e.g., the Pyrenees and the Cantabrian Mountains).Values between 32 and 36 km are obtained in the Cenozoic Ebro, Duero, Tajo, and Guadalquivir basins, while we observe a relatively constant thickness in the Iberian Massif, ranging between 30 to 34 km, except in the SW sector, where there is some crustal thinning in an NE-SW direction that continues toward the Gulf of Cádiz.Offshore, the crustal thickness ranges from 22 to 26 km on the continental shelf, where it is characterized by a relatively shallow bathymetry, between 0 and 200 m, and 10-18 km on the abyssal plains with deeper bathymetries (> 2500 m depth).For a thorough discussion on Iberian Moho details we refer the reader to Díaz et al. (2016 and2021).
It is worth highlighting two regions that diverge from the crustal thickness-elevation correlation that is to be expected according to the principle of crustal isostasy.In the first place, the crustal thinning of the SW of the Iberian Peninsula (crustal thicknesses of 26 to 30 km), is characterized by positive regional anomalies of both gravimetric and geoid height (e.g., Fernàndez et al. 2004), coinciding with an average topographic height of about 200 m (Vergés and Fernàndez 2006).This crustal thinning, could partially explain the presence of both maxima in the absence of an appreciable topographic depression in the area.Second, we have found little or no crustal expression in some reliefs of the Iberian Massif (Sierra Morena or Montes de Toledo) and the Catalan Coastal Range (Figs. 3 and 4).Studies of the elastic thickness of the lithosphere (T e ) in different areas of the Peninsula (e.g., García-Castellanos et al. 2002;Jiménez-Díaz et al. 2012) obtain values between 10 and 30 km, which explain that part of these reliefs are supported by the rigidity of the lithosphere and, therefore, with no (or, little) expression at the base of the crust.Unlike the onshore areas, the oceanic domain is characterized by significant crustal thinning (between 12 and 15 km), similar in absolute terms but with a different geometry and extension, than that observed throughout the Iberian Atlantic Margin, Cantabrian Margin, and the Western Mediterranean.
The compilation of seismic data and the results of the thermal modeling of the Iberian Peninsula (Fig. 5) shows that, in general, both the Variscan and Alpine crust are structured in three layers: the upper and middle crust of variable thickness with velocities between 5.4-6.2 and 6.2-6.5 km/s, respectively, and the lower crust with velocities between 6.5 and 7.2 km/s, velocities that are slightly higher on average (6.9-7.0 km/s) in the Variscan crust (Díaz and Gallart 2009).Rabbel et al. (2013) suggest that the increase in V p may be related to a change in the composition of the lower Variscan crust with respect to the Alpine crust, although a structural origin cannot be ruled out.In addition, the lower crust is generally more reflective than the upper and middle in the Variscan and Alpine domains, (Ayarza et al. 2021), with the exception of the Western Mediterranean basins, where this reflective character is lost in the areas of greater amounts of thinning (Watts et al. 1990, Torne et al. 1992;Collier et al. 1994).The reflectivity of the lower crust is either related to lithological stratification produced by intrusions of mafic rocks with compositional stratification within the intrusion itself (Condie 2015) or to the presence of shear zones (e.g., Reston 1988;Clerc et al. 2015, among others).
Thermal modelling results show that this stratification is also reflected in variations of the values of radiogenic heat production and thermal conductivity obtained from the best fit model that are summarized in Fig. 5.In this figure, we observe how, along the Pyrenean orogen, the conductivity is slightly higher in the middle/lower crust of the Alpine zone (3.1 and 2.5 W/mK) than in the Variscan zone (2.1 and 2.0 W/mK), while for the rest of the Peninsula, the conductivity values remain constant (2.4,2.1, and 2.0 W/m K), except for the SW Iberian Massif, which shows slightly higher values (Fig. 5).On the contrary, in the northern margin, radiogenic heat production is lower in the Alpine upper crust than in the Variscan one (1.0 and 1.65 µW/m 3 ), while this trend is reversed for the lower crust (0.3 and 0.2 µW/m 3 , respectively).We also highlight the differences observed at the level of the middle/lower crust between the northern half of the Iberian Peninsula (Duero Basin) and its southern half (Tajo Basin) with values of radiogenic heat production for the middle/lower crust being slightly higher in the northern part (1.0 and 2.4 vs 0.5 and 0.2 µW/m 3 , respectively).With regard to the upper crust, a nearly constant value is obtained throughout the area (1.65 µW/m 3 ), except for the Axial Zone of the Pyrenees and the SW Iberian Massif, which we comment on below.
The stratification of the crust is locally disturbed by the presence of high-velocity layers/bodies, as in the case of the Iberian Reflective Body (IRB) that Simancas et al. (2003)  relate to an intrusion of a mafic rock sill (Fig. 5) with velocities of 7.0 km/s and radiogenic heat production and thermal conductivity of 0.3 µW/m 3 and 2.5 W/m .K, respectively.Another example is found in the Cantabrian Mountains with the presence of an eclogitized lower crust penetrating the lithospheric mantle, with velocities between 6.9 and 7.2 km/s, as well as a heat production and thermal conductivity of 0.33 µW/m 3 and 2.0 W/m .K, respectively (Carballo et al. 2015a;Pedreira et al. 2015).This eclogitized lower crust is probably present throughout the Pyrenees, as suggested by the ECORS profile (Daignières et al. 1981).Furthermore, on the northern and southern margins of Iberia and along the Western Iberian Atlantic Margin, mantle rocks are found at different crustal levels, emplaced during the Mesozoic extensional events.On the southern margin (Betics and Rif ), mantle rocks emplacement is related to their subduction under the African margin and their subsequent retreat toward the Iberian margin, which caused these rocks to rise along the subduction channel, being emplaced at different crustal levels and even outcropping at the surface, e.g., the Ronda Massif (Betics) or Beni-Bousera (Rif ) peridotites (Melchiorre et al. 2017b).
Some caution should be taken with the thermal parameters obtained from thermal modelling.Typically, thermal conductivity and radiogenic heat production are inferred from modeling by best-fitting the observables, in particular surface heat flux.However, obtaining reliable data on surface heat flux is challenging.Thus, the results listed above should be taken as proxies until more precise lab measurements are available.

The thermal LAB
Unlike the Iberian crust, which is relatively well-understood based on multiple studies, the base of the lithosphere (LAB) and, in particular, the structure and composition of the lithospheric mantle, have historically received less attention.The first data reporting on the lithospheric mantle of Iberia was obtained at the end of the 1980s, through the Iberian Lithosphere Heterogeneity and Anisotropy project (ILIHA), in which arrivals were detected from around 60 km depth (Díaz et al. 1993).Subsequently, the thermal modeling of different geophysical observables, either in 1D, 2D or 3D, has allowed us to obtain the depth of the thermal LAB and propose variations in the composition of the lithospheric mantle, which were discussed in Fullea et al. (2021).To these results, we must add the seismic tomography studies that provides a detailed image of the distribution of velocity anomalies in the upper mantle of Iberia, which have been used in the thermal models presented in this study to further constrain the density distribution of the lithospheric mantle.In addition, lateral variations (that is, at constant pressure) of the seismic velocity are fundamentally related to variations in the thermal structure of the lithosphere and asthenosphere.The most relevant tomographic studies used in the thermal modeling presented in this study are surface wave seismic data from Palomeras et al. (2017), global and regional travel-time and multifrequency tomography models from Bezada et al. (2013), Bonnin et al. (2014), Civiero et al. (2018), andVillaseñor et al. (2015).For the Western Mediterranean, the S-wave full waveform inversion of Fitchtner and Villaseñor (2015), and for the Pyrenees, the teleseismic P-to-S converted model of Macquet et al. (2014) and Chevrot et al. (2015) in addition to the compilation of passive seismic data from Chevrot et al. (2022).
Figure 4a, d shows the depth of the thermal LAB in the Iberian Peninsula.We observe that, although there are differences regarding the absolute values, which are mainly related to the different thermal parameters used in the thermal modeling, all the results coincide with the lithospheric thickening from the Variscan terranes of Western Iberia to the Alpine chains.Maximum LAB depths, which range between 200 and 170 km depending on the modelling used, are located in the north and northeast of Iberia, throughout the Pyrenees, in the Cantabrian Mountains, and in the Iberian Range, with local maxima of 200 and 190 km located in the Axial Zone of the Pyrenees and Cantabrian Mountains, respectively.On the contrary, the western and southwestern zone of Iberia is associated with significant thinning that extends both toward, the west in the Gulf of Cádiz and the Iberia Atlantic Margin, with values of 100 km and toward the Western Mediterranean Neogene basins, where values lower than 70 km are found.This trend changes completely for the Gibraltar Arc, which shows an important lithospheric thickening with values greater than 190 km.
Broadly speaking, we see that in the interior of Iberia the Alpine lithosphere is thicker than the Variscan lithosphere.In the oceanic part, the thinnest lithosphere is located in areas that have experienced recent extensional episodes, as is the case of the back-arc basins of the Western Mediterranean (Valencia Trough, Algerian Basin, and the eastern sector of the Alboran Basin).For this regional trend, the minimum located at the SW of the Iberian Peninsula, in the Ossa-Morena and South Portuguese Zones, and the maximum associated with the Gibraltar Arc clearly stand out.The minimum of the SW zone largely derives from the observed high in the geoid height (7 m) located in a partially elevated area with an average topography of 385 m (Casas-Sainz and de Vicente 2009) that implies a density deficit in the lithospheric mantle.According to Fernàndez et al. (2004), this deficit can be interpreted either as due to a thinning of the lithospheric mantle or to a decrease in its average density on the order of 25 kg/m 3 .The first hypothesis is consistent with the presence of a thermal anomaly related to the geodynamic processes that affect the entire Gibraltar Arc, while the second one may be related to a depletion of the mantle during the Variscan orogeny.
Instead, the lithospheric thickening that extends from the Western Betics to North Africa, with maximum values sof 250 km in the Rif area (e.g., Fullea et al. 2010 andJimenez-Munt et al. 2019) contrasts with the significant thinning observed in the eastern part of the Alboran basin (70 km).This lateral variation of lithospheric thickness, occurring across a distance of about 300 km, reflects the great complexity that is observed at all levels along the present boundary of the Eurasian and African plates.The slow (currently about 3-4 mm/year) and prolonged convergence between Africa and Eurasia has resulted in a sharp boundary in the oceanic domain that is limited by the E-W trending Azores-Gibraltar Fault Zone that extends from the triple point of the Azores to the Gorringe Bank.This is in contrast with the continental domain, which is characterized by a diffuse limit, where transpressive deformation in an NW-SE direction predominates.This regional field is superimposed on the opening of the Alboran back-arc basin, which was as active from 24 to 7 Ma (Late Oligocene-Early Miocene to Late Tortonian) behind the Gibraltar Arc migration (e.g., Rosenbaum et al. 2002;Faccenna et al. 2004;Jolivet et al. 2009;Fullea et al. 2010;Vergés and Fernàndez 2012;van Hinsbergen et al. 2014;among others).
Differences in the thermal LAB depth determined by the models discussed above amount to values of ± 40 km on average, with local anomalies that diverge in their location and intensity.Comparing the results from Fullea et al. (2021) and Torne et al. (2015), we observe that the LAB of Fullea and co-authors is shallower than that proposed by Torne and co-authors for the majority of the Iberian Peninsula, with the exception of the Pyrenees and SW region (Fig. 4a, d).Maximum differences are observed along the North Iberian Massif and in the SE and NE Volcanic Provinces, with local anomalies above 50 km.In general, 2D profiles show a deeper LAB, which best match the results of Torne et al. (2015) with the exception of the Gibraltar strait region and the Western Mediterranean basins.In these areas, 2D profiles show a much deeper (> 190 km) and shallower (< 70 km) LAB, respectively, than that proposed by Fullea et al. (2021) and Torne et al. (2015) (Fig. 4a, d).

Temperature distribution at depth
One of the advantages of the thermal modeling presented in the previous sections is that, in addition to the details of the structure and distribution of densities in the lithosphere and upper mantle, it allows obtaining the temperature distribution at different depths with greater or lesser detail, depending on whether the estimation is based on 2D or 3D models.The greater detail of the 2D models is implicit in the modeling itself, as the 2D models allow for more refinement in most of the defined thermo-chemical parameters.6,7,8,9 show and compare the temperature distribution at 5, 10, and 20 km and at Moho depth; these were obtained from the regional and 2D thermal models and also from Andres et al. (2018) at Moho depth from Curie Depth Point modelling.
At crustal depths of 5 km, the regional models show temperatures in the range of 75-150 °C, although they differ in their location of the maxima.In the model of Fullea et al. (2021) two regional maxima are observed in the central part of the Iberian Massif and in the Campo de Calatrava and SE Volcanic Provinces in the Almeria, Murcia, and Levante fields, while in the model of Torne et al. (2015) these maxima are associated, in most cases, with the Alpine topographic reliefs (Fig. 6).(2015a) vs the regional model of Fullea et al. (2021).n: number of analyzed points.r = correlation coefficient At 10 km depth, basal temperatures are above 150 °C, being slightly higher (75-100 °C) in the model of Fullea et al. (2021) (Figs. 8b and 9a).Also noteworthy are the maximum > 250 °C of Fullea et al. (2021) located in the central regions of the Iberian Massif and the volcanic zone of Murcia (Cartagena-Mar Menor), and the high values (> 250 °C) obtained from 2D profiles in the southern Betics-Alboran region (Fig. 6).In the Alboran region, the temperature increase is supported by the measured high SHF values; however, it is not consistent with the low values measured in the Betics region (< 60 mW/m 2 ) (Fig. 3a).Differences between the regional models are mostly in the range of 75-100 °C (Fig. 9b).2D profiles show similar results to those obtained by Torné et al. (2015), with the exception of the Betics-Alboran region, and mostly lower temperatures when compared to the model of Fullea et al. (2021) (Fig. 9b).
At 20 km depth, the model of Torne et al. (2015) shows values above 350 °C that coincide in the northern half with results from 2D modelling, whereas Fullea et al. (2021) show values slightly higher, above 400 °C, with high intensity local temperature anomalies (> 500 °C) of the SE and NE volcanic provinces (Fig. 7).In the Betics and Alboran regions, 2D models show temperatures above 550 °C, much higher than those obtained by the regional models.
The comparison of the crustal temperature at 5 and 10 km depth obtained by the regional models (Fig. 9a) shows that the model of Fullea et al. (2021) predicts higher temperatures than those obtained by Torne et al. (2015), for the interval of 0 to 100 °C.At 20 km depth, the observed differences increase to 100-150 °C.The correlation coefficients between both models are rather poor, particularly at 10 and 20 km depth, − 0.254 and 0.098, respectively, which contrasts with the correlation coefficient of 0.915 obtained at Moho depth (Fig. 9a).
Comparing results from the 2D N-S profile of Carballo et al. (2015a) (Profile 4 of Figs. 6 and 7) with the results obtained from the model of Torne et al. (2015), we observe a good correlation from 10 km down to the Moho, with correlation coefficient values ranging from 0.616 at 10 km depth to 0.771 at Moho depth.
On the contrary, when comparing the results from the model of Fullea et al. (2021) to the N-S 2D profile of Carballo et al. (2015a) (Fig. 3b), the obtained correlation coefficients are below 0.2, with the exception of the value at Moho depth, where it increases to 0.72.The model of Fullea et al. ( 2021) also shows higher temperatures (in the range of 0 to 100 °C) than those obtained from the N-S 2D profile (Fig. 9c).
The temperature at the Moho derived from the Curie-Point Depth (Andres et al. 2018) is higher than that obtained by the thermal models.Figure 7 shows that the hottest areas (> 850 °C) are located along the Pyrenees-Cantabrian system and the NE corner of the Iberian Chain.Lower temperature values (< 750 °C) are found along the Iberian Massif with the exception of its SW most corner.

Discussion
In this study, we have summarized the regional trend of temperature distribution at crustal depths in the Iberian Peninsula, based on the results of the thermal modelling performed so far.Above 10 km depth, there is coincidence in the temperature distribution of thermal anomalies with some differences in their intensities and exact location.The analyzed models coincide in the sense that, at 5 km depth crustal temperatures range from 75 to 150 °C with local maxima located in the Iberian Massif and SE and NE Volcanic Provinces.A similar pattern is observed at 10 km depth, where temperatures range from 200 to 275 °C, with local anomalies up to 290 °C.
Results from 2D models mainly differ from 3D regional models in the Betics and Alboran Sea, where they register temperatures above 275 °C.These high temperatures are likely related to the lithospheric thinning obtained in the transition from the Iberian Margin to the Alboran Basin, which is not so well-constrained in the regional models.The temperature anomalies of the Iberian Massif and SE Volcanic Province proposed by Fullea et al. (2021) mainly trace the outcrops of the Variscan granitoids of the Iberian Massif and the volcanic rocks of the Calatrava and SE Volcanic Provinces (Figs. 1b  and 7).In the volcanic provinces they are supported by the measured surface heat flux (Fig. 3a).
Some differences are seen in the temperature distribution at Moho depths between the 3D regional thermal models.The model of Fullea et al. (2021), which is based on the Bayesian geophysical-petrological inversion of surface waves, surface heat flux, elevation, and geoid anomalies obtains higher temperatures (75-100 °C on average), than those deduced by the model of Torne et al. (2015), which is based on geoid and elevation inversion integrated with 3D inversion gravimetric modeling.Temperatures based on the heat flow deduced from CPD (Curie-Point Depth) (Andrés et al. 2018) show that there is a clear difference in the temperature at the Moho between the Variscan and Alpine terranes, the latter ones being much hot as 100 to 150 °C (Fig. 8).
The differences observed between the model of Fullea et al. (2021) and the model of Torne et al. (2015) and 2D profiles are related to the different information they use and their lateral resolution, which results in differences in the topography of the Moho and LAB depths, e.g., at the western regions.These differences are reflected in the calculated correlation coefficients (Fig. 9), with the lowest values at shallow depths, which increase with depth to values greater than 0.7 at Moho depth.
In summary, the model of Torne et al. (2015) and the 2D profiles vary in a similar way (correlation coefficient above 0.6), except at upper crustal levels (5 km depth) (Fig. 9).This is mainly related to the enhanced definition of the thermal parameters at upper crustal levels in the 2D profiles.The results from Fullea et al. (2021) show greater differences compared to those of Torne et al. (2015) and the 2D profiles, with higher temperatures at the depths analyzed (Fig. 9).
The differences and similarities observed between the analyzed thermal models are largely controlled by the depth of the thermal LAB, the assumed isotherm, and to a lesser extent, by the values of the thermal parameters used.2D thermal models permit a more detailed definition of the structure of the crust and the lithospheric mantle, and in turn, of their thermophysical properties.Fullea et al. (2007) assess the effect of varying the thermal parameters within geologically meaningful ranges and conclude that the lithospheric thickness, and hence the lithospheric geothermal gradient, is mostly affected by changes in the linear thermal expansion coefficient in the mantle, and to a lesser extent by changes in crustal thermal conductivity and radiogenic heat production.Changes in thermal conductivity and radiogenic heat production are more relevant at upper crustal levels, since they have a greater influence in the shallow geotherm.
The densities inferred from inversion of geoid and elevation methodology (Torne et al. 2015) and from the 2D approaches, are highly dependent on the chosen reference column.In the 2D models, the assumed composition of the mantle, which indirectly will control the temperature distribution in the lithosphere, is also essential.Although the approach of Fullea et al ( 2021) is based on similar assumptions, the inclusion of surfacewave data, which are more sensitive to the lithospheric geotherm is an additional constraint.This explains the smoother topography of the Moho and LAB geometry obtained by these authors compared to the results from 2D and potential field and elevation inversion, in addition to the different lateral resolution of the models discussed herein.
All of the presented models assume a thermal steady state, and although this assumption is valid for ancient stable regions, as could be the case for the Iberian Massif, it could be at odds in recently deformed Alpine regions, which are likely affected by transient perturbations in the temperature distribution.In these regions, deformed during the Cenozoic, the thermal steady state assumption may overestimate the actual lithospheric thickness in thinned tectonic domains, while it may underestimate it in cases with significant lithospheric thickening.Thus, the temperature distribution presented in this work should be taken with some caution, particularly in the Gibraltar Arc segment of the Betic-Rif orogenic system, which is presently being deformed by ongoing NW-SE transpressive tectonic stresses.
The compilation of thermal conductivity and radiogenic heat production at crustal levels, resulting from the best fit model, shows that on northern margin of the Iberian Peninsula the conductivity is slightly higher in the middle-lower crust of the Alpine domain (3.1 and 2.5 W/mK) than in it is the Variscan zone (2.1 and 2.0 W/mK), while for the rest of the Iberian Peninsula, conductivity values remain constant (2.4,2.1, and 2.0 W/mK).On the contrary, radiogenic heat production is lower in the upper Alpine crust than in the Variscan crust (1.0 and 1.65 µW/m 3 ), while this trend is reversed for the lower crust (0.3 and 0.2 µW/m 3 , respectively).There are also some differences in the middlelower crust between the northern half of the Iberian Peninsula (Duero Basin) and its southern half (Tajo Basin) with values of 1.0 and 2.4 vs 0.5 and 0.2 µW/m 3 , respectively.With regard to the upper crust, a practically constant value is obtained throughout the study area (1.65 µW/m 3 ), except for the Axial Zone of the Pyrenees and the SW Iberian Massif.We raise some caution with these results, since the thermo-physical parameters obtained from thermal modelling should be taken as proxies that result from the best fit model, but they cannot replace in-situ or lab measurements.Site-specific values of crustal radiogenic heat production and thermal conductivity are essential to constrain the shallow geotherm and, hence, to evaluate regions of elevated geothermal gradients with the potential for geothermal energy exploitation.
Development in thermal modeling that integrates different geophysical observables and petrological data has made it possible to study and propose temperature, composition, and density distributions in the Iberian lithosphere.However, much effort is needed to improve the knowledge of the petrophysical properties of the crustal rocks to better assess the thermo-physical parameters of the different crustal domains.
Unlike the crust, the thermal and compositional structure of the lithospheric and sublithospheric mantle of the Iberia Peninsula remain less well-known.The general pattern of lithospheric thickness, a parameter that exerts considerable influence over the crustal geotherm, is relatively well-known from geophysical models constrained mainly by seismological and gravity data (Fullea et al. 2007(Fullea et al. , 2010(Fullea et al. , 2021;;Torne et al. 2015;Carballo et al. 2015a, b;Pedreira et al. 2015;Jiménez-Munt et al. 2019;Kumar et al. 2021).The present-day consensus is that in the interior of Iberia the Alpine lithosphere is thicker than the western Variscan lithosphere, with thinned lithosphere also present in the SE Mediterranean Iberian margin (the SE Volcanic Province).
Understanding the thermal lithosphere in detail has practical implications for the evaluation of the geothermal potential of medium-high enthalpy sites.The first conclusion drawn from the models of the Iberian Peninsula is that there is no evidence of lithospheric anomalies that can give rise to high temperatures at shallow depths, but they do show good potential at intermediate depths.Taking the available temperature data from existing oil wells in Spain as a reference (IGME 1987), we consider less than 2000 m as shallow depth, and around 5000 m as high depth (Fig. 10).From the graph, we can deduce that temperatures above 120 °C at depths between 3000 and 5000 m can be frequent in sedimentary basins, reaching an exceptionally high value of 200 °C at 4750 m on the Alboran margin, likely related to the volcanic activity in this region.
Therefore, with regard to the direct use of geothermal energy for district and greenhouse heating, and for industrial processes, the potential seems to be excellent throughout the Peninsula, the main challenges being the availability of groundwater and the drilling costs.On the other hand, as for the efficient production of electricity with conventional geothermal energy, at least 150 °C is usually required, and they could be found in different areas of the Peninsula at medium-great depths.Of course, the thermal gradient can significantly vary locally, and thus, for a better evaluation of the resources, it is necessary to refine the models.Considering the main geothermal play systems for electricity production, both hot dry rock (HDR) and hot sedimentary aquifer (HSA) models would be located on the Iberian Peninsula.HDR systems would be found mainly to the west, in the Variscan zone, associated with deep well-rooted granites with a sedimentary cover acting as a heat trap.On the other hand, the HSA systems with the greatest potential would be found in the east of the peninsula associated with thinned areas of the European Rift and the Mediterranean margin.Special attention will be required in the Alboran back-arc basin, which experiences extremely complex geodynamics and with evidence of geothermal anomalies.

Concluding remarks
In this study, we have examined the regional temperature distribution trend at crustal depths in the Iberia Peninsula, based on the results of the thermal modelling performed so far.Integrated geophysical-petrological thermal modeling allowed us to study and propose a first approximation of the temperature distribution, composition, and density distribution in the Iberian Peninsula.There is consensus that the Alpine lithosphere in the east of Iberia, is thicker than the Variscan lithosphere, in the west.Thinned lithosphere is also present in the SE Mediterranean Iberian Margin (the SE Volcanic Province).
Lithospheric thermal models also show no evidence of lithospheric anomalies giving rise to high temperatures at shallow depths, but they do show a suitable geothermal potential at intermediate depths.
The potential for direct use of geothermal energy for district and greenhouse heating, and industrial processes, seems to be excellent throughout the Iberian Peninsula, the main challenges being the availability of groundwater and drilling costs.
Regarding electricity production, hot dry rock (HDR) systems would be located mainly in the Iberian Massif, while hot sedimentary aquifer (HSA) systems, with the greatest potential, would be located in the eastern areas with a thin lithosphere and on the Mediterranean Margin of Iberia.
Site-specific values of crustal radiogenic heat production and thermal conductivity as well as an even distribution of surface heat flux measurements are essential to better constrain the shallow geotherm on the Iberian Peninsula and, hence, to evaluate regions of elevated geothermal gradients with the potential for geothermal energy exploitation.
The Iberian Massif and the Mediterranean Margin show the highest potential geothermal resources and should be, therefore, prioritized in future resource estimation investigations.

Fig. 4
Fig. 4 Results of the 2D lithospheric-thermal models superimposed on the thermal model of Fullea et al. (2021) (top a-c) and of Torne et al. (2015) (bottom d-f).a, d Base of the thermal lithosphere, b, e Moho depth and c, f comparison of the Moho obtained by thermal modeling and that observed from seismic Deep Seismic Sounding (triangles) and Receiver Functions (squares).Seismic results from Díaz et al. (2016).Thick colored lines show results from 2D thermal models referenced in Fig. 3 caption and along the text

Fig. 5
Fig. 5 Crustal profiles showing the distribution of V p , thermal conductivity and radiogenic heat production in the crust and uppermost lithospheric mantle.The map in the lower right corner shows the location of the profiles.Díaz and Gallart (2009) Adapted from, Pedreira et al. (2007), Carballo et al. (2015a, b) and Palomeraset al. (2011)

Figures
Figures 6, 7, 8, 9  show and compare the temperature distribution at 5, 10, and 20 km and at Moho depth; these were obtained from the regional and 2D thermal models and also fromAndres et al. (2018)  at Moho depth from Curie Depth Point modelling.At crustal depths of 5 km, the regional models show temperatures in the range of 75-150 °C, although they differ in their location of the maxima.In the model of Fullea et al. (2021) two regional maxima are observed in the central part of the Iberian Massif and in the Campo de Calatrava and SE Volcanic Provinces in the Almeria, Murcia, and Levante fields, while in the model ofTorne et al. (2015) these maxima are associated, in most cases, with the Alpine topographic reliefs (Fig.6).

Fig. 9
Fig.9Temperature differences and correlation coefficient between the regional thermal models ofFullea et al. (2021) andTorne et al. (2015) and the N-S 2D profile ofCarballo et al. (2015a) at 5 km, 10 km, and 20 km and at Moho depth.See Fig. 3b for location and text for details.Comparison of: (a) the regional model of Fullea et al. (2021) vs the model of Torne et al. (2015); (b) comparison of the N-S 2D profile of Carballo et al. (2015a) vs the regional model of Torne et al. (2015); (c) comparison of the N-S 2D profile of Carballo et al. (2015a) vs the regional model ofFullea et al. (2021).n: number of analyzed points.r = correlation coefficient