Distribution of Temperature and Strength in the Central Andean Lithosphere and Its Relationship to Seismicity and Active Deformation

We present three‐dimensional (3D) models of the present‐day steady‐state conductive thermal field and strength distribution in the lithosphere beneath the Central Andes. Our primary objective was to investigate the influence that the structure of the Central Andean lithosphere has on its thermal and rheological state, and the relationship between the latter and the active deformation in the region. We used our previous data‐driven and gravity‐constrained 3D density model as starting point for the calculations. We first assigned lithology‐derived thermal and rheological properties to the different divisions of the density model and defined temperature boundary conditions. We then calculated the 3D steady‐state conductive thermal field and the maximum differential stresses for both brittle and ductile behaviors. We find that the thickness and composition of the crust are the main factors affecting the modeled thermal field, and consequently also the strength distribution. The orogen is characterized by a thick felsic crust with elevated temperatures and a low integrated strength, whereas the foreland and forearc are underlain by a more mafic and thinner crust with lower temperatures and a higher integrated strength. We find that most of the intraplate deformation coincides spatially with the steepest strength gradients and suggest that the high potential energy of the orogen together with the presence of rheological lateral heterogeneities produce high compressional stresses and strong strain localization along the margins of the orogen. We interpret earthquakes within the modeled ductile field to be related to the weakening effect of long‐lived faults and/or the presence of seismic asperities.

There has been debate over recent decades regarding the rheological stratification of the lithosphere and its correlation with the thickness of the seismogenic layer T s (Burov, 2010;Chen et al., 2012;Jackson et al., 2008). The two main models of possible rheological profiles (also known as yield strength envelopes or YSEs) for the lithosphere that have been published are commonly known as the jelly sandwich model and the crème brûlée model; the former is characterized by the presence of a weak lower crust between a strong upper crust and a strong lithospheric mantle, while the latter consists of a strong upper crust with a weak lower crust and a weak lithospheric mantle (e.g., Burov & Watts, 2006;Jackson, 2002). Comparative analyses of intraplate seismicity and rheological stratifications have revealed a correlation between frequency-depth distributions of earthquakes and lithospheric strength profiles, with peak seismicity occurring close to brittle-ductile transitions and cut-off depths of 10−20 km in most regions (e.g., Doser & Kanamori, 1986;Ito, 1990). However, the occurrence of intermediate to deep crustal earthquakes beneath the brittle-ductile transition in some regions has cast doubt over this correlation, leading some authors to propose variations in slip behavior with depth and different mechanisms for seismogenesis (e.g., Deichmann, 1992;Foster & Jackson, 1998;Hobbs et al., 1986;Jackson & Blenkinsop, 1993;Lamontagne & Ranalli, 1996;Petley-Ragan et al., 2019;Prieto et al., 2017;Scholz, 1988;Tse & Rice, 1986).
The variability in the modeled lithospheric rheological profiles and their unclear correlation with respects to s T require the separated investigation of each region in order to understand its present deformation processes. The area of the southern Central Andes of South America is a particularly interesting region in which to explore the different hypotheses regarding the driving factors of mountain building (e.g., Allmendinger et al., 1997;DeCelles, Ducea, et al., 2015;Jordan et al., 1983;Oncken, Chong, et al., 2006;Ramos, 2008;Ramos et al., 2002). Although a large variety of models have previously been used to investigate the structures, seismicity, thermal field, and lithospheric composition of this region (e.g., Babeyko et al., 2002Babeyko et al., , 2006Eichelberger et al., 2015;Gerbault et al., 2003;Hindle et al., 2005;Metcalf & Kapp, 2015;Mulcahy et al., 2014;Ouimet & Cook, 2010;Prezzi et al., 2009;Salomon, 2018;Tassara et al., 2006;Yang & Liu, 2003), none have been used to either fully investigate rheology and thermal field as a function of lithology distribution, or to analyze the relationship of rheology and thermal field with seismicity.
Research into crustal seismicity and neotectonic structures in the southern Central Andes has shown that most of the recent activity has been concentrated along the margins of the orogen (e.g., Devlin et al., 2012;Oncken, Hindle, et al., 2006;Wigger, 1988). In a previous publication we suggested the possibility of a link between the localization of active deformation and the strength of the lithosphere as a subject that required further investigation (Ibarra et al., 2019). In this contribution, we therefore set out to investigate the rheological state of the southern Central Andes, how it relates to the lithospheric composition of the region, and whether or not the heterogeneous distribution of seismicity and neotectonic structures is controlled by regional variations in the structure, temperature, and rheology of the lithosphere.
To that end we present 3D lithology-derived numerical models of the present-day steady-state conductive thermal field and lithospheric strength. We have further improved our previous data-driven and gravity-constrained 3D density model of the region, by introducing thermal and rheological properties based on published average values for each lithology (Afonso & Ranalli, 2004;Čermák & Rybach, 1982;Goetze & Evans, 1979;Ranalli & Murphy, 1987;Vilà et al., 2010;Wilks & Carter, 1990). We used the seismicity of the region and the neotectonic structures compiled from global and local databases as the active deformation against which to compare our results.

Geologic Setting and Major Geophysical Features
The region of the southern Central Andes is part of a subduction-related orogen on the western margin of the South American plate (Figure 1). Subduction processes in South America have been almost continuous since the early Paleozoic, causing alternating periods of compression and extension (e.g., Bock et al., 2000;Coira et al., 2009;Maloney et al., 2013;Ramos, 2009Ramos, , 2010. The compressive Andean tectonic processes that have resulted in the present-day morphotectonic provinces, however, did not start until the Middle Cretaceous (e.g., Amilibia et al., 2008;Bascuñán et al., 2016;Rossel et al., 2013). It has been inferred that the acceleration of the spreading rate of the southern Mid-Atlantic Ridge during the Cretaceous produced an increase of the trench-ward velocity of the South American plate over the retreat velocity of the subducting Nazca plate, resulting in a compressive stress field of the upper plate (e.g., Pardo-Casas & Molnar, 1987;Somoza, 1998). However, the N-S segmentation of the Andes and the spatio-temporal variations in Cenozoic deformation have been attributed to either the shallowing-steepening cycles of slab segments and changes in the state of interplate coupling (e.g., Horton, 2018;Jordan et al., 1983;Martinod et al., 2010), or the differential reactivation of basement heterogeneities of the South American crust (e.g., Hongn et al., 2010Hongn et al., , 2007Kley et al., 1999;Kley & Monaldi, 2002;Pearson et al., 2013).
IBARRA ET AL.   Burns et al., 2015, and Kay et al., 2010, regions with anomalous low crustal velocity and resistivity (dashed black lines-after Beck et al., 2015, andBianchi et al., 2013), and a region with high seismic velocity (dashed light blue lineafter Schurr & Rietbrock, 2004). The large red rectangle in broken lines outlines the region that was modeled and which is shown in Figure 2. APVC: Altiplano-Puna Volcanic Complex; APMB: Altiplano-Puna Magma Body; SPMB: Southern Puna Magma Body; AB: Atacama Block.
The region of the southern Central Andes that we modeled consist of (from west to east) forearc, volcanic arc, and backarc areas, which are in turn sub-divided into different morphotectonic/geologic provinces or units (i.e., regions characterized by particular geomorphologic features, structural style, and stratigraphy) on the basis of their tectonic evolution Ramos, 1999). The forearc comprises the Coastal Cordillera (CC), the Longitudinal Valley (LV), and the Chilean Precordillera (PC), and the volcanic arc is represented by the Western Cordillera (WC). The backarc includes the Andean Plateau (or Altiplano-Puna Plateau; AP and PN, respectively), the Eastern Cordillera (EC), the Subandean Ranges (SR), the Santa Bárbara System (SB), the Pampean Ranges (PR), and the undeformed Chaco-Parana foreland basin (FB; Figure 2).
The Coastal Cordillera is a prominent morphological feature of the forearc that primarily consists of Jurassic to Lower Cretaceous basic to andesitic volcanic and plutonic rocks associated with a thinned continental crust (González et al., 2003;Lucassen et al., 2006;Rossel et al., 2013). To the east, the Chilean Precordillera is a thick-skinned thrust belt formed by Paleozoic igneous-metamorphic rocks, Mesozoic and Tertiary sedimentary and volcanic rocks, and Late Cretaceous to Paleogene plutons (Amilibia et al., 2008;Lucassen et al., 2001;Mpodozis & Ramos, 1989;Scheuber et al., 1994). The Western Cordillera corresponds to the present-day volcanic arc and consists of Mesozoic sedimentary rocks and voluminous Tertiary volcanic rocks and plutons (Scheuber & Reutter, 1992). Farther east, the Andean plateau (Altiplano-Puna) is an intraorogenic and internally drained, high-elevation region comprised of numerous fault-bounded basins developed on Late Neoproterozoic to Paleozoic igneous and metamorphic rocks, filled with Cretaceous and mostly Cenozoic clastic sedimentary rocks, evaporites, and volcanics, locally reaching thicknesses of IBARRA ET AL.   Jordan et al., 1983, andRamos 1999), the measured surface heat flow (colored triangles-after Hamza & Muñoz, 1996, Henry & Pollack, 1988, and Springer & Förster, 1998 and the focal mechanisms of crustal earthquakes. Focal mechanisms shown in orange, red, and black were taken from Devlin et al. (2012), Mulcahy et al. (2014), and the Global Centroid Moment Tensor Catalog (Dziewonski et al., 1981;Ekström et al., 2012), respectively. The sizes of the focal mechanisms correspond to the magnitude of the earthquakes, but note that the scale of the red mechanism symbols has been doubled for the sake of clarity. The boundaries between the morphotectonic provinces are generally transitional; the solid black lines thus do not represent any geological feature in nature. CC: Coastal Cordillera; LV: Longitudinal Valley; PC: Chilean Precordillera; WC: Western Cordillera; AP: Altiplano; PN: Puna; EC: Eastern Cordillera; SR: Subandean Ranges; SB: Santa Bárbara System; PR: Pampean Ranges; FB: Chaco-Parana foreland basin. >6 km (Alonso et al., 1991;Jordan & Alonso, 1987;Siks & Horton, 2011). The Eastern Cordillera is a thickskinned thrust belt mainly composed of Late Neoproterozoic to Paleozoic metamorphic rocks, covered by Cretaceous and Cenozoic sediments and volcanics (Mon & Salfity, 1995). In transition to the foreland, the Subandean Ranges are a thin-skinned thrust belt involving Paleozoic, Mesozoic, and Cenozoic sedimentary rocks (Dunn et al., 1995;Echavarria et al., 2003;Mingramm et al., 1979). Farther south the Santa Bárbara System is a compressionally inverted, tectonically active Cretaceous rift province that exposes Paleozoic metasedimentary rocks and overlying Cretaceous and Cenozoic sedimentary and volcanic rocks (Kley & Monaldi, 2002;Marquillas et al., 2005;Mon & Salfity, 1995); and the Pampean Ranges are tectonically active reverse-fault bounded basement blocks mainly composed of Late Proterozoic to Paleozoic metamorphic and igneous rocks, often associated with asymmetric uplift along reactivated Paleozoic shear zones (González-Bonorino, 1950;Jordan & Allmendinger, 1986;Ramos et al., 2002;Strecker et al., 1989;Toselli et al., 1978).
A number of previous investigations have revealed important W-E heterogeneities within the lithosphere (e.g., Lucassen et al., 2001;Meeβen et al., 2018;Prezzi et al., 2009). Most of the forearc is characterized by a relatively thin mafic crust beneath the Coastal Cordillera, thickening toward the east; this layer exhibits high seismic velocity and low attenuation (e.g., Assumpção et al., 2013;Rossel et al., 2013;Schurr & Rietbrock, 2004). In contrast, the adjacent volcanic arc (Western Cordillera) and the orogenic backarc are characterized by a thick felsic crust with low seismic velocity and a thin lithospheric mantle (e.g., Beck & Zandt, 2002;Koulakov et al., 2006;Yuan et al., 2002). Farther east again, the foreland is characterized by a thinner crust with a felsic upper layer, a mafic lower crust, and slightly higher seismic velocity than within the orogen (e.g., Feng et al., 2007;Rosa et al., 2016).
Three major anomalies have been identified within the crust by means of geophysical imaging ( Figure 1): These are the Altiplano-Puna Magma Body (APMB; Chmielowski et al., 1999), the Southern Puna Magma Body (SPMB; Bianchi et al., 2013), and the Atacama Block (Schurr & Rietbrock, 2004). The APMB lies between 10 and 40 km beneath the plateau area and between −20° and −23°, and is characterized by low resistivity, low seismic velocity, a high Vp/Vs ratio, and high attenuation (Brasse et al., 2002;Comeau et al., 2015Comeau et al., , 2016Díaz et al., 2012;Koulakov et al., 2006;Ward et al., 2013). This anomaly is interpreted to represent a region that is undergoing metamorphism and partial melting (e.g., Beck & Zandt, 2002;Zandt et al., 2003). The SPMB has been defined by Bianchi et al. (2013), who detected a low-velocity zone in the southern Puna crust between −26° and −27.5° (using P-wave tomography), which they suggested could be due to another mid-crustal region with active partial melting. The Atacama Block is located within the Chilean Precordillera area, beneath the Atacama Basin; it is characterized by high seismic velocity and low attenuation (Schurr & Rietbrock, 2004).
The heat flow distribution correlates fairly well with the aforementioned W-E heterogeneities of the lithosphere, with low to intermediate heat flow in the forearc (20-30 mW m −2 ; Figure 2) and foreland (50-70 mW m −2 ; Figure 2) areas, and increased heat flow values within the orogen (over 100 mW m −2 ; Figure 2). The surface heat flow is extremely high (up to 320 mW m −2 ) throughout the plateau area, coinciding with the volcanic arc and the two low-velocity anomalies (Hamza & Muñoz, 1996;Henry & Pollack, 1988;Springer & Förster, 1998).
Information on seismicity over recent decades, as well as on sites of active deformation, has been compiled and analyzed in various projects, revealing that most of the activity has been focused on the Coastal Cordillera, the Chilean Precordillera, the Eastern Cordillera, the Subandean Ranges, and the broken foreland of the Santa Bárbara System and the Pampean Ranges (e.g., projects PUDEL, PUNA '97, PISCO '94 and ANCORP '96; Armijo et al., 2015;Arnous et al., 2020;Costa et al., 2018;García et al., 2019;Graeber, 1997;Heit, Sodoudi, et al., 2007;Meigs & Nabelek, 2010;Rietbrock et al., 1997;Santibáñez et al., 2019;Schurr et al., 1999;Siame et al., 2015;Strecker et al., 1989;Weiss et al., 2015, Figure 3). Although tectonic structures indicating active normal and strike-slip faulting in the Andean Plateau exist (e.g., Allmendinger et al., 1989;Montero-López et al., 2010Schoenbohm & Strecker, 2009), the region is generally characterized by a low level of seismicity compared to the other Andean sectors (Mulcahy et al., 2014;Oncken, Hindle, et al., 2006). Generally, first-order structures and deformation zones are parallel to the N-S orientation of the orogen or strike NE-SW (Figure 3), with a dominance of reverse faulting and strike-slip components as indicated by focal mechanisms (Figure 2) and fault-slip analysis (e.g., Riller & Oncken, 2003). The Coastal Cordillera has recorded the deepest events (up to depths greater than 50 km; Figure 3), and also the largest range in depths because the earthquakes occur at the subduction interface. Most of the seismicity to the east of the plateau area (i.e., the Santa Bárbara System and the Pampean Ranges) is restricted to depths between 10 and 30 km (Figure 3). Within the Puna area the activity is much shallower (within the upper 15 km; Figure 3) with coexisting compression, trans-tension, trans-pression and extension, as indicated by focal mechanisms (see symbols in Figure 2) and geological structures (Allmendinger et al., 1989;Montero-López et al., 2010;Schoenbohm & Strecker, 2009).

Methods and Data
For this research we modeled the present-day thermo-mechanical state of the lithosphere in the Central Andes. The methodology involved calculating the 3D steady-state conductive thermal field and computing the maximum differential stresses for both brittle and ductile deformation mechanisms, taking into account lithology-dependent thermal and mechanical properties. All calculations were performed on the basis of our previously published data-derived and gravity-constrained 3D density model of the region (   Armijo & Thiele, 1990, Armijo et al., 2015, Oncken, Hindle, et al., 2006, Riller & Oncken, 2003, Santibáñez et al., 2019, and Weiss et al., 2015. Dashed light blue lines outline the boundaries between morphotectonic units as in Figure 2. the western part of deformed foreland; while the Eastern Domain corresponds to the foreland and eastern parts of the Subandean Ranges, the Santa Bárbara System, and the Pampean Ranges. The other two domains are the Atacama Block and the APMB-SPMB, which are restricted within the Central Domain, and beneath the Atacama Basin and the Andean Plateau, respectively. Our model consists of (i) a continental crust divided into five domains (Western Domain, Central Domain, Eastern Domain, Atacama Block and APMB-SPMB), each of which are sub-divided vertically; (ii) a single sedimentary layer; (iii) an oceanic crust, and (iv) the mantle. All divisions are based on seismic studies and previous models of the region, with uncertainties inherent to the gravity modeling approach, particularly along the vertical direction. The model box extends between longitudes −71° and −60°, and latitudes −29° and −19°, and is restricted to a maximum depth of 220 km; however, the vertical extent of the model in this study was set to a depth of 100 km in order to include only the lithospheric part of the mantle (Figure 4). By using this model, we can take into account the lateral and vertical variations in the radiogenic heat production, the thermal conductivity, and the rheological properties of the crust.
The results of the models were compared to and analyzed with existing data on neotectonic deformation and seismic activity. The location of first-order neotectonic structures and deformation zones ( Figure 3) was compiled from previous publications (Armijo & Thiele, 1990;Armijo et al., 2015;García et al., 2019;Oncken, Hindle, et al., 2006;Riller & Oncken, 2003;Santibáñez et al., 2019;Weiss et al., 2015), while the hypocenters of seismic events ( Figure 3) were taken from the ISC-EHB Bulletin (International Seismological Centre, 2020) and published data from a temporary local seismic network (PUDEL; Mulcahy et al., 2014). The ISC-EHB Bulletin contains teleseismically well-constrained seismic events with magnitudes >3.75 from the ISC Bulletin for the period 1964−2016, which were relocated using the EHB algorithm (Engdahl et al., 1998;Weston et al., 2018) to minimize errors in location (particularly depth) due to an inferred 3D Earth structure. Based on depth quality, the database is divided into three categories − L1, L2, and L3 − with free depth errors <5 km, between 5−15 km, and >15 km, respectively. The data published in Mulcahy et al. (2014) includes earthquakes recorded by the Southern Puna Seismic Array between March 2008 and September 2009 that met good quality criteria, with maximum calculated magnitudes of 3.2, and vertical and horizontal errors generally smaller than 15 and 10 km (for 80% of the data), respectively.

Thermal Model
We derived the temperature distribution in the Central Andean lithosphere by calculating the conductive thermal field for steady-state conditions, based on the assumption that conduction is the dominant mechanism for heat transport within the lithosphere. The conductive heat equation, also known as Fourier's law, for steady-state conditions is as follows: where ∇ is the Nabla operator, S is the radiogenic heat production, T is the absolute temperature, and λ b is the bulk thermal conductivity.
To solve the equation in 3D for a thermally equilibrated system we used GOLEM, a MOOSE-based application . GOLEM is a numerical simulator for modeling coupled thermo-hydro-mechanical processes (for further information on the mathematical and computational aspects see . The simulation requires the assignment of thermal properties (thermal conductivity and radiogenic heat production) to every division or unit in the density model, as well as the definition of upper and lower thermal boundary conditions.
Each unit was assigned fixed bulk thermal conductivity and radiogenic heat production values on the basis of its prevailing lithology and published laboratory measurements of rock samples with comparable lithologies (Table 1; Čermák & Rybach, 1982;Vilà et al., 2010). We set the annual mean surface temperature distribution derived from global climatological models as the upper boundary condition (Banzon et al., 2016;New et al., 2002;Reynolds et al., 2007, Figure 5a), and the temperature distribution at 100 km depth within the mantle as the lower boundary condition (Figure 5b).
We decided to crop the density model to 100 km depth in order to ensure that our thermal model included only the lithospheric part of the mantle (of both, the upper and the subducting plates), in which the initial assumption of conduction as the main mechanism for heat transport was likely to be valid. We based our decision for this particular depth on receiver function studies across the region that imaged the lithosphere-asthenosphere boundary at depths as shallow as 100 km (Heit et al., 2008;Heit, Sodoudi, et al., 2007).
In order to maintain coherence between the density and thermal models, the temperature distribution at the lower boundary condition was derived by converting seismic velocities from the same global tomographic IBARRA ET AL.  Ranalli and Murphy (1987). e Wilks and Carter (1990). f Afonso and Ranalli (2004). g Goetze and Evans (1979). * brittle behavior was imposed on the sediments.  model that Ibarra et al. (2019) used to calculate the densities in the mantle. The method used (Meeßen, 2017) is a modified version of the approach developed by Goes et al. (2000) to convert seismic velocities into temperature and density for a given rock composition, taking into account the effects of anharmonicity and anelasticity. We used the global shear-wave tomography model of Schaeffer and Lebedev (2013) because of its higher resolution compared to other models (0.5° in latitude and longitude, and 25 km in the z direction).
In order to take into account the rock compositions in this conversion we relied on recorded compositions of mantle xenoliths of Cretaceous-Cenozoic age from within the region . Further details of the method can be found in Text S1 of the Supplementary Information.

Rheological Model
To assess the strength of the lithosphere throughout the modeling region we calculated the yield strength for both brittle and ductile deformation mechanisms. The calculations were based on our density and thermal models, since the structural interfaces defining different rock units within the lithosphere and the temperature distribution were pre-requisites for the rheological calculations.
The strength of rocks is defined as the maximum differential stress (Δσ max ) they can resist without experiencing permanent deformation (either brittle or ductile), and is highly dependent on temperature and pressure as well as on the rock composition, strain rate and state of existing deformation (Ranalli, 1995(Ranalli, , 1997: where σ 3 and σ 1 are the minimum and maximum principal stresses, respectively. As a general rule, rocks experience brittle deformation at shallow depths and ductile deformation at greater depths, depending mainly on the temperature distribution; the prevailing mechanism of deformation will be the one that requires the least differential stress at a given depth (Goetze & Evans, 1979). Brittle behavior is governed by the Coulomb-Navier failure criterion, described by an empirically determined temperature-independent law (Byerlee, 1968), that is usually referred to as Byerlee's law: where Δσ b is the brittle yield strength, ρ b is the bulk density, g is the acceleration due to gravity (g = 9.81 m s −2 ), z is the depth below topography, f p is the pore fluid factor (ratio of pore fluid pressure to lithostatic pressure; a typical hydrostatic value is f p = 0.36; Ranalli, 1995), and f f is a coefficient that depends on the friction coefficient of rocks and the fault type. Since the Andean tectonic regime is overall compressive, a f f = 3.00, representative of thrust faulting and a general friction coefficient of 0.75, yield a good approximation for deformation in general (Ranalli, 1995).
Ductile deformation is controlled by solid-state creep mechanisms (dislocation or glide). In the crust and upper mantle, dislocation creep is the dominant deformation mechanism (described by the power law in Equation 4); however, when differential stresses of dislocation creep are greater than 200 MPa in the mantle, Dorn's law describing dislocation glide for olivine (Equation 5) yields a better approximation for deformation of mantle rocks (Goetze, 1978;Kirby & Kronenberg, 1987): where Δσ d is the ductile yield strength,  is the reference strain rate (= 6.00E−15 s −1 at the Andean margin; Kreemer et al., 2014), A p is the pre-exponential scaling factor, n is the power-law exponent, Q p is the power-law activation energy, σ D is the Dorn's law stress (σ D = 8.5E9 Pa), Q D is the Dorn's law activation energy (Q D = 535 kJ mol −1 ), A D is the Dorn's law strain rate (A D = 5.7E11 s −1 ), R is the universal gas constant (R = 8.314 J K −1 mol −1 ) and T is the absolute temperature.
We used the code presented in Cacace and Scheck-Wenderoth (2016) to calculate the maximum differential stress for brittle and ductile behavior at each X-Y-Z position (Equations. 3, 4, and 5); the lowest of these two defines the yield strength at that particular point (Goetze & Evans, 1979, Equation 6). Integrating the yield strength over depth at each X-Y position for the entire lithosphere and the entire crust provides the integrated lithospheric strength and the integrated crustal strength, respectively (Equation 7): Integrated strength maps are useful for analyzing horizontal variations in rheology, while the establishment of yield strength envelopes, that is, plotting the yield strength against depth for a particular X-Y position, allows us to identify vertical heterogeneities.
To perform the calculations the code requires physical properties for each unit (Table 1) and the temperature distribution throughout the model. The physical properties were taken from published laboratory measurements on common rock types, corresponding to the prevailing lithology in each unit (Table 1, Afonso & Ranalli, 2004;Goetze & Evans, 1979;Ranalli & Murphy, 1987;Wilks & Carter, 1990). Since the properties reported for similar lithologies by different authors are not always consistent, we took care when selecting secondary creep parameters to ensure that felsic units were weaker than mafic ones. A detailed description of the selection procedure can be found in Text S2 of the Supplementary Information.

Thermal Model
Depth slices from the model are presented in Figures 6a-6e. The same general trend can be seen at all depths, with the orogen being hotter than the adjacent domains. Since the plotted temperature maps correspond to depth below the surface, we can exclude the possibility of this being a topographic effect. Temperatures increase from the Eastern Domain toward the Central Domain and then decrease again to the Western Domain (toward the trench). There are high temperatures throughout the Western Cordillera and the Altiplano-Puna plateau, particularly beneath the southern Puna (spatially correlated with the SPMB), and beneath the boundary area between the Altiplano and the northern Puna (spatially correlated with the APMB). The high-temperature anomalies become more conspicuous with depth toward the middle crust, but then become progressively less well defined within the lower crust. The southern temperature anomaly covers a smaller area than the northern one, as well as being less distinct and more irregular.
At depths of 5-10 km in the upper crust the temperature range of the anomalies corresponding to the APMB and the SPMB is 325°C−475°C, compared to average temperatures in the Central Domain of 275°C-400°C (Figures 6a and 6b). In the middle crust, at depths of 25-40 km, the average temperature for the same region is between 675°C to 875°C, with the anomalies reaching up to 780°C-960°C (Figures 6c and 6d). The average temperature in the lower crust of the Central Domain is close to 1,000°C at 60 km depth (Figure 6e). The average temperatures in the upper crust (5−10 km) decrease toward the Western and Eastern Domains to 160°C−260°C (Figures 6a and 6b). In the Western Domain, temperatures increase with depth to up to 450°C in the middle crust and 600°C in the lower crust (Figures 6c and 6d), while in the Eastern Domain, there is a progressive increase in average temperature with depth, reaching 600°C at the Moho (Figure 6f).
Although surface heat flow is not a highly reliable parameter to compare (Hamza & Muñoz, 1996;Scheck-Wenderoth & Maystrenko, 2013), it does provide information on the general trend of the thermal field throughout the region. A comparison of surface heat flow derived from the model with that measured in the field (Hamza & Muñoz, 1996;Henry & Pollack, 1988;Springer & Förster, 1998

Rheological Model
In Figure 8 we present our modeled integrated strength for the lithosphere ( Figure 8a) and crust (Figure 8b), together with the percentage of the crustal contribution to the integrated lithospheric strength (Figure 8c). Both strength maps show the same pattern, with a relatively low integrated strength in the Central Domain increasing toward the Western Domain and the Eastern Domain. The lowest integrated strengths for both the lithosphere and the crust lie beneath the Altiplano-Puna plateau where most of the strength resides in the crust, as shown by the high proportion of crustal strength (more than 80%, Figure 8c). Although the Western Domain, the Eastern Domain, and the Atacama Block present areas with high integrated lithospheric strength, the highest integrated crustal strength is modeled in the Western Domain and the Atacama Block. The sharp rectangular shape of the latter is a model artifact due to the coarse resolution and the structured type of the model grid. Figure 8d shows the horizontal gradient of the integrated lithospheric strength, together with epicenters of crustal seismicity and neotectonic structures in the region (Figure 3). Most of the epicenters and active structures coincide with regions in which steep gradients in strength are modeled, that is, they are located at the transitions between weak and strong domains, particularly in the Pampean Ranges, the Santa Bárbara System, the Subandean Ranges, the Chilean Precordillera, and the eastern Puna.
In Figure 9 we present plots of yield strength envelopes for those locations that are structurally representative of the morphotectonic units within the modeled area (see Figure 8d for the locations). In the Eastern Domain the mantle is relatively strong, accounting for an important part of the total strength and conforming to the jelly sandwich model for the lithosphere (Figures 9a and 9b & 9c). In the Subandean Ranges and the Santa Bárbara System, the lower crust is relatively weak and decoupled from the mantle (Figures 9a and 9b), in contrast to the Pampean Ranges where there is a strong lower crust and crust-mantle coupling (Figure 9c). The yield strength envelope for the Puna plateau reveals that the strength there resides mainly in the upper crust, conforming to the crème brûlée model for the lithosphere (Figure 9d). The Atacama Block presents a strong upper crust, a very strong middle crust, and a weak mantle (Figure 9e), which is also regarded by Burov and Watts (2006) as conforming to the crème brûlée model.
In order to test how the configuration of the seismogenic layer relates to the modeled brittle-ductile transitions (BDTs), we compared the distribution of seismicity with the vertical distribution of strength in the model. For this we selected only those hypocenters with less than 5 km of vertical uncertainty and plotted depth-frequency distribution histograms for the morphotectonic units that contained most of the seismic activity (the Pampean Ranges, the Santa Bárbara System, and the Puna plateau; Figure 10). As depth proxies to compare against BDTs, we used the depth at which the frequency of hypocenters starts to decrease (D dec ) and the maximum depth of the hypocenters, equivalent to T s (D max ). These proxies were used previously by different authors who found correlations in other regions between those depths and BDTs (e.g., Chen & Molnar, 1983;Meissner & Strehlau, 1982;Sibson, 1982).
The BDTs in the crust beneath the Puna plateau, the Santa Bárbara System, and the Pampean Ranges lie at depths of 7−8 km, 16−17 km and 21−22 km, respectively (Figures 9b and 9c & 9d), and are thus shallower than D dec and D max in the corresponding morphotectonic unit (particularly, D max is consistently ∼15 km below the modeled BDTs; Figure 10). However, they follow the same trend as the hypocenter depth-frequency distribution, with greater depths toward the Pampean Ranges. Moreover, the discrepancies between the BDTs and D dec are between 1 and 5 km, which is within the uncertainty range for the depth of hypocenters.   Figure 11). Solid black lines outline the morphotectonic units (as in Figure 2).
We examined the sensitivity of our model with respect to rock-type rheology by using different rock types for the least well-constrained unit (the lower crust in the Eastern Domain). The similar trends in crustal and lithospheric strengths ( Figure S2 in the Supplementary Information) suggest that the correlation between the highest gradients of strength and the location of active deformation (Figure 8d) is robust. The correlation between brittle-ductile transitions and hypocenters in the Santa Bárbara System remains the same for each of the different rock types, because there is no brittle behavior in the lower crust ( Figure S3 Armijo & Thiele, 1990, Armijo et al., 2015, Oncken et al., 2006b, Riller & Oncken, 2003, Santibáñez et al., 2019, and Weiss et al., 2015 and epicenters of crustal seismic events (red circles-taken from the ISC-EHB Bulletin and Mulcahy et al., 2014). Numbers 1 to 5 mark the position of the yield strength envelopes shown in Figure 9 a to e, respectively. Dashed black lines outline the morphotectonic units (as in Figure 2).

Thermal Model
Our modeled thermal state of the lithosphere is based on the assumptions that heat is transported solely by conduction and that the system is in thermal equilibrium. Published research has demonstrated that heat conduction is the prevailing mechanism for heat transport within the lithosphere (e.g., Jaupart & Mareschal, 2007;McKenzie et al., 2005), and that steady-state conditions are valid for subduction systems old enough for the thermal structure in the mantle to have equilibrated (i.e., at least 200 m.y. old), which is the case for the Andes (e.g., Allmendinger et al., 1997;Hall, 2012;Isacks, 1988). However, the crust also exerts some control over the thermal field and the Andean crust is known to have thickened over the past 40 m.y.  may be too short for thermal equilibrium to have been achieved. There has also been extensive magmatic activity in the volcanic arc as well as in the backarc (particularly during the past 15 m.y.; Kay et al., 2010;Kay & Coira, 2009;Trumbull et al., 2006), associated with hydrothermal activity and exceptionally high surface heat flow, reflecting mechanisms of advective/convective heat transport in the uppermost crustal levels (Hamza & Muñoz, 1996;Henry & Pollack, 1988;Springer & Förster, 1998).
The advection of cooler crustal material to deeper levels through crustal stacking would initially reduce the temperature at those levels. However, an increase in temperature would be expected with time in response to high radiogenic heat production from the stacked upper crustal units and progressive thermal equilibration through heat conduction. Our steady-state model would therefore overestimate temperatures in the middle crust of the Central Domain if thermal equilibrium had yet to be achieved in the region.
Nonetheless, the results obtained from geothermobarometric investigations in the Puna (Burns et al., 2015;Schmitt et al., 2001) suggest otherwise. Geothermobarometry relies on different chemical equilibria in minerals to estimate the depth and temperature at which the magma from which the rocks derived was last in equilibrium. Crustal temperatures calculated for recent eruptive deposits (less than 2.5 m.y. old) on the Altiplano-Puna plateau indicate the presence of a deep magma reservoir at 17−21 km depth, with temperatures between 925°C−973°C (Burns et al., 2015;Schmitt et al., 2001). Moreover, the integrative analysis of seismic velocity and electric conductivity anomalies together with petrological constraints, suggests the presence of a crustal layer with andesitic partial melts, exhibiting temperatures between 930°C−970°C, as shallow as 18−19 km depth (Pritchard et al., 2018). These estimates exceed our modeled temperatures by at least 50°C, considering common uncertainties of 50°C in thermobarometric calculations and 25°C in our modeled temperatures; this suggests that our model does not overestimate the temperatures at this depth in the Central Domain (Puna region), which would be expected if the system were far from thermal equilibrium.
For the shallow upper crust of the Central Domain, the geothermobarometric calculations indicate a shallow magma reservoir at 4−8 km depth with temperatures of 770°C−890°C (e.g., Burns et al., 2015;Schmitt et al., 2001), while the comprehensive analysis of geophysical anomalies suggests the existence of shallow dacitic partial melts (top at 9 km) with a temperature of 660°C (Pritchard et al., 2018). It should be noted that our model may underestimate the local temperature by at least 200°C. Furthermore, the modeled IBARRA ET AL. values for surface heat flow within the Central Domain rarely exceed 90 mW m −2 , providing a good fit in the Altiplano, but not in the Puna and the volcanic arc, where most measurements are higher (up to 320 mW m −2 ; Figure 7).
Discrepancies between the conductive thermal model and the measured surface heat flow and estimated temperatures can likely be explained by convective and advective processes. Anomalously high heat flow values have been reported for locations with active volcanism and/or hydrothermal activity (Hamza & Muñoz, 1996;Henry & Pollack, 1988;Springer & Förster, 1998). Furthermore, low resistivity and velocity anomalies have been detected near the surface, suggesting the presence of shallow aquifers, hydrothermally altered rocks, and brines in the Andean plateau (Comeau et al, 2015(Comeau et al, , 2016Pritchard et al., 2018); within such hydrothermal systems, heat transport would be enhanced by processes of coupled fluid and heat circulation. Further deep in the crust, geochemical and geophysical investigations suggest that magma reservoirs at ∼8 km depth are being fed by the advection of hot material from deep reservoirs at ∼20 km depth (e.g., Delph et al., 2017;Kay et al., 2010;Pritchard et al., 2018;Schmitt et al., 2001). Such a feeding system could be at the origin of the elevated temperatures calculated for pre-eruptive conditions with respect to our modeled conductive thermal field.
There remains debate as to the origin of the large mid-crustal magma bodies. Published conceptual models can be divided into two opposing hypotheses with respect to the contribution of the mantle and the crust. One of the models proposes that the increased heat flow and magma production in the mantle triggered by delamination processes produced a MASH zone (melting, assimilation, storage, homogenization;Wörner et al., 1992) at the base of the crust that feeds the mid-crustal reservoirs and induces anatexis (e.g., Babeyko et al., 2002;Delph et al., 2017;de Silva & Kay, 2018;Kay et al., 2010;Perkins et al., 2016). This model is supported by the presence of mafic volcanic rocks in the plateau with mantle signature (e.g., Kay et al., 1994), the 1:1 ratio of mantle and crustal contribution inferred from the analysis of ignimbrites (Kay et al., 2010 and references therein), and the reported low-velocity zones in the lower crust beneath the APMB and at the base of the crust beneath the southern Puna (Delph et al., 2017;Kukarina et al., 2017).
However, other geochemical analyses of ignimbrites from the plateau suggest that even if the MASH zone were present, its contribution would be minor compared to the crust, given the estimated low percentage of mantle-derived magma (less than 30%;de Silva et al., 2006;Francis et al., 1989). Furthermore, the analysis of mafic samples suggests that lithospheric removal occurred as small-scale Rayleigh-Taylor type instabilities, with low volumes of melts associated with foundered lithospheric material and no asthenosphere melting (Drew et al., 2009;Murray et al., 2015). Considering mass balance calculations (Ducea & Barton, 2007) and the timing of geodynamic models for Rayleigh-Taylor instabilities and magma ascent (Polyansky et al., 2016;Wang & Currie, 2015), it is difficult to establish a clear causal relationship between lithospheric removal and the partially molten middle crust. To address these issues, the opposing model for evaluating the origin of the large mid-crustal magma bodies proposes that both lithospheric removal and crustal melting are coeval and a consequence of the accelerated thickening of the Andean crust (DeCelles et al., 2009, DeCelles, Zandt, et al., 2015. Crustal processes are more relevant in this model, and melting conditions are attained due to the increased radiogenic heat production of the thickened crust.
Although our steady-state model is unable to discriminate between the alternative hypotheses of formation of the mid-crustal magma bodies, its results are in agreement with the model of DeCelles,  in terms of the significance of the large amount of radiogenic heat produced by a thickened crust. Our results show that the increased thickness of the crust and consequent elevated production of radiogenic heat have a major influence on the thermal field in the Central Domain and explain the high temperatures compared to the Eastern and Western Domains. Moreover, our model predicts positive thermal anomalies in those domains where crustal melts are inferred from seismic velocity and electric conductivity (e.g., Bianchi et al., 2013;Chmielowski et al., 1999;Ward et al., 2013Ward et al., , 2014. Modeled temperatures in these zones are close to the geothermobarometric estimations, suggesting that significant advection of material from a MASH zone is not required to explain the anomalies. Transient feedback effects between the different proposed mechanisms can, however, only be evaluated with dynamic models, which need to be developed in future research endeavors. In the Western Domain there are only heat flow measurements available for comparison with our model. Although the general trend is reproduced and the lowest heat flow values are modeled in this domain, the heat flow modeled for the Coastal Cordillera (60 mW m −2 ) is not as low as the available measured data, which can be as low as 20 mW m −2 (Figure 7). Low temperatures and heat flow in this domain are consistent with the mafic nature of the crust and the progressive crustal thinning toward the trench, which increases the thickness of the mantle within the model space. Crustal mafic rocks and mantle rocks are characterized by low and extremely low radiogenic heat production, respectively (Table 1), and the radiogenic contribution to the heat budget is therefore greatly reduced compared to that in the Central Domain. The presence of basaltic rocks in the upper crust of the Coastal Cordillera that are fractured and permeable allows infiltration and gravity-driven groundwater flow from the continent to the ocean (Burns et al., 2016;Flóvenz & Saemundsson, 1993), which would reduce the temperatures at shallow depths, and may explain the differences between modeled and measured heat-flows.
The temperatures and heat flow in the Eastern Domain are intermediate to those in the Central and Western Domains. Although the crust is as thick as in the Western Domain, the rocks are more felsic and produce more radiogenic heat (Table 1), thus increasing the heat input into the system. Our model predictions agree with bottom-hole temperatures recorded from depths ranging between 1,300 and 3,800 m (Collo et al., 2018) to within 8.2°C on average, with modeled temperatures in most cases cooler than those measured in the wells (Figure 11). The best fit for the heat flow data is achieved in this domain, although there is a local opposing trend between modeled and measured heat flow, with modeled heat flow decreasing toward the east but measured heat flow increasing in the same direction ( Figure 7). As in the Western Domain, convective processes in the form of groundwater flow could explain small differences between modeling results and observations. The Subandean Ranges are an area of recharge in which infiltration of cold meteoric water occurs (Larroza & Fariña, 2005), which would reduce shallow subsurface temperatures. As it infiltrates, the meteoric water is heated and, due to the high elevation and active thrusting in the ranges, it flows driven by pressure and gravity toward the foreland, possibly increasing the temperature of the foreland aquifers (Husson & Moretti, 2002). Since our conductive model does not take into account these aspects of coupled heat and fluid transport, it may tend to overestimate temperatures in the shallow subsurface of the Subandean Ranges and underestimate temperatures in the foreland.
In spite of the limitations of our thermal modeling approach, we propose to take the next step and investigate how variations in the thermal field, as controlled by the well-constrained geological structure (Ibarra et al., 2019), affect the rheological state of the lithosphere.

Rheological Model and Active Deformation
The modeled lithospheric strength configuration strongly depends on the modeled temperature distribution whose validity we discussed above. The other controlling parameters are the background strain rate and the lithology-dependent mechanical rock properties. For our modeling we used a uniform strain rate for the entire region, although observations suggest a reduction in the present-day strain rate from north to south and strain localization along the edges of the orogen (Gerbault et al., 2003;Schemmann, 2007). Since we were interested in the long-term strength variations resulting from thermal and compositional heterogeneity, the application of a spatially invariant strain rate allowed us to isolate the effect of those heterogeneities on the lithospheric strength distribution.
The low integrated strength levels in the crust and lithosphere of the Central Domain (Figures 8a and 8b) and the vertical distribution of the yield strength beneath the plateau (Figure 9d) suggest a significantly weakened lithosphere within the orogen, particularly beneath the Altiplano-Puna plateau where the upper Figure 11. Modeled temperature plotted against bottom-hole temperature measurements (red circles). The dashed black line represents a 1:1 relationship between X and Y axes. The geographical location of the corresponding wells is shown in Figure 7. crust is the only strong layer. In the Eastern Domain, the Western Domain, and the Atacama Block the integrated strength is higher and the yield strength envelopes indicate a stronger lower crust and/or mantle (Figures 8a and 8b, and Figures 9a and 9b, 9c & 9e). These results suggest a clear first-order spatial correlation between thermal field and rheology: Regions with high temperatures coincide with regions of low integrated strength.
The effect of rock-type rheology on strength is less evident but is still discernable. Figures 6, 8 and 9 illustrate that although the crust in the Atacama Block has higher temperatures than the foreland (Figure 6), the crustal integrated strength is higher in the former than in the latter (Figure 8b). The yield strength envelopes also show a considerably stronger crust in the Atacama Block (Figure 9e) than in any of the morphotectonic units in the deformed foreland (Figures 9a and 9b & 9c). This lack of correlation between thermal field and integrated strength can be explained by the different rock-type rheology of the units and the different crustal thicknesses; the mafic rocks in the Atacama Block are stronger than the rocks in the foreland ( Figure S1 in the Supplementary Information), and have a high yield strength, even at high temperatures ( Figure 9e).
In view of the observed clear dependence of strength on temperature and the strong control of crustal thickness and radiogenic heat production on the thermal field, we suggest that the crust (its thickness and composition) exerts a first-order control on the thermal field within the studied area, and consequently also on the strength distribution. This is to some extent in agreement with the results obtained by McKenzie et al. (2005), who modeled the thermal structure of oceanic and continental lithospheres and suggested that Moho temperatures beneath the continents are primarily dependent on the thickness of the crust and its radiogenic heat production.
Our results also indicate that the strength of the crust and lithosphere may influence the location of deformation processes. Active deformation in the form of recorded seismicity and neotectonic structures is preferentially located in areas with high strength gradients, that is, in the transition zones between weak and strong domains (Figures 8d and 12). Previous investigations have also reported concentrations of seismicity and active deformation in transition zones and zones of weakness: Fernández-Ibáñez and Soto (2008) investigated crustal rheology and earthquake distribution in the Gibraltar Arc and found that most of the crustal seismicity in the region occurred in domains of medium to high strength; Ito (1990) determined the regional variations of the seismogenic layer and the brittle-ductile transition zone in a region of Japan and noted that large earthquakes appeared to nucleate where there were abrupt changes in the depth of the seismogenic layer and the BDT, associated with variations in the thermal field; Sloan et al. (2011)  that crustal seismicity, and in particular seismicity in the lower crust, was focused on the transition zone from a cold, thick lithosphere to a hot, thin lithosphere. Similar research with global and regional foci by Mooney et al. (2012) and Assumpçao et al. (2004) found that intraplate seismic activity concentrated around steep lateral gradients in lithospheric thickness, and in lithospheric thin spots. More recently, Magrin and Rossi (2020) found that most of the seismic events in the northern tip of the Adria microplate occur in the front of the Southern Alps facing south, where there is a sharp transition from high to low values of seismic velocity, density, rigidity and Young's modulus. South of our study region, Rodriguez Piceda et al. (2020) found that the highest horizontal surface strain rates in the southern Central Andes occur in areas with abrupt variations in crustal density, suggesting an influence of crustal composition and associated rheological properties on short-term deformation. Furthermore, Barrionuevo et al. (2021) have shown by means of geodynamic modeling how such variations in crustal composition influence the deformation mode the orogenic wedge.
These observations could be related to lateral variations of gravitational potential energy (GPE) as suggested by Ibarra et al. (2019). The GPE of a mountain range grows with the square of both its average elevation and the thickness of its crustal root (Stüwe, 2007); as a result, the forces required to produce further vertical growth increase exponentially, preventing strong compressional deformation (e.g., Molnar & Lyon-Caen, 1988). Regions with high GPE exert a net force on those with low GPE, imposing a compressional stress field on the transition zone (e.g., Stüwe, 2007). We computed the GPE from the density configuration of our model to calculate the horizontal forces arising from differences in GPE between the region and a mid-ocean ridge ( Figure S4a in the Supplementary Information). The net forces were then compared to the modeled integrated strength of the lithosphere to identify areas prone to deformation ( Figure S4b in the Supplementary Information). The resulting forces vary in the same order of magnitude as the modeled strength, suggesting that the static density configuration of the system is able to contribute to deformation. However, the areas that seem to be most prone to deformation do not coincide spatially with most of the regions with observed active brittle deformation, which strongly points to the effects of tectonically induced stresses. Although tectonic stresses are not accounted for in our calculations, local and global estimations of the horizontal deviatoric stress field arising from buoyancy and tectonic forces show indeed that the Central Andes are characterized by extensional deviatoric stresses within the orogen interior, and compressive deviatoric stresses along its margins (e.g., Flesch & Kreemer, 2010;Ghosh et al., 2009).
A different perspective, particularly in extensional and cratonic settings with smoother topography where GPE variations would not produce such a pronounced build-up of stresses, is that of Cacace et al. (2008), who used a finite element viscous thin sheet model to investigate the influence of contrasting rheological structures on deformation and stresses in the Central European Basin System. Their results showed strong strain localization along lateral heterogeneities and variations in the direction of the principal stress axes. In line with these results, the compressive dynamic simulation of the Central Andean 3D density structure (Ibarra et al., 2019) showed strain localization along the boundary between the felsic crust of the Central Domain and the more mafic crust of the Eastern Domain.
The observed localization of active deformation along high strength gradients could influence the position of the volcanic arc as well. Although the position of volcanic arcs in subduction systems is known to be controlled by the dehydration of the subducting slab at approximately 100 km depth (Tatsumi, 1986), the marked diversion of the volcanic arc (i.e., the Western Cordillera) toward the east, around the Atacama Block (Figure 1), is not associated with any change in the dip of the slab. Since the volcanic arc lies approximately along the boundary between the strong crust in the forearc and the weak crust in the orogen (Figure 8b), we suggest that the ascent of magma in this zone was facilitated by pathways formed by localized opening of fractures, thus causing the diversion of the volcanic arc.
The seismic activity and neotectonic structures in the Coastal Cordillera and the Puna do not show a clear correlation with the high strength gradients. Most of the seismicity in the Coastal Cordillera is confined near the interface between the subducting slab and the upper plate (e.g., Allmendinger & González, 2009;Bloch et al., 2014, Figure 12). The Puna plateau presents a complex and areally extensive deformation system that has contributed to the Late Cenozoic uplift history (e.g., DeCelles, Ducea, et al., 2015;Pingel et al., 2020), with N-S striking reverse faults and fold axes (e.g., Montero-López et al., 2021), WNW-ESE striking left-lateral transtensional fault systems, NE-SW striking right-lateral transpressive fault systems (e.g., Norini et al., 2013;Riller et al., 2001;Zhou & Schoenbohm, 2015), and superposed normal faulting (e.g., Allmendinger et al., 1989;Montero-López et al., 2014). Although the widespread extensional processes in the southern Puna have been inferred to be associated with delamination based on geophysical and geochemical studies (e.g., Beck et al., 2015;Calixto et al., 2013;Drew et al., 2009;Murray et al., 2015), Riller et al. (2012) have shown by means of analog modeling that a N-S transition from weak to strong crust may produce complex deformation features as found in the Puna. Our results support the existence of such N-S transitions within the Central Andes, from a weak crust beneath the Puna plateau to a strong crust beneath the Pampean Ranges (Figure 8).
Some of the early investigations into the distribution of seismicity with depth suggested that brittle failure was restricted to depths above the 600°C isotherm in the mantle and the 300°C−450°C isotherms in the crust because, as a first approximation, mantle rheology is governed by the rheology of olivine and crustal rheology is governed by the rheology of quartz or plagioclase (e.g., Chen & Molnar, 1983;McKenzie et al., 2005). However, within our study area deep crustal seismic activity has been observed well below the depths of the 300°C−450°C isotherms (Figures 10 and 12). Similar observations from different tectonic settings in western North America, the Dead Sea basin, the East African and Baikal rift systems, the Tien Shan foreland, and the northern Alpine foreland of Switzerland, have challenged the early ideas and called for a more thorough investigation on both the seismic process and the rheological stratification of the lithosphere (Aldersons et al., 2003;Bryant & Jones, 1992;Camelbeek & Iranga, 1996;Deichmann, 1992;Déverchère et al., 2001;Nyblade & Langston, 1995;Sloan et al., 2011;Wong & Chapman, 1990).
A number of authors have found a good correlation between the location of brittle-ductile transitions and various proxies for the depth-frequency distribution of hypocenters, based on the establishment of rheological profiles of the lithosphere (e.g., Fernández-Ibáñez & Soto, 2008;Meissner & Strehlau, 1982;Sibson, 1982). The locations of the brittle-ductile transitions in our model show a correlation (within the range of uncertainty) with the depth at which the frequency of hypocenters starts to decrease with depth (D dec ). This implies that a proportion of earthquake hypocenters lies within the ductile regime. Although there are uncertainties in the depth of the BDTs inherent to the modeling approach, independent observations such as radial anisotropy that suggest crustal flow as shallow as 10−20 km in the Altiplano-Eastern Cordillera region (Lynner et al., 2018) are consistent with our modeled BDTs. The occurrence of earthquakes below brittle-ductile transitions has been explained as being due to the propagation of aftershocks, which is not the case for the region given the time-space distribution of events (declustered catalog in Petersen et al., 2018), or the presence of strong lithologies that remain brittle at higher temperatures (e.g., Fernández-Ibáñez & Soto, 2008;Sloan et al., 2011). Likewise, the Quaternary mafic volcanism in the Puna, the Ordovician magmatism in the Pampean Ranges, and the mafic extrusions of the Cretaceous Salta Rift in the Santa Bárbara System (e.g., Drew et al., 2009;Rapela et al., 2018;Viramonte et al., 1999) may have left behind strong mafic residual rocks in the crust that could favor brittle behavior at greater depths and act as focal points (or asperities) for the initiation of seismic slip. Interestingly, most of the hypocenters in the region are located above the 600°C isotherm (as in Figure 12) that bounds areas inferred to correspond with sectors that exhibit brittle behavior of olivine.
The long-lived deformation along the flanks of the orogen and the development of associated shear zones could also explain the observed deep seismicity, particularly in the Santa Bárbara System and the Pampean Ranges where there has been repeated reactivation of deep crustal structures (e.g., Hongn et al., 2007;Iaffa et al., 2013;Kley & Monaldi, 2002;Ramos et al., 2002). In Figure 12, some of the hypocenters within the Eastern Domain appear to lie on a west-dipping plane that could in fact represent a shear zone formed during the activity of the Cretaceous Salta Rift (e.g., Grier et al., 1991;Kley et al., 2005), as suggested by the seismotectonic studies of Cahill et al. (1992). The presence of such deformation zones could have a weakening effect on the lithosphere, reducing strength and favoring strain localization. For example, a reduction in the friction coefficient from 0.75 to 0.2, as reported for natural fault systems (e.g., Zoback et al., 1987), would result in brittle-ductile transitions as deep as the Moho in the Pampean Ranges. Unfortunately, the existence of localized compositional heterogeneities and deformation zones is not resolved in our density model (nor consequently in the rheological model) because of its regional character and coarse resolution. Further models would help to improve our understanding of the feedback effect between lithospheric strength and fault growth, in which strength controls the location of faulting, and faulting reduces strength, thereby favoring strain localization.
With respect to the relationship between T s and the rheological stratification of the lithosphere in our model, we note that T s appears to be greater in those regions that have a jelly sandwich type of strength distribution (Pampean Ranges and Santa Bárbara System; Figures 9b and 9c and Figures 10b and 10c), compared to areas with crème brulée lithospheric profiles (Puna; Figures 9d and 10a). This observation is consistent with the variation in thermal gradients between the different regions: relatively low thermal gradients associated with jelly sandwich strength profiles result in increased ductile strengths, and therefore deep brittle-ductile transitions and high mantle strengths, compared to those in regions with high thermal gradients and crème brulée strength profiles.
We have so far only discussed the earthquake distribution within the crust. With regards to the mantle, the only seismic activity in our modeled region is located within the subducting Nazca plate. However, some earthquakes in the Pampean Ranges occur close to the Moho (Figures 9 and 10); if possible errors in the vertical positions of hypocenters and the Moho depth are taken into account, the mantle could possibly present some intraplate seismicity. Moreover, the lack of dense seismic networks in the region needs to be considered as a bias in data availability. The possible presence of earthquakes in the mantle could be associated to those long-lived faults in the Andes that are likely rooted in the mantle, where different weakening mechanisms and brittle seismogenic asperities would allow seismic slip to initiate and propagate (Inbal et al., 2016;Prieto et al., 2017). In the case where all earthquakes occur in the crust, the absence of intraplate seismicity in the upper mantle of the Central Andes can be explained either by the mantle being too weak (beneath the orogen) or by the level of stress being insufficient to produce failure in a strong mantle (beneath the foreland).

Conclusions
We have modeled the steady-state conductive thermal field for the Central Andes on the basis of our previously published data-derived 3D density model of the region. The thermal model was subsequently used to derive the strength distribution considering brittle and ductile mechanisms of deformation in the lithosphere. The results were compared with the distribution of earthquakes and neotectonic structures. We found that: (1) The deep temperature distribution is primarily controlled by the thickness and composition of the crust. The thick and felsic crust with high radiogenic heat production in the orogen presents higher temperatures than the forearc and foreland, despite these two areas having higher temperatures at the lower boundary condition (2) Although convective/advective processes are active in some parts of the Central Andes and thermal equilibrium has probably not yet been attained, the thermal model reproduces measured surface heat flow and shallow crustal temperature reasonably well in those regions where advective processes are less pronounced (3) The integrated strength of the lithosphere is mainly controlled by the temperature distribution with the hot orogenic lithosphere having a lower integrated strength than the cooler forearc and foreland areas (4) Active deformation in the form of seismic activity and neotectonic structures occurs preferentially in areas that present high horizontal gradients of integrated strength. These areas also exhibit lateral variations of crustal densities, seismic velocities, and lithospheric thickness, as has also been reported from a variety of different tectonic settings. We suggest that the presence of such transitions zones, together with the high gravitational potential energy within the orogen, induce strong strain localization (5) Some earthquakes occur beneath the modeled brittle-ductile transitions, suggesting that possible seismic asperities and weakening mechanisms associated with long-lived structures in the Andes could exert control on the hypocenter depth distribution and the strength of the lithosphere