Using a Paleosurface to Constrain Low‐Temperature Thermochronological Data: Tectonic Evolution of the Cuevas Range, Central Andes

Dispersion of low‐temperature thermochronologic data from nine samples collected on a deformed paleosurface preserved on the Cuevas range (Central Andes) can be exploited to unravel complex thermal histories. The nine samples yielded data that have both intersample and intrasample dispersions; the data set includes apatite fission‐track ages (180–110 Ma), mean track lengths (11–13 μm), apatite helium (10–250 Ma), and zircon helium ages (180–348 Ma). We ran inverse thermal history models for each sample that reveal spatial variations of the Miocene reheating along the paleosurface. Next, we ran a multiple‐sample joint model to infer a common form for thermal history for all samples. Our results suggest that initial exhumation during the Famatinian orogeny was followed by a residence between ~2.5 and 7.0 km depth during the Paleozoic and the Triassic. The onset of Mesozoic rifting was responsible for an increase of the geothermal gradient and extensive horst exhumation, which brought the basement of the Cuevas range close to the surface (~1–2 km) in the Late Jurassic. Between the Late Cretaceous and the Paleocene, the combination of low relief, a humid climate, and low erosion rates (0.006–0.030 km/Ma) facilitated the development of the Cuevas paleosurface. During the Miocene, this paleosurface experienced differential reheating with a high geothermal gradient (>25 °C/km) due to the sedimentary cover and local magmatic heat sources. During the Andean orogeny, in the Pliocene, the Cuevas paleosurface was deformed, exhumed, and uplifted.


Introduction
Ancient landscapes preserved in mountain belts and cratonic areas record past surface conditions and represent a structural constraint for understanding subsequent tectonic cycles (Everglades, 2013;Ollier, 1991;Phillips, 2002;Rabassa & Ollier, 2014;Thiry et al., 2014). Low-temperature thermochronologic systems (<200°C) can be used to reconstruct long-term landscape evolution and thus the evolution of ancient landscapes (Braun et al., 2012;Jordan et al., 1989).
Multiple-sample thermal history models (MSMs) use thermochronologic data from samples collected in different positions within a rock body, considering variations in the temperature experienced by the samples (Gallagher, 2012). Ancient landscape features such as paleosurfaces have the potential to provide an independent structural constraint for multiple-sample thermal histories of rock bodies, since paleosurfaces represent a laterally continuous past surface level. MSMs from samples collected along a continuous paleosurface offer the opportunity to do spatial thermal history modeling from the "x-y" domain; this modeling approach presents an alternative to traditional MSMs in which only the structural depth is considered ("z" domain) (Gallagher et al., 2005). Moreover, a MSM from a paleosurface has the potential to bracket when the ancient landscape was formed and spatially constraint the subsequent thermal evolution.
The Sierras Pampeanas province is located in the Central Andes of Argentina between 26°S and 33°S ( Figure 1). Several paleosurfaces called the "Pampean Peneplains" have been described on top of the basement highs. The age of these surfaces has been poorly constrained. Moreover, the lack of temporal controls complicates understanding the processes involved in the formation of these paleosurfaces and their relation with the tectonic evolution of the Sierras Pampeanas Dávila & Carter, 2013;Jordan et al., 1989;Rabassa & Ollier, 2014;. ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
We acquired apatite fission track (AFT) and apatite and zircon (U-Th-Sm)/He thermochronology (AHe and ZHe, respectively) data from nine samples collected on a paleosurface on the Cuevas range, in the Sierras Pampeanas ( Figure 2). Although the samples were collected along a continuous paleosurface, AFT and AHe data exhibit intrasample and intersample dispersions. We performed a joint MSM which succeeds in explaining the observed dispersion. The MSM exploited observed dispersion to reconstruct the Paleozoic to Pliocene thermal evolution of the Cuevas range, reconstruct the exhumation history of the range, and quantify the thermal evolution of the Miocene basin which partially overlies the paleosurface.

Tectonic and Thermal Evolution of the Sierras Pampeanas
The Sierras Pampeanas is a tectonomorphic province located between 28°and 33°S in the Central Andes. This province is composed of several broken foreland basins characterized by thick-skinned deformation ( Figure 1) (Caminos, 1979;Jordan et al., 1989;Kley et al., 1999;Ramos, 1999). To the north, this province  Jordan and Allmendinger (1986). The red square denotes the study area presented in Figure 2. CR = Cuevas Range; QR = Quilmes Range; CAB = Campo-Arenal Basin; CCR = Cumbres Calchaquies; AR = Aconquija Range. Dashed black line denotes the extent of the Salta rift basins (after Marquillas et al., 2005;Schmidt et al., 1995). is bounded by the elevated Puna plateau and the inverted Salta rift basin (Ramos, 2017;Salfity & Marquillas, 1994). The Sierras Pampeanas are characterized by Paleozoic to Pliocene ZHe, AFT, and AHe ages Coughlin et al., 1999;Jordan et al., 1989;Löbens et al., 2013;. This thermochronological record has been interpreted to result from several contrasting tectonic and erosional processes including postorogenic cooling and relief reduction, extensive exhumation related to compressional phases, and horst exhumation during several rifting cycles (Bense et al., , 2017Coughlin et al., 1999;Löbens et al., 2013). These contrasting tectonic settings also imply significant spatial and temporal variations of the geothermal gradient in the Sierras Pampeanas (Dávila & Carter, 2013;. The basement of the Sierras Pampeanas is mainly composed of Proterozoic-Paleozoic crystalline rocks (Ramos, 2008;Rapela et al., 1998). The Famatinian orogeny (550-450 Ma) was an accretionary event characterized by extensive ductile deformation, basement exhumation, and arc volcanism (Bense et al., 2013;  (Coughlin et al., 1999;Löbens et al., 2013;Mortimer et al., 2007;Zapata et al., 2018), the black ellipse contains the data presented in this contribution. Color code denotes ages; squares represent AFT ages and circles AHe ages. The black arrow denotes the location and the viewpoint of the Figure 3a landscape image. (b) Simplified stratigraphic section measured in Section 1 modified from Muruaga (2001), the thickness for these units in Section 2 are also presented (Bossi & Muruaga, 2009). The location of both sections is presented in (a). (c) Modified cross section after Bossi and Muruaga (2009); the location of the section is presented in (a). Ramos, 2008;Rapela et al., 1998). After the Famatinian orogeny, different tectonic settings have been proposed for the Sierras Pampeanas. As the SP is characterized by an unconformity between the basement and Mesozoic to Cenozoic strata, the observations supporting the inferred tectonic settings are based on the geology preserved along strike. To the south of the Sierras Pampeanas (>30°S), the Carboniferous has been interpreted as a period of orogenic collapse and relief reduction Jordan et al., 1989). During the Permo-Triassic, the termination of the magmatism and basement cooling inferred from inverse thermal modeling have been interpreted as a result of flat-slab subduction and associated compression in the southern Sierras Pampeanas e.g., Ramos & Folguera, 2009). Farther north in the Velasco and the Aconquija ranges, a different Paleozoic evolution has been proposed. In this region, based on the lack of sedimentary cover, deformation, and basement exhumation, the Carboniferous to Triassic has been interpreted to mark a period of tectonic quiescence .

Tectonics
During the Mesozoic, the Sierras Pampeanas province was characterized by the development of several rift basins related to the opening of the South Atlantic Ocean (Horton, 2018;Oncken et al., 2006). Late Triassic to Cretaceous rifting cycles were responsible for extensional magmatism and basement exhumation in local depocenters, such as the Salta, Sierras Chicas, and El Gigante basins ( Figure 1) (Bense et al., 2017;Marquillas et al., 2005;Schmidt et al., 1995;Viramonte et al., 1999). In the Velasco range, Early Cretaceous high thermal gradients (35-50°C/km) related to the rifting phases have been inferred using thermal modeling of low-temperature thermochronologic data . During the Late Cretaceous and the Paleogene (100-40 Ma), postrift thermal relaxation facilitated the accumulation of clastic sediments in the Sierras Pampeanas Cretaceous rift basins (Marquillas et al., 2005;Schmidt et al., 1995). During the Mesozoic and the Paleocene, several episodes of humid climatic conditions were responsible for the development of paleosols described on the sedimentary cover and the basement of the Sierras Pampenas (Caminos, 1979;Rabassa & Ollier, 2014). Late Cretaceous to Paleocene paleosols resulting from warm and humid conditions have been documented and dated in the Salta basin (Andrews et al., 2017;M. Do Campo et al., 2007).
The Eocene and Oligocene record of the Sierras Pampeanas evolution is scarce; whereas foreland basin depocenters associated with the early stages of the Andean orogeny have been reported in this province (Bossi & Muruaga, 2009;Fosdick et al., 2017). During the Miocene, several fine-grained lacustrine and fluvial successions were deposited and subsequently overlain by coarse sediments related to the later phases of the Andean orogeny (Bonini et al., 2017;Bossi et al., 2001;Dávila & Astini, 2007;Strecker et al., 1989).
Despite the high relief of the Sierras Pampeanas, this province is characterized by old (>40 Ma) AFT and AHe ages, with the exception of Miocene ages from the Aconquija and Famatina ranges ( Figure 2). Bense et al. (2013) have suggested that most of the relief formation and rock exhumation took place during the Permo-Triassic based on thermal modeling of these old thermochronologic ages. In contrast, several studies have proposed Miocene relief development based on the Miocene stratigraphic record and thermal modeling (Bense et al., 2017;Davila et al., 2012;Mortimer et al., 2007;Zapata et al., 2019).
To explain these old thermochronological ages, authors have also suggested that there was a low Cenozoic geothermal gradient (18-23°C/km) in the Sierras Pampenas due to the reduction of the mantle wedge during flat subduction (e.g., Collo et al., 2017;Dávila & Carter, 2013); this means that there has not been sufficient exhumation to expose the relatively deep-lying AFT and AHe closure isotherms so that Cenozoic exhumation has been underestimated. However, thermochronological studies of the Mio-Pliocene sedimentary fill and the crystalline basement have suggested geothermal gradients between 25 and 35°C/km in the Sierras Pamepeanas, suggesting that flat-slab subduction has not reduced the geothermal gradients during the Miocene .

Ancient Landscapes in the Sierras Pampeanas
Several paleosurfaces have been described on top of the Pampean ranges and called the "Pampean peneplains." These paleosurfaces are often dismembered, tilted, and folded as a result of the Andean orogeny, which is responsible for uplifting these basement ranges (Carignano, 1999;Jordan et al., 1989; 10.1029/2019TC005887 Tectonics . The use of the term peneplain has genetic implications, such as being formed by subaerial erosion, close to the sea level, and having a continental-scale extent (Phillips, 2002). The term peneplain has been used carelessly in the literature, especially in the Sierras Pampeanas where the genetic origin of these paleosurfaces is still obscure, as has been discussed by several authors Fairbridge & Finkl, 1980;Phillips, 2002;Rabassa et al., 2010). Another type of ancient planations is a preserved etchplain (Cui et al., 1999;Phillips, 2002), which is a surface defining the limit between fresh and weathered rocks in weathering profiles; these surfaces can form at up to 1 km depth as a result of chemical denudation under tropical and subtropical climates (Rabassa & Ollier, 2014;Thomas, 1968).
Based on field observations, some authors have considered these ancient landforms in the Sierras Pampeanas to be part of a single Paleozoic surface dismembered by tectonic forces during the Miocene (Beltramone, 2007;Criado Roque et al., 1981;González Bonorino, 1972;Sayago, 1986;Schmieder, 1921). Subsequently, other authors have proposed that these surfaces were formed at different times Jordan et al., 1989;Rabassa et al., 1996Rabassa et al., , 2010. Ages between the Paleozoic and the Cenozoic have been inferred for these paleosurfaces based on field relations and thermochronological data, which constrain the time in which the basement of the surfaces was close to surface temperatures (Bense et al., 2017;Jordan et al., 1989;Rabassa et al., 2010;. Additionally, different mechanisms have been proposed for the origin of these landforms in the Sierras Pampeanas, including prolonged subaerial exposures (Fisher et al., 2002;González Bonorino, 1972), the product of horst exhumation during extensional tectonism (Bense et al., 2017;, and as etchplains formed in humid tropical and subtropical climates (Carignano, 1999;Rabassa et al., 2010;Wayland, 1933). In the tectonic model presented by Bense et al. (2013), these surfaces have been considered to have formed after the creation of relief in the Sierras Pampeanas during long periods of very slow cooling. The lack of refined temporal constraints on different paleosurfaces has complicated correlating between them and establishing their relationship with climate and tectonics.

Geological Setting of the Cuevas Range
The Cuevas range is 3,000 m above sea level with 800 m of local relief and a NNE trend ( Figure 2). This range is located on the northwestern margin of the Campo-Arenal basin in the northern part of the Sierras Pampeanas province. The granitic basement of the Sierra Cuevas range is equivalent to Ordovician granitoids formed between 550 and 450 Ma during the Famatinian orogeny (Ramos, 2008;Rapela et al., 1998). Approximately 1.4 km of Cenozoic clastic strata, grouped into the Hualfin, Las Arcas, Chiquimil, and Andalhuala formations, were unconformably deposited on top of the paleosurface (Figure 2b) (Bossi & Muruaga, 2009). This sedimentary basin fill is associated with the development of foreland basins related to the Cenozoic Andean compressional phases (Bossi et al., 2001;Mortimer et al., 2007).
The age of these Cenozoic units has been constrained south of the Cuevas range where the Cenozoic basin fill is~1 km thicker (Section 1 in Figure 2). The age of the Hualfin Fm. is bracketed between Eocene and Miocene based on stratigraphic relations and regional correlations (Bossi & Muruaga, 2009). The Las Arcas Fm. was unconformably deposited on top of the Hualfin Fm.; the lack of syndepositional volcanism suggests that deposition of this unit occurred before of magmatism started at 12 Ma. Miocene volcanism in the Campo Arenal basin was related to the Farallon Negro volcanic complex, which is located~15 km south of the Cuevas range ( Figure 2a).
The Chiquimil Fm. was intruded by subvolcanic rocks with an Ar-Ar age of 9.14 ± 0.02 Ma (Figure 2b) (Bossi & Muruaga, 2009;Sasso, 1998). A tuff from the Andalhuala Fm. west of the Cuevas-Hualfin range has an Ar-Ar age of 7.14 ± 0.02 Ma (Latorre et al., 1997;Marshall & Patterson, 1981). Coarse conglomerates of the Corral Quemado Fm. and the "Punaschotter" were deposited during the Pliocene and the Quaternary (Bonini et al., 2017;Bossi et al., 2001;Bossi & Muruaga, 2009). A tuff collected in the base of the Corral Quemado Fm. west of the Cuevas range has an Ar-Ar age of 3.60 ± 0.02 Ma (Figure 2b) (Latorre et al., 1997;Marshall & Patterson, 1981). Depositional patterns of the Corral Quemado Fm. suggest that the deposition of this unit was coeval with deformation of the Cuevas range (Bossi et al., 2001;Mortimer et al., 2007). Seismic records from the Campo Arenal basin suggest homogeneous Cenozoic sedimentary cover in the segment of the paleosurface where the samples were collected.

Sampling
Nine samples were collected from a paleosurface carved into granitic basement rocks on the west side of the Cuevas range ( Figure 2). The samples were collected at elevations between 2,400 and 3,400 m above sea level, in different positions along strike ( Figure 3a). All of the samples were collected from the top of a continuous surface within an area of~20 × 5 km, with the exception of Sample A2, which seems to be part of a faulted segment of this paleosurface. Regardless of the mechanisms involved in the formation of the Cuevas paleosurface, this geomorphic feature constitutes a reference structural level for the basement of the range.

AFT and U-Th/He Methods and Causes of Dispersion
Thermochronology is based on the quantification of the products of radioactive decay after the closure of a radioisotopic system (Dodson, 1973;Harrison, 2005;Reiners, 2005). The closure temperature marks the retention of radiogenic products, which varies depending on the cooling rates, mineral kinetic variability, and other conditions intrinsic to each system (Ehlers & Farley, 2003;Harrison, 2005;Lisker et al., 2009;Reiners et al., 2002). Slowly cooled samples residing at partial annealing or partial retention temperatures for long periods of time, where daughter decay products are lost by diffusion or annealing, exhibit more age dispersion due to the small changes in the capacity of each mineral to retain daughter products (e.g., Barbarand et al., 2003;Shuster et al., 2006).
The AFT method is based on the quantification of damage to the crystal lattice (tracks) that results from the spontaneous fission of 238 U. Fission tracks can be partially or fully annealed at temperatures between 60 and 120°C; this interval is known as the partial annealing zone (Wagner et al., 1989). Fission tracks shorten within the partial annealing zone; therefore, the track length distributions and the mean track length (MTL) can be used as proxies to quantify the amount of annealing experienced by the apatite crystal (Green et al., 1985). The apatite annealing resistance depends strongly on the chlorine content, which controls the internal mineral kinetics. An equally useful kinetic indicator is the crystal resistance to the acid used to reveal the tracks (etching). Therefore, the C-axis parallel length of the etch pit (Dpar) can be used as a proxy to quantify resistance to annealing (Barbarand et al., 2003;Carlson et al., 1999;Ketcham et al., 1999). The procedures involved in sample preparation and age determination in the AFT system can have a significant influence on data dispersion. Documented sources of dispersion are related to the etching procedures, equipment, and the method used for the age determination, the presence of mineral defects can be confused with tracks, and thus, the experience of the operator is important (Enkelmann et al., 2012;Hurford, 1998;Ketcham et al., 2009;Sobel & Seward, 2010).
Apatite U-Th-Sm/He thermochronology is based on the production and accumulation of helium that forms from the alpha decay of U, Th, and Sm. This method provides information on the thermal history between 40 and 100°C in apatites; this temperature interval is known as the apatite partial retention zone (Farley, 2002). Zircon U-Th/He thermochronology quantifies the accumulation of helium from the decay of U and Th to obtain thermal history information between 140 and 220°C; this interval is known as the zircon partial retention zone (Guenthner et al., 2013). The helium retentivity of the crystal can be affected by different factors including effective uranium content (eU = U + 0.235*Th), size, geometry, and cooling rate (Brown et al., 2013;Flowers, 2009;Guenthner et al., 2013;Shuster et al., 2006).
Radiation damage reflects deformation of the internal lattice associated with spontaneous fission and alpha decay. Alpha damage creates zones of damage in the crystal lattice which trap helium, thereby increasing gas retentivity in the mineral; after the accumulation of significant radiation damage, these helium traps can interconnect and increase helium diffusivity (Flowers, 2009;Flowers et al., 2007;Gautheron et al., 2009;Guenthner et al., 2013;Shuster et al., 2006). For apatite, radiation damage-induced age dispersion is controlled by three factors. (1) The amount of effective uranium; crystals with more effective uranium produce more alpha particles and therefore more lattice damage. (2) Long residence at low temperatures (below the partial retention zone) allows more damage to accumulate. (3) Crystals reheated to temperatures within the partial retention zone experience helium loss correlatable with the accumulated damage (Ault et al., 2009;Flowers, 2009;Flowers et al., 2007;Gautheron et al., 2009;Shuster et al., 2006;). Available models for radiation damage annealing assume that alpha damage anneals analogously to fission tracks during the reheating events (Flowers, 2009;Gautheron et al., 2009;Guenthner et al., 2013).

10.1029/2019TC005887
Tectonics There is a correlation between the size of the crystal and helium diffusivity and hence the U-Th/He age (Brown et al., 2013;Reiners & Farley, 2001). In the case of radiation damage and size controls in the AHe and the ZHe systems, the shape of the nonlinear crystal size and eU versus age correlation curves depends on the thermal path experienced by the sample (Brown et al., 2013;Flowers et al., 2007). Therefore, the size and eU trends can be inverted to reconstruct and refine a specific thermal path (Ault et al., 2009;Brown et al., 2013;Flowers, 2009;Flowers & Kelley, 2011;Guenthner et al., 2013).
Age dispersion can also be related to internal crystal eU zonations responsible for retentivity zoning and inaccurate alpha ejection corrections (Ft) (Ehlers & Farley, 2003;Flowers, 2009;Hourigan et al., 2005;Vermeesch, 2008). Finally, additional sources of dispersion are related to helium implantation from neighboring high eU grains and microinclusions with high U and Th content (Fitzgerald et al., 2006;Vermeesch et al., 2007). Detailed procedures of the ZHe, AHe, and AFT methods are presented in supporting information Texts S1 and S2.

Thermal History Modeling
Thermal history modeling was conducted with the QTQt software, which uses a Bayesian transdimensional statistical approach. This allows the extraction of thermal histories using multiple thermochronometers from multiple samples located in different structural positions within a rock body (Gallagher, 2012). We used the radiation damage model from Flowers et al. (2009) for the AHe models and from Guenthner et al. (2013) for the ZHe models. The fission track data were modeled with the annealing model from Ketcham et al. (2007). Only U-Th/He single grain aliquots with reproducible ages or dispersed ages with possible size or eU age controls were included in the models. We consider ages to be reproducible when the 1σ SD is <20% of the mean age Flowers & Kelley, 2011). Unfortunately, the angle between the C axis and the confined tracks was not measured for every sample; therefore, uncorrected length distributions were incorporated in the models. Four single-sample models (SSMs) including the C-axis track length distributions were performed for comparison ( Figure S3).
In QTQt multiple-sample thermal history models (MSMs), the samples need to be arranged according to the relative temperatures experienced during the thermal history; therefore, one must identify the samples with the coldest and the hottest thermal histories (cold and hot samples, respectively). The temperature difference between the cold and the hot samples at any specific point in the thermal history is defined as the temperature offset (Gallagher, 2012).
Two geological constraints were included in the models. Since the crystallization ages of the granites have been constrained to be between 550 and 450 Ma (Ramos, 2008 and references therein), we set a high-temperature constraint to represent the time of mineral crystallization. The models were allowed to find solutions between 0 and 550 Ma. A constraint to force the model to be close to surface temperatures between the Eocene and the Miocene (50-10 Ma) was also included; this constraint is based on the depositional age of the Hualfin Fm., which was deposited on top of the Cuevas paleosurface ( Figure 3a) (Bossi & Muruaga, 2009). The data used as constraints and the parameters incorporated in the thermal models are presented in Table S3 and the predicted values in Table S4 following the recommendations presented in Flowers et al. (2015) and Gallagher (2016). The likelihood chain is presented for each model in Figure S6. Additionally, all QTQt files containing the raw data and the results of the models are presented in the Data Set S1. To visualize the modeling result, we present the expected model and the associated uncertainty. The expected model is a mean of the sampled models weighted by the posterior probability (Gallagher, 2012). We present the data predicted from the expected model and the 2 sigma interval of the predicted data.

Low-Temperature Thermochronology Results
We used the U-Th-Sm/He method to date 19 single apatite grains from five samples and four zircon grains from Sample C3. Summarized U-Th-Sm/He data are presented in Table 1. The complete data set is presented in Table S2. AHe ages vary between 8 and 280 Ma. Apatite grains had an equivalent spherical radius between 36 and 62 μm and eU contents between 38 and 114 ppm (Figure 3b). The five samples with AHe data exhibit intrasample dispersion (1σ SD > 20% of the mean age). In Samples A2 and C8, two single grain ages older than the AFT age and with relatively low eU were considered outliers. Most aliquots in Samples C3, C6, and H2 exhibit positive relations between age and eU; therefore, two single grain age from the Samples H2 and C6, which do not fit the intrasample eU versus age or grain size versus age trends, were excluded from the models (Figure 3b). In total, four aliquots were excluded. These outliers can be related to crystal zonation which was not controlled in this study (Ault & Flowers, 2012;Flowers & Kelley, 2011); additional causes of error could be related to high eU microinclusions and helium implantation (Vermeesch et al., 2007). Sample C3 has ZHe single grain ages between~370 and~180 Ma, with eU values between 350 and 1,000 ppm, exhibiting a negative correlation between age and eU ( Figure 3b).
AFT data are summarized in Table 1. The complete data set is presented in Table S1, and detailed counting data are presented in the QTQt text files in Data Set S1. All samples passed the P (χ2) test (>5%). Therefore, the distribution of single grain ages in each sample is consistent with a single age population (Galbraith & Laslett, 1993;Green, 1981). Ages vary from 181.5 ± 14.5 to 111.5 ± 7.4 Ma, and unprojected confined fission MTLs vary between 11.3 and 12.7 μm ( Figure S2). Uncorrected Dpar values vary between 1.7 and 2.6 μm. No evident trend can be observed in the fission-track age versus elevation plot ( Figure S1). Some of the slightly older samples (A2 and C7) have higher Dpar values; however, the oldest AFT samples (C3 and C2) have lower Dpar values (Figure 3c).

SSMs
Although samples were collected along the same paleosurface, the AFT and AHe ages exhibit significant intrasample and intersample dispersions; therefore, it is not obvious if these samples experienced the same thermal history. Intersample dispersion can be related to differences in mineral kinetics (e.g., Dpar, eU, and grain size) and/or to differences in the thermal histories along the surface. Thus, simple examination of the data is not enough to evaluate the causes of dispersion. Each sample was modeled individually and considered as an independent data set in order to compare the thermal histories along the paleosurface.
The expected SSMs are presented in Figure 3d, and the complete individual models are presented in Figure  S4. The comparisons between the observed and the predicted data are presented in Figure 4c and Table S4. All models included an AFT age, associated Dpar, and an uncorrected confined track length distribution. ZHe data were not included in the SSMs since these results come from a single sample. Most of the expected paths cool below 120°C between 200 and 150 Ma (Figure 3c). Model C2 has older exhumation times (250-300 Ma). After the Jurassic cooling, all expected models paths exhibit Cretaceous to Paleogene cooling followed by Miocene reheating between 40 and 100°C.

Tectonics
The expected SSMs exhibit overlapping Mesozoic thermal histories. However, the SSMs exhibit differences in the Miocene reheating peak presenting values between~45 and 95°C (Figures 3d and S4). The comparison between the probability distributions of the Miocene (3-6 Ma) maximum reheating temperatures from the SSMs suggests that the reheating temperatures span between 40 and 95°C; samples with higher and lower amounts of reheating are clearly distinguishable. Sample H2 exhibits the highest reheating temperature (~80 to 93°C) and sample C3 the lowest (~50 to 85°C) (Figure 3e).
The SSMs poorly predict the AHe data with R 2 = 0.34 in the observed versus predicted plot and with residual values up to 30 Ma (Figure 4c). The AFT data are well reproduced by the models, with R 2 = 0.99 and residual values lower than 13 Ma (Figure 4c). The predicted MTLs are higher than the observed with residuals between −0.5 and −1 μm and an R 2 = 0.93.  Figure 3e). These differences in the Miocene thermal histories could be related to variations in the sedimentary cover thickness, Miocene magmatism, and/or differential basement exhumation. Considering a geothermal gradient of 40°C/km, variations in the amount of Miocene reheating due to differential burial would require an increase of at least 1 km of the Cenozoic strata along the surface; a lower gradient would require thicker strata. However, the seismic record does not record significant changes in the sedimentary thickness along the sampled segment of the paleosurface (Mortimer et al., 2007). The preservation of a continuous paleosurface suggests that there were no significant variations in the amount of basement erosion between the samples. Thus, the samples experienced roughly the same amount of Miocene exhumation. Miocene volcanism is the most plausible explanation for the differential Miocene reheating. Several studies have shown how the emplacement of magmatic bodies can create partial resetting aureoles of the low-temperature thermochronological systems in the host rock; the size of these aureols is directly related to the duration of magmatism, the depth of emplacement, and the size of the intrusives (e.g., Karlstrom et al., 2019;Murray et al., 2018).
QTQt and HeFTy only allow 1-D modeling; these platforms do not allow the temperature offset between samples to vary in the x, y, and z directions; to overcome this limitation, several approaches have been proposed in order to understand spatial variations in rock thermal histories (Gallagher et al., 2005;McGregor et al., 2013). In this contribution, we use a QTQt pseudo-elevation profile to perform a 1-D model along the Cuevas paleosurface to test the hypothesis that the data from all of the samples can be explained by a single thermal history and to quantify possible lateral variations of Miocene reheating.
In traditional unidimensional thermal history modeling, the samples are collected at different structural depths, and the assessment of the relative temperatures experienced by the samples is typically based on a positive correlation between structural depth and cooling ages (Figure 4a). Therefore, the sample with the shallowest structural depth with respect to a reference level (e.g., sedimentary unconformity, elevation above sea level, or borehole depth) is considered to have a colder (cold sample) thermal history with respect to the sample with the greatest structural depth (hot sample) (Figure 4a). However, intersample dispersion along the Cuevas paleosurface can be related to Miocene magmatism, which is not necessarily controlled by simple spatial relations between the samples (Figure 4a). In consequence, the samples cannot be arranged according to their position along the surface. Since the hot and cold samples were clearly identified in the SSMs, we arranged the samples in a pseudo-elevation profile according to the relative Miocene reheating temperatures (Figure 3e). The sample with the coldest thermal history was located on top of the profile, and the sample with the hottest thermal history was located at the bottom. Pseudo-elevations were arbitrarily assigned between 0 and 800 m (Figure 4b), with the temperature offset allowed to vary between 0 and 60°C (Table S3).
In QTQt, the assigned elevation does not control the cold and the hot thermal histories or the thermal offset, since inverse modeling extracts thermal histories in a time-temperature space independent of the assigned elevation (Gallagher, 2012). In the SSMs, the samples with maximum and minimum reheating (cold and hot samples) are clearly distinguishable. Conversely, differences in the amount of reheating between some of the intermediate samples are not clear, especially in the models without AHe data. Thus, the position assigned to these samples may be ambiguous. However, thermal histories from these samples will not exceed the values of the cold and hot samples; therefore, the variations in the position of the intermediate samples will not significantly affect the overall thermal history or the temperature offset.
Sample C2 was excluded because we consider this sample to have an anomalously old AFT age and low Dpar values that cannot be explained by differences in mineral kinetics or spatial variations along the paleosurface. ZHe data from the cold sample were added to the MSM. The expected MSM (Figure 5a) exhibits Cambrian cooling at~540, followed by reheating during the Late Jurassic. Afterward, a fast cooling event to below 65°C took place between 160 and 140 Ma. After 140 Ma (Figure 5b), the model shows reheating above 65°C at~100 Ma, followed by slow cooling until the Miocene. A final reheating event started in the Miocene and reached a peak of~64°C for the cold sample and~84°C for the hot sample expected models. The thermal offset increased during the Miocene to a maximum of 50°C (Figures 5b and S4). The model suggests a final fast cooling path in the Pliocene (Figure 5a).

Tectonics
To evaluate the quality of the model, the fit between the observed and the predicted data in the MSM can be compared with those obtained in the SSMs (Figure 4c). The expected MSM improves the fit between the observed and the predicted AHe single grain ages with respect to the SSM, with R 2 = 0.83 and lower residual values for 13 of the 15 modeled AHe aliquots (Figure 4c). The MSM predicted AFT ages maintain the data fit (R 2 = 1) and the residual values with respect to the SSMs (Figure 4c). The expected MSM predicted MTL values have a lower correlation of R 2 = 0.64 with respect to the SSM, but the residual values are lower for seven of the eight samples (Figure 4c). The expected MSM poorly predicts the ZHe data; however, the 2sigma interval of all the accepted models overlaps with the observed data ( Figure 4b). The expected model is an average of the accepted models but not necessarily the thermal path which reproduce better the observed data.  Text S4 and  Table S5. The black boxes denote the thermal gradient estimates in the Sierras Pampeanas from other studies (Dávila & Carter, 2013;. Calculations of the exhumation history and the thermal gradients are discussed in section 7.2 and presented in detail in Text S3 and Table S5. 10.1029/2019TC005887 Tectonics 7. Discussion

Validity and Limitations of the Joint MSM
The joint MSM succeeds at extracting a common thermal history of the Cuevas paleosurface and reproducing the observed data. This result validates the assumptions and premises used to construct the MSM. This model is able to explain the observed dispersion as a result of radiation damage effects combined with differential Miocene reheating. The imposed arrangement of the samples allows the model to reconstruct and quantify the Miocene differential reheating along the paleosurface. The MSM reproduces the observed AFT data without a thermal offset, during the Cretaceous and the Paleogene, suggesting that during this period of time, the paleosurface experienced a similar thermal history. Therefore, the imposed arrangement of the samples according to the Miocene reheating did not affect the Cretaceous to Paleogene thermal history. The pre-Jurassic thermal history was constrained by ZHe data from only the cold sample. Thus, it is not possible to reconstruct the spatial evolution of the paleosurface during this period of time.
Compared to the SSMs, the MSM better reproduces the observed data and exhibits a more constrained thermal history, especially the AHe data. The increase in the number of AHe aliquots allows the MSM to refine and reconcile the thermal histories between 100 and 40°C, resulting in a more constrained thermal history which better predicts the observed data. The MSM profits from the eU controlled intrasample and intersample age dispersions. This improvement in reproducing the observed data shows that larger data sets are required to unravel complex thermal histories. The efforts to reconstruct this thermal history using smaller data sets (SSMs) resulted in highly variable and less constrained thermal histories.
We are aware that the increase in the amount of data may lead to more constrained models as suggested by Vermeesch and Tian (2014). However, more constrained models do not necessarily improve the fit between observed data and predicted values, which is the natural way to evaluate inverse models. Therefore, since the MSM maintains or improves the fit between the predicted and the observed data using a significantly larger data set ( Figure 5), we consider the MSM rather than the SSMs to be the more reliable representation of the thermal history of the Cuevas range.
Our new multiple sample modeling approach presents an innovative way to model data sets in which past thermal histories cannot be directly related to spatial relationships between the samples. This approach relies on the proper identification of variations in the single-sample thermal histories, which, in some cases, can be challenging and subjective. However, inadequate sample arrangements will result in poor data fit, providing a way to assess the model validity.

Inferring an Exhumation History for the Basement of the Cuevas Range
The MSM presents a tightly constrained Paleozoic to Pliocene thermal history of the Cuevas range. This refined thermal history combined with the seismic data and the surface stratigraphy offers the opportunity to bracket the basement exhumation history and the past thermal gradients of the range. The expected MSMs suggest two different reheating events at~160 and~100 Ma. However, no Mesozoic sedimentary rocks have been reported in the Campo Arenal basin (Bossi & Muruaga, 2009;Mortimer et al., 2007). Conversely, significant changes in the geothermal gradients during the rifting phases have been documented in the Sierras Pampeanas . Therefore, we reconstructed an exhumation history considering that the documented Mesozoic reheating events were mostly related to increases in the geothermal gradient during Mesozoic rifting, which implies that during these events, the paleodepths of this rock body remained roughly the same, as opposed to being buried.
We used the cold sample expected path from the MSM to extract several time-temperature predictions for the Paleozoic to Miocene thermal evolution of the basement of the Cuevas range ( Figure S6). This path requires the basement of the paleosurface to remain at temperatures below 120°C from the Early Paleozoic to the Triassic. In order to fit the negative correlation between the ZHe ages and eU, the cold sample expected path requires reheating between 160 and 200°C at 160 Ma. Cooling below 65°C took place~160 Ma, as the AFT data suggest. Reheating between 65 and 80°C happened~100 Ma. During the Miocene, the cold sample was reheated at temperatures between 55 and 85°C; this reheating was required to partially reset the AHe data and shorten the track length distributions.

10.1029/2019TC005887
Tectonics Using the cold sample expected path from the MSM and interpreting Mesozoic reheating to be a result of changes in the geothermal gradients, we calculate two extreme scenarios for the Paleozoic to Miocene exhumation history of the Cuevas range basement (Figure 5c). This exhumation history delimits the maximum and minimum values of the geothermal gradient and paleodepths inferred from the MSM, without considering burial. Maximum paleodepths can be inferred using the maximum temperatures from the thermal history, minimum surface temperatures, and a minimum initial geothermal gradient. Therefore, we used a reasonable minimum Paleozoic geothermal gradient of 15°C/km and a minimum mean annual surface temperature of 10°C. Analogously, minimum paleodepths can be inferred using minimum temperatures, maximum geothermal gradients, and maximum surface temperatures. Therefore, we used the minimum temperature requirements of the cold sample expected path, a reasonable maximum Mesozoic geothermal gradient of 55°C/km, and a maximum mean annual surface temperature of 20°C. The details of the inputs and the calculations are described in Text S3, and the results are presented in Table S5. The results show that burial is not necessary to explain the Mesozoic reheating events.

Tectonic Evolution of the Cuevas Range and the Campo-Arenal Basin
The granitic basement of the Cuevas range can be correlated with granitic rocks that formed during the Famatinian orogeny (550 and 450 Ma) in the Sierras Pampeanas and the Precordillera (Figure 1) (Ramos, 2008). The expected MSM shows an initial cooling event at~540 Ma, which could be related to a combination of postmagmatic cooling and basement exhumation during the Famatinian orogeny (Figures 5a and 6a).
In order to stay below 120°C until 160 Ma and explain the 160 Ma reheating peak only by changes in the thermal gradient, the Sample C3 must have remained between~2.5 and~7 km depth between~500 and 200 Ma (Figures 5c and 6a). After 200 Ma, an increase in the geothermal gradient of 11 to 16°C/km was responsible for the initial Mesozoic reheating (Figures 5c and 6b). Between 160 and 150 Ma, 1.5 to 5 km of exhumation raised the basement of the Cuevas paleosurface to between 1 and 2 km depth (Figures 5c  and 6c). This initial increase of the Triassic-Jurassic geothermal gradient and the subsequent exhumation coincide with the documented onset of the magmatic rifting phases in the Sierras Pampeanas (Oncken et al., 2006;Schmidt et al., 1995;Viramonte et al., 1999). The fast exhumation of the basement hosting the future Cuevas paleosurface can be explained by unroofing driven by a normal fault in the western horst shoulder of the Salta rift basin (Figure 6b). A humid climate in the Late Jurassic may have been the precursor to efficient erosion during this fast exhumation, creating a relatively low relief postrift landscape.
Between 120 and 100 Ma, the Salta Rift basin and associated volcanism (Bossi, 1969;Cristallini et al., 1997;González et al., 2000;Iaffa et al., 2011) were characterized by a high geothermal gradient between~33 and 55°C /km (Figures 5c and 6c). The increase of the gradient was responsible for reheating the basement of the Cuevas range above 65°C ( Figure 5b); a similar high geothermal gradient was documented farther south in the Velasco range during the Late Cretaceous . This reconstruction suggests that the Cuevas range experienced slow erosion rates (0.006 to 0.030 km/Ma) between 140 and 25 Ma (Figure 6d). Similar slow rates (0.010-0.026 km/Ma) have been documented from Mesozoic and Cenozoic paleosurfaces to the south of the Sierras Pampeanas Jordan et al., 1989).
The combination of low relief, the documented humid climate (Andrews et al., 2017;Margarita Do Campo et al., 2018), and the absence of deformation and sedimentation may have allowed the low erosion rates documented in the Cuevas range during the Cretaceous and the Paleogene. These conditions were ideal for the development of etchplains between weathered and fresh rocks. We interpret the Cuevas paleosurface to be a remnant of these etchplains, exposed after the removal of weathered rocks (Figure 6d). Modern low relief stable cratons under subtropical or tropical climates with low erosion rates (<0.01 km/Ma) (Chardon et al., 2018;Fairbridge & Finkl, 1980;Regard et al., 2016) could be a modern analog of the Late Cretaceous to Paleogene Sierras Pampeanas.
Approximately 1.4 km of Paleogene to Miocene sediments lay directly on top of the Cuevas paleosurface ( Figure 2b) (Bossi & Muruaga, 2009;Muruaga, 2001). The expected MSM model suggests differential reheating along the paleosurface between 50 and 85°C for the cold sample (C3). The hot sample (H2) experienced a more prolonged residence at higher temperatures between 80 and 95°C (Figure 5b). The temperature offset exhibits maximum values between 40 and 60°C at~10 Ma (Figures 5b and S4). As we already discussed, this 10.1029/2019TC005887 Tectonics differential reheating results from one or more magmatic heat sources to the west of the outcropping paleosurface (Figure 6e). This magmatism was related to the 12 and 5 Ma activity of the Farallon Negro volcanic complex. Moreover, several 9.5 Ma subvolcanic units intrude the Miocene sedimentary rocks deposited on top of the Cuevas paleosurface; this subvolcanic magmatism could also be the expression of the heat sources responsible for the interpreted differential reheating (Bossi & Muruaga, 2009;Mortimer et al., 2007;Muruaga, 2001). Thermal numerical models have shown how magmatic bodies can create a pattern of fully to partially reset low-temperature thermochronologic ages, as a function of the size and position of the magmatic heat source (Murray et al., 2018). The inversion of the Cretaceous and Miocene basins in the northern Sierras Pampeanas was the result of the eastward propagation of deformation related to the Andean orogeny (Mortimer et al., 2007;Pearson et al., 2012). The expected MSM shows a final fast cooling event in the Pliocene; this event is the result of the unroofing of the unconsolidated Miocene sedimentary cover (Figure 6d). The exhumation event was coeval with the deposition of conglomeratic strata of the Corral Quemado Fm. in the Campo Arenal basin after 3 Ma (Mortimer et al., 2007;Muruaga, 2001).

Implication for the Sierras Pampeanas Tectonic Evolution
The cold sample in our MSM exhibits the minimum Miocene reheating values along the paleosurface and thus the minimum geothermal gradients, which are calculated between 25 and 54°C/km (Figure 5c). Conversely, the hot sample requires minimum geothermal gradients of 42°C/km during the Miocene. These values agree more with the values presented by Stevens  and are higher than the values suggested by Dávila and Carter (2013) (Figure 6d). In our model, the absence of Cenozoic AHe and AFT ages is explained by the thin sedimentary cover, the lack of basement erosion, and the increase of the AHe system retentivity due to the accumulation of radiation damage during long residence at low temperatures (<120°C).
Our model suggests that no significant relief developed since the Famatinian orogeny until Pliocene time. In contrast, based on low-temperature thermal modeling, Permo-Triassic relief development has been suggested for the region south of the Sierras Pampeanas (Bense et al., , 2017. In the models presented in Bense et al. (2013) is not clear how the radiation damage was treated for the ZHe and the AHe systems. Moreover, these models do not account for possible partial Miocene reheating despite the documented Miocene-Pliocene strata in the Sierras Pampeanas (Davila et al., 2012;Schmidt et al., 1995;. Radiation damage controls and Miocene reheating can significantly modify the thermal modeling results, as has been shown here and in previous contributions . Additionally, the model presented in Bense et al. (2013) implies that the Pampean paleosurfaces formed after the development of the relief without explaining the formation mechanisms. Our tectonic model suggests that these paleosurfaces are the record of low relief humid landscapes prior to the development of the Sierras Pampeanas relief.
We propose that the Sierras Pampeanas paleosurfaces record the superposition of humid and dry climatic cycles (Rabassa et al., 2010;Rabassa & Ollier, 2014) and multiple rifting phases, during periods of relatively low relief. During each humid period, an associated etchplain may have formed. Subsequently, each Mesozoic rifting phase exhumed fresh basement rocks, where new etchplain levels were developed. However, more refined time constraints on these paleosurfaces are necessary to reconstruct this landscape response to these combined tectonic and climate cycles.
The prolonged tectonic evolution of the Sierras Pampeanas presents challenges for the modeling and interpretation of low-temperature thermochronological data. The long residence at low temperatures and the slow cooling rates of the Sierras Pampeanas basement rocks allowed ample time for radiation damage, causing dispersed U-Th/He ages. Additionally, the relatively thin Cenozoic sedimentary cover (1-3 km), the increase in retentivity of the AHe system, and the absence of significant basement erosion (as suggested by the paleosurfaces) during the Cenozoic explain the lack of fully reset AHe and AFT ages. These conditions complicate the interpretation of the exhumation phases and the calculation of the paleo-geothermal gradients of this province. Conversely, in this contribution, we have shown how the dispersed data and the preservation of ancient landscapes present an opportunity to refine the long-term evolution in the Central Andean province.

Conclusions
We exploit intersample and intrasample dispersions to reconstruct a refined thermal and tectonic evolution of the Cuevas paleosurface. Ancient landscapes offer an independent constraint for interpreting lowtemperature thermochronology and constraining thermal histories. SSMs based on two to four AHe aliquots and AFT data can be an acceptable approach to determining the thermal history of a sample. However, the MSM profits from a more robust and disperse data set to refine the thermal history of the paleosurface. Therefore, we can specify the temporal, tectonic, and thermal evolution of the Cuevas paleosurface, constrain the time and the erosion rates for the formation of this paleosurface, and quantify the Miocene thermal evolution of the Campo-Arenal basin.
The Cuevas basement experienced a cooling event at~540, interpreted as postmagmatic cooling and exhumation during the Famatinian orogeny. During the Mesozoic, geothermal gradients increased due to Mesozoic rifting, which reheated the rocks of the Cuevas range at~160 and 100 Ma. The basement of the Cuevas range experienced~1.5 to 5 km of normal fault footwall exhumation between 150 and 140 Ma, placing it between~1 and 2 km deep with respect to the surface. This event generated a low relief postrift landscape which was exposed to a humid climate, leading to slow erosion rates between~0.006 and 0.030 km/Ma. These conditions facilitated the development of an etchplain on top of the Cuevas range. This surface was reheated beneath~1.4 km of sedimentary cover during the Miocene; Miocene magmatic heat sources locally increased the amount of reheating up to 95°C. Finally, this paleosurface was tilted, exhumed, and uplifted in the Pliocene. The Sierras Pampeanas paleosurfaces may be the result of the alternation of humid and dry climates in a low relief landscape associated with anorogenic and active rifting periods.
The low erosion rates and the extended periods of tectonic quiescence responsible for the formation of these ancient landscapes offer ideal conditions for the accumulation of radiation damage in the AHe and ZHe systems. Therefore, eU and grain-size trends in the AHe system can be a powerful tool to refine the long-term thermal evolution of these ancient, slowly forming landscapes.