Unravelling the lithospheric-scale thermal field of the North Patagonian Massif plateau (Argentina) and its relations to the topographic evolution of the area

The North Patagonian Massif (NPM) area in Argentina includes a plateau of 1200 m a.s.l. (meters above sea level) average height, which is 500–700 m higher than its surrounding areas. The plateau shows no evidence of internal deformation, while the surrounding basins have been deformed during Cenozoic orogenic events. Previous works suggested that the plateau formation was caused by a lithospheric uplift event during the Paleogene. However, the causative processes responsible for the plateau origin and its current state remain speculative. To address some of these questions, we carried out 3D lithospheric-scale steady-state and transient thermal simulations of the NPM and its surroundings, as based on an existing 3D geological model of the area. Our results are indicative of a thicker and warmer lithosphere below the NPM plateau compared with its surroundings, suggesting that the plateau is still isostatically buoyant and thus explaining its present-day elevation. The transient thermal simulations agree with a heating event in the mantle during the Paleogene as the causative process leading to lithospheric uplift in the region and indicate that the thermo-mechanical effects of such an event would still be influencing the plateau evolution today. Although the elevation related to the heating would not be enough to reach the present plateau topography, we discuss other mechanisms, also connected with the mantle heating, that may have caused the observed relief. Lithosphere cooling in the plateau is ongoing, being delayed by the presence of a thick crust enriched in radiogenic minerals as compared to its sides, resulting in a thermal configuration that has yet to reach thermodynamic equilibrium.


Introduction
The North Patagonian Massif (NPM) is located in the backarc of the Southern Andes in central Argentina. It includes a plateau of an average altitude of 1200 m a.s.l. (meters above sea level; Coira et al. 1975) that is between 500 to 700 m higher in altitude than its surrounding basins. Specifically, the plateau's boundaries can be traced following four lineaments, the Gastre, Limay, Gualicho and Los Chacays (each of them comprising a similar length of ~ 320 km and restricted to the plateau area), with its western border being located approximately 400 km east to the subduction trench (Fig. 1). The NPM plateau occupies an area of ~ 100,000 km 2 and its rocks appear mostly undeformed . It is composed of Paleozoic crystalline basement and mid-Jurassic volcanic rocks terminating upward in a regional planation surface (Aragón et al. 2010;Gonzalez Díaz and Malagnino 1984) covered by shallow marine sedimentary rocks of Danian-Maastrichtian age. This sequence, in turn, is partially overlain by Oligocene to early Miocene mafic volcanic fields, the basalts of Somún Curá (Kay et al. 2007), that cover the southern portion of the plateau but also erupted over the adjacent region.
The present-day elevation of the plateau, as well as the presence of the Somún-Curá basalts, are difficult to explain in the geodynamic context in which they were formed. Geological observations and previous works (Aragón et al. 2010 and references therein) indicated that the current elevation of the plateau was established in response to a lithospheric uplift event during the Paleogene. This conclusion was based on two observation: (1) that the marine sedimentary strata of the Cretaceous/Paleogene boundary are found both on top of the plateau and in the surrounding basins, indicating that at times of their deposition, the area of the present plateau must have been below sea level; and (2) that Somún Curá basalts are thought to have spilled downwards from the highest elevation of the plateau to the surrounding areas, which suggests that the area was already elevated at the time of their eruption in the late Oligocene. The lithospheric uplift in the Paleogene occurred under regional horizontal extension in northern Patagonia, predating the late Miocene Andean horizontal shortening stage (Giacosa and Heredia 2004;Giacosa et al. 2005;Aragón et al. 2011;Fernández Paz et al. 2018;Iannelli et al. 2017;Iannelli et al. 2018). Horizontal shortening during the Andean orogenic phase was responsible for the build-up of the southern Andes and led to extensive deformation of the surrounding areas around the plateau. However, the marine Cretaceous/Paleogene sedimentary rocks within the plateau remained undeformed . In summary, based on these observations, we could interpret lithosphere uplift in the plateau area to date back to the Paleogene within a horizontal extensional regime, where the plateau was uplifted as an epeirogenic block without significant internal deformation. Moreover, horizontal shortening during Andean orogeny in the late Miocene deformed the surroundings and the plateau block differentially , likely indicating the presence of a rheologically differentiated lithosphere beneath the plateau compared to the surrounding Patagonian crust.
The tectonic configuration of the area during the Paleogene uplift was rather complex. Cande and Leslie (1986) and Somoza and Ghidella (2005) discussed how late Cretaceous subduction of the Farallón/Alluk ridge beneath the South American plate resulted in a parallel motion between the Farallon and the South American plates, thus  Amante and Eakins 2009). The location of the modelled area is delimited in white, the NPM plateau in red and the outcropping basalts in grey. The black stippled line denotes the limits between provinces. Province names are written in green, basins in blue and lineaments in white (modified from Gomez Dacal et al. 2017) hindering subduction in northern Patagonia during most of the Paleogene. According to Aragón et al. (2011), this new plate configuration would have led to the opening of a slab window of the "fraternal type" (Thorkelson 1996), which would have resulted in regional mantle upwelling and horizontal extension in the area of northern Patagonia. This is expressed in the syn-extensional intra plate volcanism caused by decompressional melting of the upwelling mantle and in the development of the Patagonides rift (Aragón et al. 2018). According to this geodynamic model, the NPM plateau, being surrounded by the Patagonides rift, would have been uplifted as part of an extensional regime, during a plate re-organization caused by the subduction of the Farallon/Aluk active ridge. Other studies have proposed alternative hypotheses to reconcile the onset of a thermal anomaly in the mantle, which also led to the extrusion of the Somún Curá basalts (De Ignacio et al. 2001;Kay et al. 1993Kay et al. , 2007Muñoz et al. 2000). Kay et al. (1993Kay et al. ( , 2007 considered the presence of a plume-like upwelling that intersected the subducting plate right below the NPM in the Oligocene. Other studies related mantle upwelling either to a roll back of the subducting plate (Muñoz et al. 2000) or to corner flow as a response to the plate geometry (De Ignacio et al. 2001).
The exact configuration of the lithosphere beneath the plateau is still a matter of debate, given the paucity of geophysical data available. Recent efforts to collect and integrate all available information have been summarized by Gomez Dacal et al. (2017), who developed a 3D lithospheric-scale geological model of the area, consistent with gravity data and constrained by other independent observations including seismological data. Their 3D model is by far the most complete reconstruction of the lithosphere of the area in terms of present-day geometry, density and rock lithology. One main feature of the 3D model is the presence of a thicker crust below the NPM plateau as compared to its surroundings, which suggests the plateau to be isostatically more buoyant. Gomez Dacal et al. (2017) also showed that a low-density anomaly in the upper mantle below the southern portion of the NPM would explain the observed gravity field and that this low density spatially coincides with a low velocity body that extends further to the south as constrained by the available tomographic data (Amaru 2007;Schaeffer and Lebedev 2013). This, in turn, would provide an additional argument for the presence of a buoyant mantle beneath the plateau. Moreover, the low velocity anomaly in the mantle, observed by both S and P waves, is indicative of high temperatures, which, in combination with a thickened crust, would also point to a buoyant lithosphere beneath the plateau.
All of these observations raise the question whether a causal relationship exists between the presence of a buoyant lithosphere and the plateau elevation. More specifically phrased, we ask whether it is possible to explain the presentday topography of the area and obtain clues to its evolution by unravelling variations of the temperature distribution in space and time, also taking into account the heterogeneities of the 3D lithosphere configuration. To address this question, we carry out a numerical study to shed light on: (i) the present-day 3D thermal configuration of the plateau and its surroundings and (ii) its evolution through time after the onset of a mantle heating event during the Paleogene, considered as possible causative mechanism for the plateau formation. In our study, we consider a heterogeneous lithosphere as constrained by the 3D geological model of Gomez Dacal et al. (2017) that undergoes conductive cooling in response to a (quasi)instantaneous thermal anomaly. We parameterize the latter by raising the topography of the 1300° C isotherm, representative of the thermal Lithosphere Asthenosphere Boundary (LAB), as a proxy for the proposed mantle dynamics. With this general model setup given, we have systematically varied the magnitude of the thermal perturbation imposed, to analyze the temporal and spatial conduction of the thermal signal through the lithosphere. This allows us to investigate the amounts of relative vertical movement and speculate about the present-day thermal stability of the lithosphere in the area.

Input 3D geological model
As an input to the thermal modelling, we use the gravityconstrained 3D structural and density model of Gomez Dacal et al. (2017). This model is the final result of a data integration procedure conducted to constrain the lithospheric configuration of the NPM plateau area and its surroundings. The model includes geological and geophysical information of different types (seismological, gravimetric, geological observations, among others) and of different extents, from local data to regional or global models Chernicoff and Zappettini 2004;Heine 2007;Kostadinoff et al. 2005;Pankhurst et al. 2006;Pavlis et al. 2012;Ramos 1988, among others). All gravity-independent information was integrated in a 3D density model, which was then used as input to assess the depth configuration of the Moho by considering the observed gravity field based on the general density contrast between the crust and the mantle (Gomez Dacal et al. 2017). The final model has limitations related to the scarcity of relevant data in the area; yet it is the most complete model of the studied domain at this scale to date. The most relevant characteristics of the model for the present study are summarized in Fig. 2a, b and c (geometries) and in Table 1 (lithologies and densities). In what follows, we briefly describe the pertinent features of the model, while 1 3 the interested reader is referred to the original publication by Gomez Dacal et al. (2017) for more details.
The model resolves four layers (a sedimentary cover, two different crustal domains and the lithospheric mantle), each of them limited by the following surfaces: topography, top to the crystalline basement and Moho (Figs. 1 and 2). It extends laterally for 500 × 500 km; thereby, including the NPM plateau and its surrounding areas. The topography of the NPM plateau ( Fig. 1) was derived from the digital elevation model (DEM) ETOPO 1 of 1 arc minute resolution (Amante and Eakins 2009), that portrays not only the high elevation of the plateau but also its little relief variations inside. The depth to the top of the crystalline basement was derived from the Mobil isopach data of the ICONS Atlas (Heine 2007), a global compilation of sediment thickness of intra-continental sedimentary basins. Accordingly, the only domain where the overlying sedimentary cover exceeds 200 m of thickness is the sector of the Neuquén basin ( Fig. 2a). For their 3D gravity modelling, Gomez Dacal et al. (2017) assigned a constant density of 2400 kg/ m 3 to this sedimentary cover considering the main composition of the different lithologies found in Neuquén basin and the minor effects of compaction on the effective porosity. This density value was tested in a sensitivity analysis which finally showed that such near-surface variations would not change the main conclusions regarding the deep lithosphere density configuration. The depth to the Moho (Fig. 2c) is another result of the density and gravity modelling as detailed in Gomez Dacal et al. (2017). After integrating the available gravity-independent constraining data in the wider region, the authors found a large negative gravity residual (reaching -100 mGal) just in the plateau area. Considering typically large density differences between the crust and the mantle (> 200 kg/m 3 ) and the seismically constrained Moho depths north of the plateau, the depth of the Moho beneath the plateau was well  Dacal et al. 2017). Boundary between crustal domains is shown in green. c Depth to the Moho (Gomez Dacal et al. 2017). d Depth to the thermal LAB constrained by the observed gravity field. Accordingly, the depth to the Moho below the plateau surface varies between 40 and 50 km; thereby, deepening from the plateau boundaries to its center (Fig. 2c). The steepest Moho gradient is observed below the northwestern limit of the NPM plateau, i.e. along the border with the Neuquén basin. The smallest Moho depth gradient, is featured along the western border where the plateau Moho is bounded by the crustal root of the Andes (Gomez Dacal et al. 2017). As a main result of the modelling process of Gomez Dacal et al. 2017, the NPM plateau reveals a larger crustal thickness than its surroundings (Fig. 2b).
Information on the density of the deeper crust was derived from the S-wave tomographic model SL2013Sv (Schaeffer and Lebedev 2013), which is a global model that assures data coverage in every knot of the grid (lateral resolution of 0.5° and vertical of 25 km). Seismic velocities of the tomographic model were converted into densities using Birch's law (Birch 1961), which was slightly modified to fit the densities calculated from lower crustal xenoliths (granulites) found in Paleogene basalts of Paso de los Indios . The derived crustal density distribution depicts a lateral differentiation in two domains that largely coincide with the location of the Paleozoic terranes Chilenia and Patagonia, proposed by Ramos (1988) and their sutures as also derived from magnetic data (Chernicoff and Zappettini 2004). Using these evidences in combination with surface geological observations (Gregori et al. 2008;Kostadinoff et al. 2005;Pankhurst et al. 2006). Gomez Dacal et al. (2017) laterally divided the crystalline crust into two domains of average density: Chilenia (2700 kg/m 3 ) in the north and Patagonia (2810 kg/m 3 ; Fig. 2b; Table 1) in the south. Finally, the model includes a mantle domain that features spatially variable densities, as obtained by a careful investigation and combination of geophysical and petrological data (Gomez Dacal et al. 2017). The data used were derived from two tomographic models: SL2013Sv by Schaeffer and Lebedev (2013) and the global P-wave velocity model UU-P07 by Amaru (2007). In spite of not having a good distribution of local seismic observation points in the area, both models assure to have enough coverage of rays crossing each cell of the grid inside the selected zone. Gomez Dacal et al. (2017) obtained different distributions of density for the mantle by converting both velocity models (Vp: Amaru 2007 and Vs: Schaeffer and Lebedev 2013), using different algorithms (An and Shi 2007;Cammarano et al. 2003;Goes et al. 2000;Karato and Karki 2001). To select the most appropriate distribution, the authors compared them with the density calculated from more than 50 samples of mantle xenoliths distributed in 4 localities inside the plateau and one to the south (Mundl et al. 2015;Ponce et al. 2015). The xenoliths from inside the plateau are mostly spinel-hazburgites and spinel or garnet-peridotites that were emplaced in alkali basalts during Miocene to Pleistocene times (Mundl et al. 2015). The xenoliths found to the south of the plateau are spinel peridotites, pyroxenites, hazburgites and lherzolites found in Eocene basalts (Ponce et al. 2015). The results from the comparison was conclusive and indicated that the best consistency between the seismological and the petrological data is obtained for the case of a conversion from P-wave tomographic model of Amaru (2007), using the software Velt (An and Shi 2007) that follows the approach of Goes et al. (2000). Therefore, this distribution of densities for the mantle was incorporated in the 3D gravity modelling (Gomez Dacal et al. 2017). The model finally adopted in this study is supported by different types of data, which confirms its liabilityalthough some parts of the study area are still poorly constrained. Therefore, even though the main result of the model, the crustal thickness, has a certain degree of uncertainty, it is also the case that its spatial variations are supported by the results of a systematic sensitivity analysis of the respective gravity responses of different input parameters, particularly by well analyzed distribution of densities in the mantle (Gomez Dacal et al. 2017). To summarize the main outcome of Gomez Dacal et al. (2017) is that the observed negative gravity anomaly beneath the plateau can only be caused by a lighter and thicker crust and neither by any observed mantle density anomaly, nor by the presence of the relative thin sedimentary cover.

Thermal modelling approach
To unravel the transient evolution of the thermal configuration of the lithosphere since Paleogene times, we solve the differential equation of conductive heat transfer by a 3D Finite Elements Method. For this we rely on the open-source software GOLEM (Cacace and Jacquey 2017; Jacquey and Cacace 2017), a numerical simulator for Thermal-Hydraulic-Mechanical modelling. Under the assumption of a conductive lithosphere, as conduction is the predominant heat transport mechanism in the lithosphere (Turcotte and Schubert 2002), the relevant equation reads as: where k is the rock thermal conductivity (considered as a bulk rock property in the following analysis), S is the rate of heat production by radioactive decay and the term on the right hand side accounts for thermal energy storage, with DT Dt being the total Lagrangian derivative, the rock density and c p the isobaric heat capacity. In this study, we solve Eq. (1) for two different purposes: (1) to conduct a steady-state conductive thermal model to simulate the present-day distribution of temperatures, therefore without accounting for the variation of the temperature through time DT Dt = 0 and (2) to evaluate how temperatures vary through time (and space) using a transient thermal modelling approach. Equation 1 is an initial and boundary value problem (IBVP), the solution of which depends on (1) the selection of the relevant thermal properties and (2) a proper choice of boundary and initial conditions.

Thermal properties
Rock thermal properties have been assigned to each unit of the 3D model, as based on their prevalent lithology Table 1). In doing so, we have used resolved compositional changes derived from the structural and density model of Gomez Dacal et al. (2017). For the sediments, the value for the radiogenic heat production was selected according to the specific composition (average value for sandstones, claystones and siltstones; Table 1) while for the thermal conductivity, we considered the influence of the porosity to compute a bulk material property as: where k b is the bulk thermal conductivity which is used in the temperature calculation, k w the thermal conductivity of the water assumed equal to 0.6 W/(m °K), k s the thermal conductivity of the rock matrix. The latter is calculated as the average of the values for sand, clay and siltstones as found in Midttomme and Roaldstet (1999). Finally, the sediment porosity ∅ was derived from the gravity constrained bulk rock density ( b ) using the following expression: and considering typical values for the matrix density ( s ) and the water density ( w ).
The thermal properties for the crystalline crust in the Patagonia crustal domain were considered as an average between (1) mean values for metamorphic rocks (Čermák and Rybach 1982;Vilà et al. 2010) as mainly found in outcrops of the domain and (2) granulites (Čermák and Rybach 1982;Vilà et al. 2010), corresponding to the composition of the lower crustal xenoliths of the area . The values for the Chilenia domain were calculated as the mean between (1) the values for a plutonic composition (Čermák and Rybach 1982;Vilà et al. 2010), the main lithology of the domain obtained from gravity constrained density modelling (Gomez Dacal et al. 2017) and (2) granulites, thereby assuming that the lower crustal composition of the two terrains is similar. For the Patagonia terrane, we also took into account chemical analyses performed on lower crust xenoliths ) and outcropping rocks (Pankhurst et al. 2006). These chemical data were used to obtain an additional independent estimate of the radiogenic heat production. Accordingly, we define the concentrations of the radioactive isotopes of Potassium, Uranium and Thorium (C K, C U and C Th, respectively), which were then used together with the specific radiogenic heat production of each element (S K , S U and S Th respectively; Turcotte and Schubert 2002) to obtain a final bulk value as: The values obtained by the two different approaches for the radiogenic heat production of the crystalline crust of Patagonia were consistent and the value selected for the thermal modelling is shown in Table 1. The mantle values for k and S were selected as the ones assigned by Čermák and Rybach (1982) and Vilà et al. (2010) respectively, for an ultramafic composition.
As based on the above described analyses, there is a small range of possible values for both the thermal conductivity and radiogenic heat production that can be assumed for each layer (Table 1). Therefore, we tested the sensitivity of the overall thermal modelling results to these ranges. The final values chosen (Table 1) were selected as the ones that best reproduce the few available well temperature measurements (courtesy of the company YPF; personal communication; Fig. 2a). We count on five measured bottom-hole temperatures, of which three are located to the North of the plateau and the other two to the South (Fig. 2a). All measurements range in depth between 730 and 1400 m and most of them (4 out of 5) had been taken from the crystalline basement. The final model (Table 1) predicts local temperatures that all deviate from the spatially corresponding measured temperatures by less than 4.5 °C, while the median of temperature differences is -0.168 °C and the standard deviation is 3.3 °C. The other performed test gave differences between measured and calculated temperatures that are higher (> 15 °C). An analysis of the impact of the selected values for the most relevant thermal properties is described in the final discussion to the manuscript.

Boundary conditions
We have chosen to use fixed temperatures for both the top and the bottom boundary of the model (Dirichlet boundary conditions), while considering all lateral boundaries as closed for heat flow. We assigned a value of 9.5 °C at the topography, which corresponds to the mean annual temperature in the area according to Del Barrio and Martin (2012). For the lower boundary condition, we considered the depth to the 1300 °C isotherm, as obtained from the conversion of P wave velocities (Amaru et al. 2007) to temperatures by means of the software Velt (An and Shi 2007). Thus, we follow the common approach that this isotherm represents the thermal LAB. We selected this specific tomographic model and conversion approach because the resulting density distribution is the one that shows the best fit with xenoliths data according to Gomez Dacal et al. (2017) (see "3D steady-state temperature distribution" for more details). It is worth mentioning that the conversion of P waves to temperatures suffers from a certain degree of uncertainty (of maximum ~ 100-150 °C), with the error increasing with increasing depth (Goes et al. 2000). However, in our study we limit to relatively shallow mantle depth levels (80% of the model has a maximum depth less than 150 km), thus minimizing the error from the conversion. Moreover, the temperatures used in this study are robust as Gomez Dacal et al. (2017) tested various conversions of different tomographic data obtaining a similar temperature trend and selected the final distribution as the one that best fitted independent petrological data. The depth to the 1300 °C isotherm varies between 110 and 210 km and increases towards the northwestern border of the model (Fig. 2d). We gridded the LAB surface with a 50 km spacing to keep the tomographic model resolution.
This LAB topography was set as a fixed boundary condition to compute the present-day thermal state of the study area under the assumption of steady-state heat conduction. The model provides a temperature distribution for a heterogeneous lithosphere being in thermal equilibrium with respect to the boundary conditions set. Therefore, these results can be used as a reference model to properly quantify the present-day degree of relaxation of the mantle thermal anomaly that we consider to have perturbed the system in the past and thus caused the plateau formation.

Initial conditions
With the transient analysis, we aim to test the degree of thermal relaxation due to a heating at mantle level during the Paleogene, as suggested by different authors De Ignacio et al. 2001;Kay et al. 1993Kay et al. , 2007Muñoz et al. 2000). For that purpose, we model such a thermal pulse in the Paleogene by an instantaneous increase in the temperature at the level of the present-day LAB and test the purely conductive response of the lithosphere to it, in terms of the evolution of the temperatures in time and space.
A basic assumption here is that the geological configuration has not changed significantly since Paleogene times, as there is no available data to constrain a possible small change in geometry and lithology. With this assumption, we also imply that the Paleogene volcanism did not affect the geometry distribution in the area, as there is no evidence of interaction with the crust in the magmatic rocks of this time (Kay et al. 2007). The latter statement is not surprising in a context of a regional horizontal tectonic extension, as the one that prevailed in the Paleogene in the study area. Since we know that our steady-state model is well constrained and parameterized, it is a good reference model to test the hypothesis of the heating in the mantle.
We designed the transient thermal simulations to cover two major phases: a first heating and a subsequent cooling phase (Fig. 3). In the heating phase, we tested a number of simulations with different setups in terms of magnitude of the imposed thermal anomaly and time duration of the perturbation. We gradually increased the temperature at the thermal LAB by 100, 300 and 500 °C, over different periods of time (from 20 to 40 Ma), using as initial condition the temperature distribution of the steadystate thermal model (Fig. 3). We have chosen to apply an increase of temperature of 300 °C, which is an average value commonly assigned for a thermal increment produced by a mantle plume in the thermal LAB (Farnetani 1997;Herzberg et al. 2007, among others). As we cannot exclude that also other mantle heating mechanisms have been active, we also tested possible extreme variations of this temperature increase (100° and 500 °C). The time duration of the thermal anomaly was also modelled following end member scenarios, as it is not yet clear when the heating would have started in nature. Kay et al. (2007) suggested that the thermal anomaly already existed in the late Eocene because of evidence of active volcanism at that time in the area. For this specific setting, we took a mid-Eocene value at 45 Ma ago as the starting time for the heat anomaly. Aragón et al. (2015) argued; however, that the anomaly should have developed earlier because there is evidence of partial melting in the mantle as registered in individual eruptions dating back to 52 Ma before present. For this case, we have set the start of the heating in the model earlier, i.e. at 65 Ma.
Surface subsidence within the surrounding areas to the plateau concomitant with a marine transgression at around 20 Ma ago can be taken as indicative of the end of surface expressions of the mantle heat anomaly . We therefore consider 25 Ma ago as the ending time for the deep mantle heating in our simulations, just after the eruption of the Somún Curá basalts.
During the second phase (the cooling phase), lasting from 25 Ma ago to the present-day, the system is let to relax the thermal anomaly and thus we imposed a value of 1300 °C to the level of the present-day thermal LAB (Fig. 3).

Thermal isostasy
We also compute the topographic response to modelled variations in the thermal field, under the assumption of instantaneous (thermal) isostasy, following the approach of Hasterok and Chapman (2007): where Δ T is the elevation change predicted for a regional geotherm T(z) compared with a reference geotherm T ref (z) , V is the coefficient of thermal expansion and Z max is the maximum depth of the region we are evaluating.
We consider the steady-state temperature as the reference and evaluated the difference with the geotherms in the same point of space but in different times. We integrated them over the full lithosphere and used two different coefficient of thermal expansion, one for the crust ( 3 × 10 −5 K −1 ) and other for the mantle ( 3.2 × 10 −5 K −1 ), as suggested by Hasterok and Chapman (2007). Figure 4 shows the depth maps (below surface) of the 100 °C and the 450 °C isotherms. The depth to the 100 °C isotherm (Fig. 4a) varies from 2.5 to 3.3 km below surface, thereby being considerably shallower below the NPM plateau than below its surroundings. The maximum depth difference of the isotherm (of ~ 0.8 km) between the plateau and its adjacent domains can be observed in the NW limit of the plateau, whereas the smallest is predicted along the southern border (~ 0.1 km). For the largest parts of Patagonia, the depth of the 100 °C is shallower than 3 km below surface. In contrast, for Chilenia a greater variation in the isotherm depth is predicted, being shallower in the western area than in the eastern parts (Fig. 4a).

3D steady-state temperature distribution
By comparing the depth map of the 450 °C isotherm (Fig. 4b) with that of the 100 °C, we notice that the deep thermal field follows a similar pattern than the shallow one, although variations in temperature are larger in amplitude and wavelength. The 450 °C isotherm varies by more than 6 km in depth, rising from ~ 19.5 km below surface in the north-western corner of the model to ~ 13.5 km below surface in the middle of the NPM plateau (Fig. 4b). This again indicates a warmer NPM plateau as compared to the surrounding domains. The largest depth difference of the 450 °C isotherm (~ 6 km) is modelled between the plateau and the highest portion of the Andes included in the model. In contrast, inside Patagonia terrane the depth difference between the plateau and its neighboring domains is notably smaller (~ 1 km).
The presented 3D thermal model allows to directly relate modelled thermal anomalies with the structural configuration and thereof to the imposed lithological variations. Domains characterized by a shallow 100 °C isotherm, located close to the northern limit of the modelled area in Fig. 4a, coincide spatially with the Neuquén basin, the only area with sediments of significant thicknesses (> 200 m; Fig. 2a). Below the NPM plateau, where the crust is thicker than in its surroundings (Fig. 2b), the isotherms are shallower (Figs. 4a and b), thus indicating higher thermal gradients. The deep thermal field (represented by the 450 °C isotherm; Fig. 4b) can also be correlated with the subdivision of the crystalline crust into different terranes (Fig. 2b), superposed with variations in the LAB topography (Fig. 2d), which was used as lower boundary condition. By considering both depth maps, we can conclude that Patagonia is warmer than Chilenia, especially in the domain of the elevated NPM plateau (with topographies of > 500 masl; Fig. 4) than the surrounding domains.
To better understand the depth variations in the thermal configuration, in Fig. 5b we plot the temperature distribution along a profile that crosses the model area from NW to SE, together with its structural configuration (Fig. 5a). The difference in temperature between the plateau and its adjacent domain at depths between 1 and 2 km below sea level reaches maxima of approximately 80 °C. At the Moho depth, temperature differences across the profile are up to 350 °C (Fig. 5). Figure 6a shows the temperature evolution extracted at a point located in the centre of the NPM plateau at the depth of the Moho (test or observation point: longitude = − 68.5°, latitude = − 41°, depth = − 48.5 km) for two simulations. Both simulations are characterized by an increment of 300 °C in LAB temperature, while they differ in terms of the temporal duration of the heating phase: model A, 20 Ma and model B, 40 Ma. In both cases, we observe an adiabatic configuration during the early phase in the system evolution (~ 7 Ma), which is related to the time that it takes for the heat disturbance to propagate from the LAB up to the depth of the Moho. This stage is then followed by a phase in which the temperature increases in a mostly linear fashion as a function of elapsed time. The observed increase in temperature is mainly due to the thermal anomaly imposed at the lower boundary, having reached the Moho. This can explain why the two simulations show a similar temporal evolution, being both characterized by the same magnitude in the thermal anomaly. Heating at the Moho level continues even after the cessation of the thermal disturbance and it takes approximately another 7 Ma for the system to start cooling down, as seen by a gradual flattening of the temperature curves. The maximum temperature reached by the system is a nonlinear function of (1) conductive cooling after the end of the heating pulse and (2) the internal heat generation by radioactive decay within the crust. A maximum in temperature marks the onset of a stage in the system evolution when the contribution of internal crustal heat generation cannot keep up with thermal relaxation, as evidenced by a change in the slope of both temperature curves. This point marks a third stage in the evolution of the system where cooling prevails. However, despite cooling rates decreasing with time, within the time frame considered both simulations portray a lithosphere that has yet to reach thermal equilibrium, meaning that the thermal disturbance has not been fully relaxed by conductive cooling till the present day.

Transient thermal field
Having considered a different time duration for the heating phase translates into the two systems having reached different magnitudes in temperature values, with the amplitude of the thermal anomaly at the Moho positively correlating with the duration of the heating phase. By imposing a heating phase lasting for 40 Ma, we arrived at a temperature increment (difference in temperature at this point of the Moho between the initial steady state and the transient simulation at present-day) of about 60 °C, while by considering half of the duration the modelled change in temperature is less than 40 °C. We can also observe that the temperature gradients (indicated by the slopes in the curve during both the heating and the cooling phase) are higher (that is heat is conducted faster) in the case of a heating phase lasting over 40 Ma. Despite these differences, both simulations depict a system that was still heating up at the time when the Somún Curá basalt erupted (27 Ma) and that, at present day, has not yet reached a thermal equilibrium. The model considering a heating phase of 40 Ma is also able to match the time of the last volcanic eruption in the area (approximately 9 Ma ago), which coincides with the time where the system reaches a maximum in the computed temperature evolution.
In a second stage, we also carried out a sensitivity analysis of the modelling results with respect to the imposed value of the thermal anomaly, tested in a range varying from 100 °C up to 500 °C. The models were compared based on the maximum decrease in temperature between consecutive time steps during the cooling phase (Fig. 6b). All model realizations portray a similar trend in the temperature evolution, even though they each have different absolute values. Despite these differences, all simulations depict systems that have not reached a final equilibrium at present day, in agreement with the main findings from the two reference model described above.
In an attempt to quantify the effects of the thickness contrast modelled at crustal level between the plateau domain and its surrounding areas, we also compared the results from the model realizations extracted at two points in the crust (at 25 km depth), one being inside and the other one outside the NPM plateau. Figure 6c portrays the temperature evolution of these two points for 300 °C thermal anomaly along a time frame of 30 Ma. It can be observed that crustal rocks are less sensitive to the imposed anomaly inside the plateau, as  Fig. 4). a Main structural units and interfaces; b modelled steady-state temperature distribution for the present day indicated by a decrease in thermal gradients and a delay of approximately 15 Ma for the peak temperature with respect to the outside domain. The observed delay in the thermal response also implies that the crust beneath the plateau is further from its thermodynamic equilibrium at present day than the surrounding domain, which implies that the plateau is still subjected to thermal buoyancy forces that could explain its present state and high elevation.
To better quantify this last observation, in Fig. 7 we plot thermally induced changes in the topographic elevation (thermal uplift and subsidence) for the test point inside the NPM plateau. The simulation consisted of a thermal anomaly of 300 °C imposed along a heating phase lasting for 30 Ma. As expected, we notice uplift during the heating phase followed by gradual subsidence during the cooling phase, though the latter does not match the amount of induced thermal uplift. From our simplistic analytical model, we obtain a maximum of approximately 375 m at around 25 Ma ago, while the residual topography at present day is around 150 m.

Present-day state: temperature distribution and topography
Based on our 3D steady-state thermal modelling, we predict a distinct thermal configuration of the lithosphere inside the NPM Plateau as compared to the surrounding domains, as the lithosphere beneath the plateau is warmer (Fig. 4a and  b). We correlate this trend to the (petro)physical configuration of the lithosphere and in particular, to a thicker than normal crust beneath the plateau, which contributes to the heat budget of the plate by internal heat production through radiogenic decay. Temperature differences between the plateau and the surrounding domains increase with depth and reach their maximum at deeper crustal and mantle levels, where additional effects start to overprint the thermal field induced by variations in the LAB topology (Fig. 2d).
An interesting result of the model is that the resolved lateral heterogeneities in the crust, adopted to reflect the current subdivision into different Paleozoic terranes, do not exert any significant influence on the thermal field. We tested this influence by conducting a model with and homogeneous crystalline crust using average thermal properties. A description of this model and corresponding figures can be found in the Supplementary Material. The intra-crustal heterogeneities of the original model can be tracked in the geometry of the depth to the 450 °C isotherm (Fig. 4b), which results in a generally colder Chilenia, compared to Patagonia. This variation in the 450 °C isotherm is caused by the different thermal properties assigned following the different average lithologies for the two domains. Accordingly, Chilenia is more felsic than Patagonia and thus characterized by a higher radiogenic heat production and higher thermal conductivity (Table 1). Nevertheless, we also found a similar trend in the temperature pattern for the model assuming a homogeneous crust (Fig. S1), thereby demonstrating that the crustal heterogeneities have a secondary impact on the thermal field instead of the variations in crustal thickness.
In the whole model area, the thickness of the sedimentary cover is very small and the sedimentary rocks are mainly concentrated in the Neuquén basin (Fig. 2a), thus influencing the shallow distribution of temperatures only north of the NPM plateau (Fig. 4a). As the sediments have a lower thermal conductivity than the crystalline crust, their blanketing effect explains the temperature configuration at shallow depth.
From Fig. 6a, we derive a maximum temperature at the Moho beneath the NPM plateau of approximately 1005 °C at present day, which we relate to the presence of a thicker crust. This value appears to be higher than the average temperature expected for the Moho but we lack observational evidence that can either confirm or contradict our estimate.
In an attempt to "validate" how robust our estimates are, we also performed a sensitivity analysis by changing the main thermal rock properties within reasonable ranges (Table 1). Targeted properties tested comprise: (1) the mantle thermal conductivity that, in the final model, has been chosen to represent an upper value in the range that can be found in the literature (Table 1). The imposed thermal conductivity contrast at the Moho level results in an increase in heat being accumulated at this level, due to an enhanced conductive transport of thermal energy from the deeper mantle. Therefore, we tested the lower extreme value of thermal conductivity for the mantle given its lithology, i.e. 3 W/m °K, as given by McKenzie et al. (2005). Based on this new parameterization, maximum temperatures in the Moho of the NPM plateau are up to 969 °C, being 36 °C lower than those predicted by the reference model. (2) The volumetric rate of heat generation within the crust beneath the plateau. We tested a value in the lower extreme of the crustal radiogenic heat production range, 2 µW/m3 (Table 1). The new computed maximum Moho temperature in NPM is 938 °C; therefore, 67° C lower than in the reference model. A direct comparison of the computed temperatures leads to the conclusion that the crustal configuration (thickness and radiogenic heat content) exerts the primary role in controlling the lithospheric thermal field. However, the imposed reduction in the heat production does not agree with the main results of the previous density model by Gomez Dacal et al. (2017). In their study, the authors predicted a felsic and thus lighter and highly heat productive composition of the crust.
When compared with the available temperature measurements, the two models of the sensitivity analysis provide a worst fit to the data; the second model fitting even less than the first one. Even if lower than the reference model, the predicted Moho temperatures are still higher than those normally expected at this depth (~ 500/600 °C). We can thus conclude that the presence of a warmer crust beneath NPM plateau than at its surroundings is a robust feature of the thermal configuration of the lithosphere in the study area, even though the exact temperature magnitudes cannot be fully constrained by the available data. In summary, a warmer plateau's lithosphere reflects an increase in its crustal thickness and, although of secondary relevance, a shallower depth to the thermal LAB.
Accordingly, the combined structural and thermal models predict that the main differences between the NPM plateau and its surroundings can be found in the crust where the plateau is thicker and warmer. This difference could partly explain the topography variations observed in the area. The NPM plateau area is characterized by an elevation of more than 500 m a.s.l. (Fig. 1), a Moho depth of more than 40 km b.s.l. (kilometers below sea level, Fig. 2c) and an average crustal thickness of 45 km ( Fig. 2b; Gomez Dacal et al. 2017). Using this configuration, our model predicts that the 450 °C isotherm is shallower than 16 km below the plateau (Fig. 4b). In particular, the southern portion of the plateau, which is additionally influenced by a shallow thermal LAB (Fig. 2d), is also characterized by the highest topography (~ 1500 masl; Fig. 1). If we analyze these results on the base of local isostasy, a thicker crust would provide the cause of the higher elevations. This hypothesis was already advanced by Gomez Dacal et al. (2017), who suggested that the high elevation of the plateau is mostly a consequence of isostatic compensation in the presence of a thickened crust. If we also consider that this thicker than normal crust is also hotter than normal, local isostatic compensation would result in even higher elevation, due to the effects of thermal expansion on crustal densities (Turcotte and Schubert 2002). Consequently, the lower density in the domain of the plateau (caused by the combined effect of thicker and warmer crust) is able to generate the buoyancy force that is likely to support its present-day topography.
The analysis of the induced thermal uplift and subsidence in the center of the plateau (Fig. 7) demonstrates that at present, the contribution from thermal isostasy alone would amount to only ~ 150 m of relief. This value should be taken as a conservative estimate, given the main assumption that entered in its derivation. Indeed, our analytical calculation takes into account only changes in temperature, but ignores difference in composition and crustal thickness, which may provide a first-order contribution to the current topography (Hasterok and Chapman 2007). Hasterok and Chapman proposed a correction to normalize the crust that consists of a simple isostatic calculation. Prezzi et al. (2014) applied it successfully to the Altiplano-Puna plateau to isolate the thermal contribution to the computed elevation. If we apply a similar correction for our model by taking a reference model of 39 km of crustal thickness and 2.75 g/cm3 of density-information that coincides with a point with seismic data within Chilenia (Gomez Dacal et al. 2017)-we obtain a correction to the elevation of 900 m. Therefore, we arrive at a maximum of 100 m (closely to 10%) of the elevation (~ 1000 m), which is accounted for thermal effects only, that is in the order of the result obtained by thermal isostasy for the present day (Fig. 7). In summary, we find indications that the presence of a thicker and warmer crust beneath the NPM plateau, together with elevated mantle temperatures, is likely to be the primary support to its present-day topography and that the area is therefore in a state of dynamic isostatic equilibrium, but has not yet reached its final thermal equilibrium.
Another consequence of the elevated temperatures in the plateau lithosphere stems from its effects on the mechanical behavior of the lithospheric plate. Rocks at temperatures as high as those derived from the mantle in our analysis would rather deform by viscous flow (Turcotte and Schubert 2002). Gomez Dacal et al. (in press) reached a similar conclusion by means of a first-order rheological model, where they parameterized the different domains of Gomez Dacal et al. (2017) model by assigning them mechanical properties and calculating the mechanical behavior of the rocks under proposed temperature and pressure conditions. Similar observations as those summarized above are reported from the Colorado Plateau, USA. There, the plateau is located in a geodynamic scenario controlled by (1) the collision of the Farallon/Pacific ridge against North America, (2) the development of a transform margin (San Andreas fault system), (3) the detachment of Farallon plate that may open a slab window as it sinks and (4) the development of the Basin and Range province (B&R) under regional horizontal tectonic extension and uplift (Atwater and Stoke 1998;Bashir et al. 2011;Putirka and Platt 2012;Qashqai et al. 2016). Accordingly, at present, the Colorado Plateau has a regional tectonic setting that would correspond to the configuration that the NPM plateau would have had at the time when it started to form . The two areas also have similar characteristics: similar to the NPM, the Colorado Plateau, with an elevation of 1500 to 2000 masl, stands up 1000 m higher than its surroundings and is mostly internally undeformed, while its surroundings are deformed (e.g. Levander et al. 2011). In addition, the Colorado Plateau, again similar to the NPM, was affected by basaltic volcanism that erupted from the top of the plateau to the lower elevated areas in its surroundings after the uplift. Supported by a far larger amount of geophysical data, the Colorado Plateau is interpreted to have a thicker crust with respect to the surrounding B&R crust (43-49 km;Keller et al. 1979;McQuarrie and Chase 2000;Sheehan et al. 1997;Zandt et al. 1995), a pronounced low velocity zone in mantle below the southern portion of the plateau (Levander et al. 2011) and elevated temperatures (Hinojosa and Mickus 2002;Levandowski et al. 2018;Reiter et al. 1979;Roy et al. 2009). These similarities with the NPM plateau point to a similar origin and similar mechanisms at work maintaining their respective elevations. It is important to mention that even though a large amount of data is available for the Colorado Plateau area, there is still an open controversial debate on the mechanisms driving its tectonic evolution (Afonso et al. 2016). Nevertheless, there are several studies (e.g. Afonso et al. 2016) suggesting that most of the elevation in the area is likely to be the result of thermochemical buoyancy.

Insights into the geodynamic evolution of the NPM plateau
According to Aragón et al. (2011Aragón et al. ( , 2015, the time frame for the uplift of the NPM plateau is restricted to the Paleogene. Several hypotheses have been put forward concerning a thermal anomaly in the mantle of North Patagonia during this time and the connection to the generation of the Somún Curá basalts. These include processes as a transient hot spot (Kay et al. 1993(Kay et al. , 2007, rollback (Muñoz et al. 2000) and corner flow (De Ignacio et al. 2001). All of these models consider that the plateau is caused by topographic inversion of the Somún Curá basalts. An alternative explanation was given more recently by Aragón et al. (2011Aragón et al. ( , 2015, who suggested that the detachment of the Farallon plate and the opening of a slab window is the main cause of the thermal anomaly and proposed that the plateau is due to uplift before the eruption of the Somún Curá basalts. In any of these cases, the shallow thermal LAB in the southern portion of the modelled area (Fig. 2d), as derived from the conversion of tomographic velocities to temperatures, could be interpreted as the present-day imprint of the Paleogene thermal anomaly on the lithosphere configuration. A controversial evidence for such a model stems from the fact that the heating of the mantle would have had a regional influence and therefore additional reasons are to be found to explain the rather limited extent of the NPM plateau. Our modelling results suggest that the difference in crustal configuration and properties of the plateau compared to its surroundings could be such a reason. The plateau lithosphere is characterized by a thicker and warmer crust (Figs. 2c and 4), indicative of a less efficient heat conduction than around it (Fig. 6c). This setting implies a lithosphere in the plateau domain that is in a thermal state lagging in time with respect to the evolution of its surroundings. Therefore, at present day, the plateau area is in a state of dynamic thermal disequilibrium. This is the final result of the interplay between the Paleogene mantle dynamics (here imposed as a temporal perturbation in the LAB) and the thick and radiogenic enriched crust of the plateau that is responsible for the observed delay in its thermal evolution. It is important to mention that, as the radiogenic heat production is a volumetric parameter, the 3D approach is crucial to be able to reproduce a more realistic effect in the thermal conduction.
An independent observation supporting the abovementioned hypothesis comes from the age of magmatism through the area. The largest volume of volcanic eruptions was registered as Oligocene in age for the plateau area (~ 27 Ma; basalts of Somún Curá), while it was significantly older outside the NPM plateau (~ 52 Ma; Piedra Parada, SW NPM; Aragón et al. 2018). This would be consistent with the interpretation that the plateau was initiated by the heating of the mantle during the Paleogene plate tectonic configuration, regardless of the exact mechanism and origin of the anomaly, in combination with the specific configuration of the crust of the NPM plateau.
The results of the thermal isostasy demonstrated that the heating alone is not enough to elevate the surface by the required amount to match the current plateau elevation. At present, the largest part of the observed elevation (~ 90%) can be explained by isostatic equilibrium in the presence of a thicker crust beneath the plateau. One limitation in our interpretation stems from assuming that the crustal configuration did not change during the time of the uplift, here assumed because there is lack of evidence for crustal thickening at that time, when an extensional regional tectonic regime prevailed. Therefore, we must look for an additional process to fully explain the onset of the uplift. One option is delamination of the lithospheric mantle, which could have been caused by the same mantle thermal anomaly heating this portion of the lithosphere and thereby leading to its final detachment. This hypothesis was already proposed for other plateaus that also consisted of a thick crust and a warm lithospheric mantle, for example the Altiplano (Garzione et al. 2017) or Tibet (Molnar et al. 1993). This premise could also explain the delay in uplift with respect to the onset of the heating, because delamination would only occur once the lithospheric mantle is warm enough (Molnar et al. 1993). Another process that may cause the delamination of the lithospheric mantle is the drag force that may have generated the detachment of the Aluk plate after the subduction, following the hypothesis of Aragón et al. (2011). These two mechanisms-heating and drag force-could also have acted in combination. There might have been some other active co-players in such a complicated setting, as a dynamic support to the topography or isostatic rebound in response to active erosion. Since we cannot rule out the possible role of any of these mechanisms and we must notice that they would need to explain the basic observation that at present day, the area is still elevated.
Our simulations show that, irrespective of the specific time when the heating of the mantle took place or the absolute magnitude of the mantle temperature perturbation, the lithosphere beneath the plateau would still be in the process to relax such thermal perturbation ( Fig. 6a and b), thus indicating that the plateau is still not in steady state at present.

Conclusions
We evaluated the relation between the present-day temperature distribution, as well as the thermal evolution of the NPM plateau with its current topography and suggested geodynamic evolution scenarios. With this aim, we conducted 3D steady state and transient thermal simulations of the whole lithosphere of the NPM plateau and its surrounding areas using a data-based geological 3D model (Gomez Dacal et al. 2017).
The simulations of the steady-state thermal field and of a transient perturbation in the thermal LAB indicate that the lithosphere of the plateau is significantly hotter as compared to the one in its surroundings. These results are robust even for lateral variations in crustal lithology and changes in thermal properties as found from sensitivity analyses. The results are mainly due to the thickness of the radiogenic crystalline crust and secondary due to the locally shallower depth of the thermal LAB. The 100° C isotherm, as a proxy for the shallow thermal field, is 0.1 to 0.8 km shallower inside the plateau than outside. Likewise, the 450° C isotherm, as a proxy for the deeper thermal field, is from 1 to 6 km shallower in the domain of the plateau as compared to its surroundings. In both cases, the shallower depth of the isotherms correlates spatially with elevations higher than 500 m a.s.l. found in the plateau. At Moho depth, the temperature is up to 300° C warmer inside the NPM plateau as compared to the surroundings.
Considering the crustal structure together with the newly obtained results of the temperature distribution, the models suggest not only a warmer but also a thicker crust below the NPM plateau than in its surroundings, which would be isostatically more buoyant and therefore explain the elevated present-day topography of the area.
The transient simulations indicate that the proposed regional heating of the mantle during Paleogene could be still affecting, to a greater extent the plateau than its surroundings, at present. Overall, the insulating effect of the warmer and thicker crust in the plateau slows down the conduction process and results in a delayed cooling of the initial thermal disturbance, giving a distinguished characteristic to the NPM plateau. This also implies that the NPM plateau has not yet reached thermal equilibrium today, whereas the surrounding areas would be closer to steady state than the plateau.
Thermal isostasy analyses show that the effect of the temperature alone could not be the cause of the uplift of the plateau. A lithospheric mantle delamination, that could also be generated by the heating in the mantle, might also have affected the area and generated the onset of the uplift. At present, the area is still elevated, which can be explained by its isostatic equilibrium and the thermal field that is still in its way to equilibrium.
Therefore, according to our results, the elevation of the plateau could be the result of the interplay between two superposed effects: (1) heating in the mantle produced by the Paleogene tectonics and (2) a thickened radiogenic crust in the domain of the plateau.
Funding Open Access funding enabled and organized by Projekt DEAL.

Compliance with ethical standards
Conflict of interest We hereby confirm that there are no known conflicts of interest associated with this publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.