Controls of the Lithospheric Thermal Field of an Ocean-Continent Subduction Zone: The Southern Central Andes

Temperature exerts a first-order control on the rheology of the lithosphere, affecting the depth of the brittle-ductile transition zone and the occurrence of thermally activated creep processes [1, 2]. In an active ocean-continent convergent plate-boundary system, the assessment of the temperature distribution within the lithosphere is crucial for understanding the mechanisms controlling subduction geometry (e.g., slab dip and subduction-zone curvature) and the localization of deformation within the upper plate, including the vergence of thrust belts [3–8]. Estimating the thermal state of the system is challenging, however, as it requires deciphering the complex and continuous interplay between different heat transport mechanisms (conduction and convection) and heat sources. These across-scales phenomena include heat conduction, advection of the oceanic plate that steadily supplies colder material, variations of thermal properties within the plates (radiogenic heat production and thermal conductivity), frictional heating along the subduction interface, latent heat due to mineralogical phase transformations within the oceanic plate, and adiabatic heating in the asthenosphere (for a review, see [9, 10]).

Many modeling studies [11–15] have considered the temperature distribution of the lithosphere as the main driver of the dynamics of a subduction system. These studies adopted a simplified configuration of the upper plate in terms of its thickness and rock composition [14]. However, the continental lithosphere is the product of a complex tectonic and magmatic history, involving terrane amalgamation, trench erosion, crustal thickening during subduction and collision, and extensional overprinting either during the final stages of orogeny or from deep-seated processes in the mantle [16, 17]. To address these considerations, other modeling studies have effectively shown how the present-day thermal field varies with respect to thickness and compositional heterogeneities within the lithosphere [18–23] and how these variations ultimately affect the long-term rheology of the lithospheric plate [5, 23–25].

The southern Central Andes (SCA, 27°–40°S, Figure 1) are a suitable region for studying the effect of both a heterogeneous configuration of the upper plate and differences in subduction geometry on the lithospheric thermal field in an ocean-continent convergent plate-boundary system for several reasons. Bordering the convergent margin between the oceanic Nazca Plate and the continental South American Plate, the SCA encompass several morphotectonic provinces, each characterized by a distinct pre-Cenozoic geological history and lithospheric configuration [26, 27]. These differences are rooted in complex tectonic and magmatic episodes of shortening and extension that span from the Neoproterozoic to the Quaternary [26, 28–35]. Major pulses of Andean deformation are thought to have occurred during the Late Cretaceous and Miocene [36, 37], when the style of deformation was significantly influenced not only by the characteristics of the subducting plate [38–40], but also by the reactivation of inherited tectonic heterogeneities, which influenced Cenozoic phases of erosion, sedimentation, and geomorphic evolution [41]. The SCA are subdivided into four first-order morphotectonic provinces: the forearc, the magmatic arc, the back-arc, and the foreland (Figure 1; [3, 26, 42]). These extensive regions are in turn subdivided into morphotectonic provinces that are characterized by distinct structural and geomorphic features. For the purposes of this study, we have defined as the Andean orogen only the magmatic arc and backarc provinces of Principal and Frontal Cordilleras, and thus not included the fold-and-thrust belt of the Precordillera and the Payenia volcanic province. The main features of the foreland are the reverse-fault-bounded basement uplifts of the Sierras Pampeanas and the Neuquén and Cuyo basins (Figure 1).

A distinct feature of the SCA is the variation in the subduction angle along the strike of the subduction zone, transitioning between 33°S and 34.5°S from subhorizontal (<5°) in the north (Chilean-Pampean flat-slab area; [43]) to relatively steep (~30°) in the south (Figure 1; e.g., [44, 45]). Even though the present-day subduction regime has been active since at least Early Jurassic times ([46] and references therein), the flattening of the slab north of 33°S presumably began at ~19 Ma [47–49], finally achieving its subhorizontal configuration at ~7–6 Ma (for a review see [31, 43, 50]). It has been proposed that this flattening event is responsible for the absence of present-day magmatism between 27° and 33°S (Figure 1) [31, 50]. Several causal mechanisms have been suggested for triggering flat subduction at these latitudes, including (i) buoyancy of the slab due to the subduction of an aseismic ridge [50–52]; (ii) fast trenchward motion of the overriding plate that inhibits slab rollback and drives the trench to retreat [53]; (iii) stagnation of the oceanic plate in the mantle transition zone [54, 55]; and (iv) enhanced coupling between the oceanic and continental plates due to the greater strength of the continental plate [56–59].

In the SCA, surface heat flow, shallow subsurface temperatures, and magmatic activity vary significantly across the subduction system, which generally has been attributed to the geometry of the subducting Nazca plate [53, 63–71]. Earlier studies [63, 64, 68, 72], based on lower surface heat flow values (20–70 mWm^-2 v_s. 50–120 mWm^-2), suggested that the flat-slab segment is colder than its steeper counterparts, a hypothesis which was also supported by geodynamic numerical modeling [53, 65, 66] and seismic tomography [67, 69–71]. According to the latter studies, low v_p/v_s ratios (P wave/S wave velocity <1.75) characterize the flat-slab segment in the northern part of the SCA, in contrast to the higher v_p/v_s ratios encountered to the south, which are within the range typically found for most subduction zones [73]. The thermal contrast between the two differently dipping segments in the SCA is commonly linked to variations in the extent of the mantle wedge and arc magmatism, both of which are drastically reduced in the flat-slab segment [31, 73–77]. All previous studies therefore seem to suggest that the forearc is an area with low surface heat flow in response to the subduction of cold oceanic crust at shallow depths.

These interpretations have recently been challenged by local-scale studies which indicate that parts of the thermal variations in the lithosphere of the SCA are not related to the effect of the subducting plate [25, 78–82]. For example, Sánchez et al. [80, 81] provided evidence for a significant difference in surface heat flow between the orogen (85–95 mWm^-2) and the foreland (~45–60 mWm^-2) at the latitudes of the flat-slab region and proposed structural and/or compositional variations within the crust and different heat flow input at the base of the lithosphere as possible reasons for this phenomenon. However, none of these hypotheses have yet been validated by a detailed study of the configuration of the lithosphere in terms of its geometry and composition.

In light of the open questions, we aim to test in this contribution how mantle thermal anomalies and first-order structural and lithological heterogeneities in the overriding plate across and along strike of the subduction system affect the thermal field of the SCA. In order to do so, we followed a data-based modeling approach. Seismic shear wave velocities [83] were converted to temperatures to obtain the deep thermal field in the mantle and across the subduction interface. In the model domains shallower than 50 km, where no mantle shear wave velocity data are available for conversion to temperatures, we calculated the steady-state conductive thermal field. To test the assumption of thermal equilibrium, we analyzed the effects of time-dependent processes related to subduction dynamics by carrying out a transient numerical analysis for the portion of the model domain with lack of constraints from the seismic tomography. Steady-state conductive temperatures in the shallow model domain were computed based on an existing 3D structural and density model of the SCA [58], which is consistent with available geological and geophysical data. We assigned thermal properties to the sedimentary cover rocks, the crystalline crust, continental lithospheric mantle, and oceanic plate according to the lithological characteristics of the units that previously have been derived from surface geology as well as from seismic velocities and modelled densities [58]. This allowed us to evaluate the control of the lithospheric structure on the resulting temperature distribution. The validity of the inferred thermal structure is assessed by comparison with temperature and surface heat flow measurements available for the studied area and a detailed sensitivity analysis of the model. One main challenge of this approach is related to the sparse coverage of thermal measurements in certain parts of the model. To address this problem, our results are qualitatively compared with other proxies of the thermal state of the area with larger spatial coverage, including seismic attenuation and elastic thickness patterns. As a result, we have obtained a 3D thermal model of the SCA and adjacent foreland regions that describes the relative temperature variations between the geological units of different composition and lateral and depth extents. Finally, this model allows us to make a qualitative analysis of the thermal feedback mechanisms between these different geological units.

The main thickness and density variations of the layers constituting the SCA lithosphere were recently described in a 3D lithospheric-scale, density, and structural model of the SCA [58]. This model was constrained by an array of geological and geophysical data, including seismic reflection and refraction profiles, seismic tomography, sediment-isopach maps, and gravimetric observations ([61, 83–86]; see references in [58]). This structural geological model covers a region of 700 km by 1100 km with a horizontal resolution of 25 km and a depth of 200 km below mean sea level (bmsl), comprising the forearc, the Andean orogen, and the foreland regions. The vertical resolution varies as a function of the thickness of the corresponding layers, which were mainly defined on the basis of density contrasts. These layers comprise from top to bottom: (1) water; (2) marine sediments; (3) continental sediments; (4) upper continental crystalline crust; (5) lower continental crystalline crust; (6) continental lithospheric mantle; (7) shallow oceanic crust; (8) deep oceanic crust; (9) oceanic lithospheric mantle; and (10) oceanic sublithospheric mantle. Figure 2 illustrates the main structural features of the 3D model (see Rodriguez Piceda et al. [58] for more details).

Overall, maximum sedimentary thicknesses occur in the Cuyo and Neuquén basins (Figure 2(a)). The Andean orogen has a thicker crystalline crust (55 km) than the forearc (~35 km) and the foreland (~30 km) (Figure 2(b)). The remaining parts of the backarc and the foreland can be subdivided into three crustal domains: (i) a thick northern domain (40–60 km); (ii) a thin southern domain (~20 km); and (iii) a central domain with intermediate crustal thickness (35–45 km). The areas with the greatest upper crustal thickness comprise the orogen (20–40 km) and the Payenia volcanic province (20 km). In contrast, in the Neuquén Basin the upper crustal thickness thins out to 5 km (Figure 2(c)). The greatest lower crustal thickness (~30–45 km) exists in the northern part of the backarc and foreland regions of the Precordillera and Sierras Pampeanas (Figure 2(d)).

The general workflow followed in this study is illustrated in Figure 3(a). Figure 3(b) shows the 3D model box, where different thermal modeling approaches and corresponding thermal boundary conditions were applied. To predict the present-day thermal configuration of the SCA and discuss its controlling factors, we subdivided the model volume into two domains: a deep domain between a depth of 50 and 200 km bmsl, where temperatures were converted from S wave seismic velocities (here referred to as “v_s-to-T conversion”), and a shallow domain, including the crust and uppermost mantle down to a depth of ~50 km bmsl, where the steady-state conductive thermal field was calculated using as input the 3D structural and density model of the area [58]. The reasons for this subdivision are as follows: the v_s-to-T conversion being developed for application to mantle rocks and the limited quality of mantle velocity data for depths shallower than 50 km [83].

To estimate mantle temperatures between 50 km and 200 km, we used the results of the S wave mantle tomography of Assumpção et al. [83]. This tomographic model is an updated version of the S wave tomography of Feng et al. [91] for the South American region, where the velocity structure of the upper mantle was constrained through the joint inversion of S and Rayleigh waveforms and fundamental mode group velocities of Rayleigh waves. The original data set is restricted to depths between 50 km and 400 km, with a horizontal resolution of approximately 25 km [91]. Our choice on the tomography of Assumpção et al. [83] rather than other global tomographic models covering the study area [92] stems from the fact that this model has a more refined lateral and vertical resolution and offers a better spatial correlation between high-velocity features and the track of the slab ([61]; Figure S1 in supporting information).

The implementation by Meeßen [88] calculates the v_s and the corresponding density at each depth in the seismic tomography for temperatures between 300 and 3000 K in steps of 1 K. For the density computation, lithostatic pressure is computed relying on the AK135 seismic model [93]. At each grid point, the algorithm compares the computed v_s with those from the tomographic model, by performing a look-up method over the table and choosing the two closest values to the velocity from the tomography. Then, the temperatures and corresponding densities are linearly interpolated to obtain the final values.

For this study, we chose different mantle compositions (spinel or garnet), listed in Table 1, according to the respective stable aluminum phase at depth [94]. For shallow depths (50–80 km), a mantle composition corresponding to a spinel lherzolite was assigned, in accordance with mantle xenoliths found in the Payenia volcanic province [95, 96]. Between 80 km and 200 km, the stable composition was assumed to correspond to garnet lherzolite [97]. Mineral properties α, ρ(P₀,T₀), M(P₀,T₀), ∂M/∂T, and ∂M/∂P were taken from Cammarano et al. [98] and Goes et al. [87]. The thermal expansion coefficient α was assumed constant for each mineral phase. The frequency exponent a and anelasticity parameters A, H, and V were taken from Sobolev et al. [99] (Table 2).

The temperature configuration derived from the v_s-to-T conversion in the parameter space of v_s and depth is shown in Figure 4, which also depicts three 1-D v_s profiles representative of the orogen at the latitudes of the flat slab, the steep slab, and the transition zone. In general, temperature increases with increasing depth and decreasing velocity. The largest temperature variations occur for depths shallower than 100 km and high v_s (>4.6 km s^-1), which is characteristic of the flat-slab domain.

Quantifying uncertainties in the v_s-to-T conversion is difficult due to the combined effects of uncertainties related to the conversion parameters (i.e., anelasticity model and mantle composition, [100]). An additional source of uncertainty is the S wave tomography, as a 0.1% perturbation in v_s, for example, can translate into temperature variations of 50°–250 °C. Previous studies [87, 100] assumed a temperature uncertainty of 150 °C in the v_s-to-T conversion at Moho depths, but recognized larger uncertainties at greater depths. We therefore reexamined the uncertainty of the conversion method of Goes et al. [87] by testing the model sensitivity with respect to mantle composition X, thermal expansion coefficient α, and attenuation Q. Regarding the conversion method, although several approaches exist [101–104], we only tested the method of Priestley and Mckenzie [102] as implemented by Meeßen [105] (Model PM). In a first step, we set up a model (“reference model” hereafter) based on the following parametrization: (i) mantle composition X corresponding to a garnet lherzolite (Table 1); (ii) constant expansion coefficient α; and (iii) anelasticity Q by Sobolev et al. [99] (Table 2). In a second stage, we tested alternative models by varying one parameter at a time: (i) mantle composition (spinel model), (ii) thermal expansion coefficient (α model), and (iii) attenuation (Q2 Model). In all cases, the conversion was limited to the depth interval of 50-200 km, as thought to be representative of the lithospheric mantle. A detailed description of the setup of these alternative models is provided in the Supporting information (Text S1). In Section 3.1, we discuss in greater details the results from the same sensitivity analysis as applied to the tomography of Assumpção et al. [83], while in Section 4.1.2, we open a discussion on the implications of this sensitivity analysis on the modeled deep thermal field.

Temperatures were calculated with the finite element code GOLEM [89, 90]. For the thermal computation, three modifications were made to the original 3D configuration of Rodriguez Piceda et al. [58]. First, the water layer was removed, thus treating the topography/bathymetry as the top of the model (cf. Figure 1). Second, the horizontal resolution was increased from 25 km in the original structural model to 5 km, and, third, the layers were vertically refined by a factor of 3 to 32 in order to ensure that (i) each layer has at least three finite elements and (ii) most of the model domain is represented by a cubic finite element to ensure faster numerical convergence (Figure 3(c)). These modifications ensured to properly solve the temperatures in each node of the mesh without significantly increasing computational time.

Each unit of the 3D lithospheric model was populated with constant thermal properties (Figure 3(a); bulk conductivity λ and radiogenic heat production S) according to its main lithology (Text S2 in supporting information). The characteristic lithologies, in turn, were selected based on the comparison between gravity-constrained densities [58] and mean P wave velocities [67, 106–109], combined with rock-property compilations [110, 111] and other seismic properties [70, 112–117]. A range of thermal properties [118–122] related to the chosen lithology for each layer was additionally tested until the best fit was achieved with a compilation of borehole temperatures mainly limited to the foreland basins [79]. Table 3 summarizes the chosen values for each layer of the final (best fitting) model. A sensitivity analysis of the model results to the tested range of the thermal properties indicate that the modeled temperatures are most sensitive to variations in the thermal conductivity of the upper continental crystalline crust and the mantle (Text S3 in supporting information).

To close Equation (6), Dirichlet boundary conditions (i.e., fixed temperatures) were assigned along the top and base of the model. The upper thermal boundary condition was set at the topography/bathymetry (Figure 3(b)), with temperatures extracted from the ERA-5 land database ([123]; Figure 3(d), Text S4 in supporting information). The lower boundary condition was set at the depth of the upper bound of the v_s-to-T conversion: a constant depth of 50 km bmsl for areas where the Moho is shallower than 50 km bmsl and at a surface 5 km deeper than the Moho where this interface is deeper than the abovementioned threshold (Figure 3(b)). The reason for choosing this depth for the boundary condition stems from the fact that at shallower depths, the seismic velocity from the tomography of Assumpção et al. [83] reflects both mantle and crustal components. The temperature distribution at this boundary was derived from the v_s-to-T conversion (Figure 3(e); [83]; Section 2.2).

From the conversion of S wave velocities taken from the tomography of Assumpção et al. [83], we obtained the lower boundary condition of the steady-state conductive model (Figure 5(a)) and the mantle-temperature distribution for the ~50–200 km depth interval (Figures 5(b)–5(d)). Across the lower boundary condition, temperatures roughly range between 600 °C and 1000 °C. Two domains with temperatures of <700 °C are identified: (i) a cold nose (CN) between 70°W and 72°W, beneath the forearc in the central and northern portion of the study area, and (ii) a domain farther to the east where the slab flattens (FS) between 29.5°S and 32.5°S. The CN extends eastward above the transition zone and the steep-slab segment and is significantly attenuated above the flat-slab segment (Figure 5(a)). With increasing depth (>75 km bmsl), this thermal feature is no longer visible, while the cold FS domain extends vertically over the entire mantle column of the overriding plate (Figures 5(a)–5(c)). At 80 km bmsl, temperature increases to the SW, with maximum values (~1100 °C) located between 37° and 39°S and 70° and 72°W (Figure 5(b)). At 125 km bmsl, the temperatures follow a similar pattern as at 80 km depth, but differ in absolute value, with the cold FS domain reaching temperatures between 850 °C and 900 °C (Figure 5(c)). Towards the marine domain and to the south of the study area, temperatures increase up to 1300 °C. At the base of the model (200 km bmsl), temperature ranges between ~1200 °C and ~1350 °C, with the lowest temperatures correlating spatially with the track of the slab ([61]; Figure 5(d)).

We evaluated the effect of uncertainties in the parametrization of the v_s-to-T conversion applied to the seismic tomography of Assumpção et al. [83] by comparing the model described above with alternative model scenarios (for details, see Text S1 in Supporting Information). Figure 6 shows these alternative models in terms of the residual temperature at the lower thermal boundary condition of the steady-state model and at depths of 80, 125, and 200 km. The S wave velocity distributions [83] at those depths are also depicted.

The change from a garnet-and-spinel model (reference model; Table 1) to a mantle composed only of the spinel variant exerts the strongest effect on the resulting temperature at shallow depths and particularly where v_s is >4.5–4.6 km s^-1. The greatest temperature difference between the reference and the spinel models is 200 °C (Figures 6(e) and 6(f)). In our study region, v_s of such high magnitudes characterizes the forearc and the flat-slab segment in the northern part of the foreland (Figure 6(a)). The temperature residuals decrease with increasing depth, where lateral variations of v_s are less distinct (Figures 6(g) and 6(h)).

The conversion method of Priestley and Mckenzie [102] yields the largest temperature differences among all the alternative models tested, predicting up to ±400 °C difference for areas with high (>4.5 km s^-1) and low (<4.5 km s^-1) v_s, respectively (Figures 6(q) and 6(r)), as compared to the reference model converted based on Goes et al. [87].

From the steady-state conductive approach, we computed the thermal field of the crust and uppermost mantle to the depth of the lower boundary condition (Figures 7 and 8) and the surface heat flow (Text S5 in Supporting Information; Figures 8 and 9). Figure 7 shows the temperature distribution at depths of 2, 5, 10, 20, 25, 30, 35, 40, and 45 km bmsl and at the Moho. In Figure 8, we plot the modeled thermal field and the configuration of the main lithospheric units together with variations in predicted surface heat flow compared to available observations [63, 68, 126], along three representative E-W cross sections of the flat slab, the transition zone, and the steep-slab subduction segments.

At 2 km bmsl, the temperatures range between 15 °C and 165 °C (Figure 7(a)). As expected, the warmest areas are those with the highest topography (4–6 km height; cf. Figure 1) and the largest upper crustal thickness (30–40 km, cf. Figure 2(c)), which correspond to the central and northern segments of the Andean orogen. In addition, the forearc is characterized by an overall lower temperature than the orogen (40–80 °C), but by a more pronounced lateral gradient with values increasing toward the Andes. The foreland and low-elevation backarc regions are characterized by a wide temperature range (60–115 °C), with warmer temperatures in the Precordillera, the Payenia volcanic province, and the Cuyo and Neuquén basins. Down to 25 km bmsl, the spatial trends of temperature distribution are similar to those observed at shallow depths, but with different absolute values (Figures 7(b)–7(e)).

The thermal contrast between the warm orogen and the relatively cold forearc, backarc, and foreland regions is more pronounced with increasing depth (e.g., differences of ~110 °C at 10 km bmsl and ~200 °C at 20 km bmsl; Figures 7(c) to 7(d) and 8). From 20 km bmsl downward, the temperature distribution partially resembles that of the lower boundary condition (Figure 5(a)). Lowest temperatures at these depths correlate spatially with the areas where the mantle is the coldest (CN and FS areas, Figure 5(a)). Here, the temperature minimum also correlates with a thick continental crystalline crust (~40 km, cf. Figure 2(b)) and thickened lower crust (>30 km, cf. Figure 2(d)).

The modelled surface heat flow varies laterally from minima of ~45–70 mWm^-2 in the oceanic domain, most of the forearc, and foreland to maxima across the orogen (80–100 mWm^-2; Figures 8 and 9). We observe a remarkable spatial correlation between the surface heat flow distribution and continental crustal features: whereas high heat flow corresponds to areas with thick upper continental crystalline crust (>25 km; e.g., within the orogen and the Payenia volcanic province; cf. Figure 2(c)), low heat flow characterizes the deep sedimentary basins (~ >3 km; e.g., the Neuquén Basin; cf. Figure 2(a)) and/or thick lower continental crystalline crust (>25 km; e.g., most of the Sierras Pampeanas; cf. Figure 2(d)).

As a first step in the validation process, we compared modeled temperature values with the published borehole measurements of Collo et al. [79], located mainly in the central and northern foreland (Figure 10(a)). The residual temperature (i.e., the difference between modeled and measured values) is shown in Figure 10(b). Figure 10(c) illustrates measured and modeled temperatures against depth, and Figure 10(d) shows the residual temperature against depth. In general, we obtain a good fit of approximately ±20 °C between the borehole data and the modeled temperatures, with the exception of few outliers (Figure 10(c)). In the Cuyo Basin (33°S–69°W), positive temperature residuals tend to occur at shallow depths (<2 km), while negative residuals are detected at greater depths (Figure 10(b)). One plausible explanation for this trend is that the thermal conductivity of the Cuyo Basin sediments is slightly underestimated in our model. A larger thermal conductivity would translate in colder temperatures close to the basement surface and warmer temperatures close to the surface topography (see also sensitivity analysis in Text S3 of the supporting information).

Albeit limited in coverage, in a second step, we used available compilations of surface heat flow values within the SCA [63, 68, 126] to validate our thermal model. These measurements are not representative for the entire model domain as they are located mainly along the orogenic axis. They show a large variation in the magnitudes of the measured heat flow values (up to 250 mWm^-2), even between spatially close measurements (Figure 10(c)). Figure 10(d) depicts the residual surface heat flow, i.e., the difference between the predicted and the measured surface heat flow, at the location of the measurements. Figure 10(e) is a histogram of residuals of surface heat flow. In general, the model underestimates the surface heat flow along the orogenic chain with respect to the measured values, with only ~25% of the predictions matching the observations (Figure 10(e)). The model does not reproduce the extremely high heat flow values (>150 mWm^-2) reported for some volcanic areas in the axial sectors of the orogen. Additionally, due to its resolution, the model is not able to reproduce the observed variations in heat flow magnitudes between spatially adjacent measurements, which likely correlate to local features not resolved in our regional study (Figure 10(d)).

One main assumption in the calculation of the shallow temperature field was to consider that the lithosphere is in steady state. However, thermal equilibrium in the overriding plate can be disturbed by the advection of the cold subducting plate (e.g., [128, 129]). For instance, 2D kinematic models that consider this effect for the SCA [130] predict up to 200 °C colder temperatures at the subduction interface with respect to our model results (Figure S7 in Supporting Information). In light of these differences, a more appropriate modeling strategy would be to additionally account for this process. The caveat here is that performing such an analysis requires a detailed knowledge of the past temperature distribution in order to properly initialize the system. Unfortunately, we lack such constraints in the SCA. Given these considerations, we relied in our study on the assumption of steady-state conduction, where the results are less affected by the choice of the initial temperature condition, but are mainly determined by the model parameterization in terms of distribution of thermal properties and imposed boundary conditions (based on available observables in our study). Nonetheless, in an attempt to quantify the validity of this approximation for the SCA, we also computed a simulation that accounts for the additional effects of advection of cold temperatures due to the motion of the subducting slab on the present-day thermal field of the shallow domain. To that end, we first computed the resulting thermal field from the advection of a cold thermal front along the subduction interface (i.e., top of the oceanic crust). In a following step, we imposed this thermal evolution as the lower boundary condition on the 3D configuration of the overriding plate and ran a transient simulation with a duration of 7 Ma, which represents the past period during which the subduction geometry remained unchanged [26, 43]. A detailed description of this analysis is provided in the supporting information (Text S6). The comparison between the initial and final time steps at representative depth slices (10 and 40 km bmsl) indicates that the largest temperature difference (up to 370 °C) is registered in a narrow band within the forearc close to the subduction interface (Figure 11). Such a difference is due to the advection of the cold thermal front along this interface. In the remaining areas, temperatures at 7 Ma are up to 10 °C higher than the initial time step due to diffusion within the thick radiogenic crust in the orogen. In view of these results and the significantly larger lateral variations in temperature at different depth levels (Figure 7), we consider the assumption of thermal equilibrium as an adequate approximation for the thermal calculations of the shallow domain of the overriding plate in the SCA.

The model results depend on the parametrization of physical properties and boundary conditions. One source of uncertainty is related to the v_s-to-T conversion. While variations in the thermal expansion coefficient α and attenuation Q have negligible effects on the inverted thermal field, the sensitivity analysis for the v_s-to-conversion shows that using the conversion method of Priestley and Mckenzie [102] yields large temperature differences with respect to the reference model (Figures 6(q)–6(t)). One has to keep in mind though that the empirical conversion method of Priestley and McKenzie [102] is associated with large uncertainties (250–360 °C) for temperatures <900 °C [100, 131], which makes the results at shallow mantle depths highly questionable and the temperature distribution at 50 km bmsl (with reference model temperatures of mostly 400-700 °C; Figure 5(a)) a much less appropriate lower boundary condition for the performed forward thermal modelling performed.

To assess the effects of compositional variations on the derived temperature variations, we have also used the conversion method of Goes et al. [87] for two scenarios of upper mantle lherzolite variants, the reference model with a deeper garnet lherzolite modification (Table) and a spinel-only-composed alternative. The temperature difference between the reference and alternative models is only significant in a limited portion of the shallow mantle (<100 km) characterized by high v_s (>4.6 km s^-1), as is the case for the flat-slab segment and the forearc (Figure 6). Moreover, the regional thermal pattern in the mantle (i.e., the thermal contrast between the flat- and steep-slab segments) is a robust feature common to all model configurations despite variations in mantle composition. Temperature maps showing the difference between the two models when used as lower thermal boundary conditions for the forward thermal model are provided in the supporting information (Figure S12). Although the model with the alternative lower boundary condition predicts temperatures that are up to 80 °C higher than the reference model described in Section 3.2, the regional thermal heterogeneity remains unchanged, with contrasts between the forearc, orogen, and foreland, and between the flat- and steep-slab segments. From these observations, we can conclude that the modeled trends in temperature variations are within the same order of magnitudes, though still within its range of uncertainty, even when considering an alternative parametrization other than the preferred v_s-to-T conversion model.

Another limitation of the v_s-to-T conversion is related to the thermal structure within the slab. Although the oceanic plate displays higher v_s and lower temperatures with respect to the surrounding mantle (e.g., Figure 5; Figure S1 in Supporting information), the modeled thermal gradient within the plate is not as large as the one predicted by other analytical or numerical approximations of the thermal structures of subduction zones [14, 141]. The strong lateral contrasts of v_s are smoothed due to the resolution of the seismic tomography. This results in lower v_s and therefore higher temperatures than the predictions of these theoretical thermal models. On the contrary, v_s and resulting temperature anomalies within the continental mantle are of larger wavelength than those of the slab; thus, they can be captured by the longer wavelength surface waves of the seismic tomography. Therefore, our discussion mostly focuses on the thermal heterogeneities of the overriding plate and subduction interface.

A third limitation is that we do not consider additional effects from partial melting and fluids in the vs-to-T conversion. Melts in the mantle can lower seismic velocities due to an increase in anelastic relaxation [132, 133]. This effect depends mainly on the interconnection between melt pockets [134]. However, a proper quantification of this effect is hindered by inaccuracies in the estimation of the effect of anelasticity in high-temperature experiments [87, 135]. Further complexities are related to the presence of water, which would also lower the solidus and therefore would lead to an increase in the melting point [136]. Previous studies have suggested that an increase of melt fraction by 0.1 would cause a Vs decrease between 0.7% and 8.5% [87, 137], owing this range of variability to variations in melt geometry. However, the effect of these velocity variations due to partial melting on the converted temperatures below 50 km appears to be small (up to 20 °C, [138]). In the SCA, although we would expect partial melting to be somehow relevant in the subarc and forearc mantle due to the release of fluids from the oceanic plate [139], the seismic and/or geochemical information on the percentage of melt and fluids in the upper mantle (e.g., [67, 140]) is too limited to attempt any quantification of their effects on converted temperatures. We should nevertheless acknowledge that, in these areas, temperatures from seismic velocities should be interpreted with caution.

Due to its setup, the thermal modelling applied to the shallow domain (0–50 km) does not account for additional, nonconductive heat transport processes acting close to or within the subduction interface, including mantle-wedge convection and frictional shear heating [14, 130, 141–145]. For the case of viscous flow in the mantle wedge, 2D thermokinematic models of Central Chile calibrated with heat flow observations [127] indicate that this process only starts to dominate at depths larger than 80 km, which is below our forward modelled domain. The subduction interface at shallower depths (<50 km), in contrast, is largely influenced by shear heating effects [14, 127, 141]. In an attempt to quantify this process, we used the analytical approximation of England [141] for the thermal structure of the subduction interface (Text S6 in Supporting Information). We found that considering the addition of shear heating leads to up to ~40 °C higher temperatures along the subduction interface compared to a model without this heat source (Figure S11 in Supporting Information). This effect of shear heating is small compared to the temperature decrease caused by the advection of the cold oceanic plate (Section 4.1.1); thus, we consider that it would not greatly affect the temperature trends observed in the lithosphere.

Finally, additional methodological uncertainties relate to the limited resolution, coverage, and lateral differentiation of the lithospheric units in the 3D structural model, as well as to imposed thermal properties. Although there is an inherent non-uniqueness in the way thermal properties influence the results, the range over which these properties can vary is limited (see Text S3 in the supporting information). In addition to testing the effect of end-member property values, the use of borehole temperatures in combination with a wide variety of independent lithology-constraining data sets, including seismic tomography, seismic reflection and refraction data, and gravity anomalies, helped reduce the range of property variability. Future improvement in the definition of higher-order temperature contrasts would require more densely spaced seismic experiments focused on the deep crustal structure of the SCA and more wells including also extensive temperature measurements to cross-check the modeling results. Overall, the identified first-order thermal effects proved to be robust even for tested variations in imposed properties.

Our results indicate that the shallow thermal field of the lithosphere (<50 km) is largely controlled by the configuration of the continental crust, with temperatures varying according to the thickness of the sedimentary rocks and crystalline crust (Figures 7 and 9). Close to the surface (< 5 km), thick sedimentary basins (main depocenters of the Cuyo and Neuquén basins) exhibit temperatures up to ~40 °C higher than the basin margins (Figure 7). This is the effect of thermal blanketing produced by the low-conductive sedimentary layers [20, 21, 147, 148]. In contrast, the presence of thermally more conductive crystalline rocks leads to a more efficient heat transport where sedimentary cover rocks are absent and to colder shallow temperature at the same depth.

The temperature differences at 2 km and 20 km are shown in Figure 12. We observe that at shallow depths of 2 km, the temperature differences are indeed affected by variations in the topography (~100 °C in areas of 3–6 km elevation). This effect decreases with greater depth, although it is still evident at 20 km bmsl. At this depth level, temperature differences are up to 50 °C between areas of different topographic elevation. Below 20 km bmsl, the influence of the upper crustal thickness outweighs the topographic effect.

In contrast to the shallow domains of our model, for which the parametrization of the forward thermal modelling gives direct insights into the causes of obtained thermal anomalies, we rely on a more qualitative interpretation for depths of >50 km where the thermal field is derived from seismological observations. To identify and evaluate anomalous thermal trends in the mantle of the overriding plate, we have therefore subtracted a geotherm regarded as typical for continental lithosphere [146] from the tomography-derived temperatures in this deeper model domain. This normalization yields the volumetric extents of two major low-temperature domains: One that is known the “cold nose” (CN) along the trench and another one in the flat slap region (FS; Figure 13).

Below 50 km bmsl, we find that the lithospheric thermal field is dominantly influenced by both the cooling effect of the subducting slab and heat input from the mantle (thermally-driven mantle upwelling). The degree to which the slab dynamics affect the temperature distribution by advective cooling varies with distance from the trench and the subduction angle. The mantle of the overriding plate exhibits the lowest temperatures close to the trench, within the cold nose of the forearc, and where the slab flattens (mostly at 85 km bmsl beneath the Sierras Pampeanas). In contrast, temperatures at depths >50 km bmsl beneath the Sierras Pampeanas increase towards areas where the slab dips steeply (~30°). These results are consistent with the lower surface heat flow values observed [63, 64, 68]. Also these results are consistent with results from 2D thermomechanical modeling efforts for the area [67]. They suggest temperatures of 600 °C at the top of the flat slab (100–120 km bmsl)—a value that is consistent with the lower temperatures modeled across the subduction interface beneath the Sierras Pampeanas (Figure 13).

The spatial correlation between the temperature distribution of the overriding plate and the subduction angle breaks within the northern orogen above the flat slab between 20 and 100 km bmsl (29°–30°S and 70°–70.5°W; Figures 5, 7, and 13). In this domain, a temperature excess of up to ~250 °C is predicted along the remaining flat-slab segment (FS warm domain; Figure 13), with values similar to those of the steep-slab segment. As the subduction angle does not change, an additional forcing factor other than the subduction angle should be considered, to explain these locally elevated temperatures, a suggestion made previously, albeit not extensively discussed [80, 81]. Between 30°S and 33°S, our results indicate a spatial correlation between the thermal heterogeneity derived from the v_s-to-T conversion [83] and the configuration of the upper crust. The lowest temperatures in the subhorizontal slab segment occur where the upper crust is thin, while the highest temperatures occur where it thickens (cf. Figure 2(c)). This indicates that, in the warmer, northern part of the orogen, heating from a thicker upper continental crystalline crust outweighs the cooling effect of the underlying flat oceanic slab. Conversely, areas in the flat segment with a thin upper crust and a thick lower crystalline crust are significantly colder due to the smaller contribution of the lower crustal unit to the internal heat budget.

While the thick upper crust provides an explanation for low v_s (high T) below the orogen north of 30°S, the foreland of the flat-slab domain is underlain by a thinner upper crust and thus relatively high mantle temperatures (due to low v_s) there are more difficult to understand. It is likely that the resolution of the S wave tomography could have influenced the results within the latter area. North of the flat slab, high-resolution tomography [149, 150] identifies a N-S increase of v_s at 27°S with a sharp transition between low v_s below the southern Puna Plateau (that have been interpreted as resulting from delamination of the lower crust and the mantle; [151, 152]) and high v_s below the Sierras Pampeanas. Therefore, we speculate that this sharp velocity transition appears smoothed in the tomography of Assumpção et al. [83] due to its coarser resolution, resulting in lower modeled v_s, and therefore higher temperatures, in the foreland north of 30°S.

To investigate the degree of influence of mantle-related temperature variations on the shallow thermal field, we compared our results to a steady-state conductive model with a simplified lower boundary condition derived from a constant geothermal gradient of 5 °C km^-1 ([130]; Figure S13 in the supporting information), thus neglecting lateral thermal heterogeneities. At depth, there is no pronounced thermal contrast between the foreland of the flat and steep subduction segments, indicating that below 50 km bmsl, the main causative factors of the thermal heterogeneity of the foreland are the mantle heat input and the cooling effect of the slab. Moreover, the fit of the model with a simplified boundary condition and observed temperatures is decreased in which the modeled temperatures are 5 °C to 80 °C colder than the borehole data (Figure S14 in the supporting information). This implies that predicting the measured temperatures requires laterally variable heat input from the slab and the lithospheric mantle in addition to the heterogeneous crustal structure.

There is an ongoing debate concerning the importance of radiogenic heat production for the thermal field and the long-term evolution of orogens where the radiogenic crust is thickened and where thermal effects of shortening, exhumation, and partial melting are observed (e.g., [1, 153–156]). Some authors suggest that radiogenic heating has less influence on the thermal field of the lithosphere in the overriding plate than either shear heating along the subduction interface [157] or episodes of magmatic underplating [158]. Yet other authors have focused on the general thermal evolution of subduction orogens and the metamorphic record, proposing that shallow asthenospheric convection is the main process responsible for elevated temperatures in the lithosphere while disregarding any significant contribution from radiogenic heat production [159]. The latter study suggests that thermal equilibrium is achieved only 50 Ma after the main shortening phase, which contradicts the observation that peak metamorphism is synchronous with thickening [160, 161]. However, the interpretations concerning thermal equilibrium conditions in orogens have recently been disputed based on results from geodynamic numerical modeling [153]. These results suggest that radioactive heating during crustal thickening is responsible for the observed marked temperature increase and subsequent partial melting within the mid-crust after 30 Ma of shortening. The present-day SCA are within this time window after the main phase of shortening, which implies that crustal thickening that has taken place over more than 30 Ma could indeed explain a significant part of excessive surface heat flow [63, 64, 72] and the low seismic velocities observed across the orogen [67, 69–71, 162].

Predicted surface heat flow varies between the high surface heat flow domains over both the orogen and the Payenia volcanic province (80–100 mWm^-2), and low heat flow domains over the forearc, the remaining backarc, and foreland regions (50–60 mWm^-2 and 40–70 mWm^-2, respectively; Figures 8 and 9). These surface heat flow variations correlate spatially with the thickness configuration of the sedimentary strata and upper crystalline crust: high heat flow in areas with thin or absent low-conductive sedimentary rocks and thick radiogenic upper crust (i.e., Andean orogen, Payenia volcanic province) and low heat flow in regions with thick sedimentary cover and/or thin upper crust (foreland basins, forearc). Mantle-source contributions are found to play a less dominant role on surface heat flow than the abovementioned shallow structures. In the foreland regions, mantle heat flow reaching the base of the crust (Moho) is ~40% of the surface heat flow, whereas this contribution decreases down to 20-30% in the orogen (Figure S15 in supporting information). In contrast, the forearc is characterized by 60% of contribution from mantle heat flow to the surface heat flow, which is likely related to the absence of a thick upper crust and to the shallow underlying cold slab.

The observed trends between low and high surface heat flow in the forearc orogen is in agreement with results from 2D generic models of subduction zones [127, 163, 164]. The model of Wada and Wang [127] across a profile lying in the steep-slab region, however, predicts an overall lower heat flow than our model (50-75 mWm^-2 vs. 75-90 mWm^-2 Figure 8(c)). Colder conditions are likely due to differences in the model assumptions. Their model considers a lower radiogenic heat production of the upper continental crust than ours (1.3 μWm^-3 vs. 2 μWm^-3), which translates in higher temperatures, thus higher surface heat flow in our model. Additionally, the 2D model of Wada and Wang [127] does not consider lateral heat input and transport perpendicular to the 2D profile which could lead to an overall lower heat budget.

The effects of the spatially variable heat input (either from the deep mantle or from radiogenic sources) and variably efficient heat conduction demonstrate that the surface heat flow alone is a poor proxy for lower crustal or lithosphere thickness [20]. Nevertheless, heat flow estimates are often compared to seismic tomography, where S wave attenuation (Qs) anomalies correlate to some extent with the thermal field [165]. While areas of high heat flow are usually related to low Qs (high attenuation of the v_s), low heat flow is commonly associated with colder areas and hence high Qs (low attenuation of v_s). These spatial correlations were also identified across the SCA. Between 31.5°S and 33.5°S, high surface heat flow predicted for the orogen coincides with low Qs (650–670; [166]). Accordingly, low predicted heat flow correlates with high Qs (800–1050) in the forearc and the foreland [166].

Discrepancies between observed surface heat flow [63, 68, 126] and the steady-state conductive model predictions in the proximity of active volcanic centers (Figures 8 and 10(d)) can certainly be related to unconsidered transient advective heat transport processes, such as fluid migration or the existence of melts (e.g., [20, 167] and references therein). For instance, unconsidered effects of magma flow beneath the active volcanic arc can result in a local temperature increase of up to 200 °C [167]. González-Vidal et al. (2018) demonstrated that such local effects of partial melting can be interpreted from negative S wave anomalies imaged below the southern volcanic arc in the SCA area.

In view of the sensitivity of rock rheology to temperature, the lithological configuration from Rodriguez Piceda et al. [58] and temperature variations derived in this study can be analyzed qualitatively in terms of their implications for the long-term strength of the lithosphere. In particular, areas that are colder and of more mafic lithology (i.e., the northern part of the forearc and foreland) are potentially stronger and can withstand higher levels of horizontal stresses before deforming viscously [168] compared to areas that are warm and more felsic in composition (i.e., the orogen, the Payenia volcanic province, and the foreland at the latitudes of the transition to the steeper subduction segment).

The general trends in the modeled temperature distribution of the SCA lithosphere are consistent with independent elastic-thickness estimates derived from flexure analysis of the gravity field, which are an alternative, indirect proxy of lithospheric strength [25, 78, 80, 169, 170]. High and low elastic thickness are indicative of a strong and weak lithosphere, respectively [171, 172]. In the SCA, areas of high elastic thickness (40–60 km) correlate spatially with the modeled cold forearc and the northern part of the foreland. Conversely, areas of low elastic thickness (< 30 km) correlate spatially with the modeled warm areas of the orogen, the Payenia volcanic province, and the foreland at the latitude of the transition to the steeper subduction segment.

The inferred trends in lithospheric strength related to variations in the upper-plate configuration have strong implications for the long-term deformation processes of the Central Andes [4, 8, 173, 174]. For example, Barrionuevo et al. [4] argued on the basis of geodynamic numerical modeling that the vergence of the orogenic wedge at 33°–36°S is mainly controlled by the E-W-oriented asymmetry of the lithosphere-asthenosphere boundary (LAB). The LAB configuration proposed by these authors agrees with our study, where a warmer lithospheric mantle, and thus shallower thermal LAB, is encountered beneath the orogen compared to the adjacent foreland and forearc, i.e., areas where the mantle is colder and the LAB is located at greater depth. Furthermore, Barrionuevo et al. [4] suggest that heterogeneities in continental crustal composition could explain the observed N-S-oriented variations in the amounts of shortening and its spatial distribution within the orogen and the foreland between 33° and 36°S. In the orogen north of their study region, the upper and lower crustal deformation maxima are aligned vertically with the strongest overall crustal thickening (pure-shear or coupled deformation mode), whereas to the south, upper crustal deformation is horizontally displaced with respect to the locus of lower crustal deformation (simple-shear or decoupled deformation mode). Pure-shear deformation at 33°S would be mainly related to a more felsic and weaker crust; in contrast, simple-shear deformation at 36°S would result from a mafic and stronger crust [4]. These results are compatible with gravity-constrained density distributions in the crust in the same region as proposed by Rodriguez Piceda et al. [58] and with our results of an N-to-S-oriented decrease in crustal temperatures.

Furthermore, our results provide new aspects for the controversial debate over the governing mechanisms responsible for the spatially disparate thick-skinned deformation in the broken-foreland provinces of the Sierras Pampeanas between 27° and 33°S and the Santa Bárbara System farther north. While some authors have proposed that the flat slab is responsible for the thick-skinned deformation [26, 175, 176], others have argued that this style of deformation is controlled by the reactivation of crustal heterogeneities prior to slab flattening such as Paleozoic sutures and associated deformation fabrics between crustal terranes or the inversion of Cretaceous normal faults [6, 38, 41, 43, 174, 177–181]. In the case, that the flat slab is responsible, the concept of “bulldozed keel” was suggested, which means that the oceanic plate continuously transfers tectonic stresses to the front of the flat segment where the slab is already steep [26, 74, 175, 182]. In this scenario, however, the role of the thermal field and the rheology of the overriding plate in the transmission and localization of stresses is not taken into account. Martinod et al. [176] argue that the lithosphere above the flat slab is a subject to minor deformation because it is colder and stronger. Therefore, deformation localizes where the slab starts to resume its steep subduction angle, triggered by slab-pull forces, rather than where the slab is already steep, as proposed by the “bulldozed-keel” models [26, 74, 175, 182]. Our results support the hypothesis of Martinod et al. [176] in that the mafic crust and cold lithosphere beneath the Sierras Pampeanas suggest that the lithosphere there is strong. This may therefore inhibit the formation of crustal-scale faults, thus dismissing the “bulldozed keel” effect as an efficient mechanism for propagating deformation. Instead, the localization of deformation in the Sierras Pampeanas could be due to either mechanical weakening within inherited basement heterogeneities and/or increased slab-pull where the slab resumes steep subduction. In summary, the results of our study support the hypothesis that inherited rheological heterogeneities and/or slab steepening controlled the geometry of faults that delimit the spatially isolated basement uplifts and the intervening sedimentary basins of the broken foreland between 27°S and 33°S, rather than slab flattening. Discriminating between these hypotheses ultimately requires additional geodynamic modeling applied to the case of the SCA.

On the basis of geodynamic numerical modeling, it has been proposed earlier that the strength of the overriding plate sustained the flat subduction setting of the oceanic plate over the last ~20 Ma [56, 57, 59]. This has been suggested to be the case for the northern part of the SCA (27°–33°S), where the subhorizontal slab segment underlies the thick and dense crust of the foreland of the South American plate [58]. Interestingly, the cold temperatures modeled for this area in our study exhibit two effects that might also favor slab shallowing: (i) cooling of the subduction interface that enhances the coupling between the continental and oceanic plates and (ii) efficient E-W stress transmission along the cold and strong overriding plate, which forces the trench to retreat. In this context, there is a positive feedback between the cold lithosphere in the northern foreland and the flat subduction setting: The shallow slab, together with a thin radiogenic crust, would cause low temperatures in the lithosphere in this domain, thereby strengthening the overriding plate, which in turn would promote slab shallowing.

By combining conversion of S wave seismic tomography to temperatures and steady-state conductive numerical modeling, we derived the 3D lithospheric-scale temperature distribution of the southern Central Andes and adjacent forearc and foreland regions and conclude the following:

• (1)

  Distinct controlling factors of the thermal field are dominant at different depths. At shallow depths (<50 km bmsl), the thermal contrast between the warm orogen and the relatively cold areas of the forearc and foreland is modulated by the thickness of the upper radiogenic continental crystalline crust, which generates lateral changes of heat production. In the uppermost levels (< 5 km), the effect of the sediment thickness is superimposed, leading to depocenters of the foreland basins being warmer than the edges due to thermal blanketing. The cool oceanic slab outweighs the heating effect of the continental crust in regions with relatively shallow slab depth (<85 km bmsl) and where the upper continental crystalline crust is thin. This occurs in the forearc and in most of the northern part of the foreland (29°–33°S) where the slab flattens. At depths of >50 km bmsl, the spatial correlation between crustal features and thermal heterogeneities is insignificant, and the mantle heat flow and the effect of the cold slab become the main controlling factors. However, down to 100 km bmsl, temperatures are additionally affected by the radiogenic contribution of the upper continental crystalline crust in the northern part of the orogen where this unit has a significant thickness (>30 km)

There exists a strong contrast in surface heat flow between the warm orogen (80–105 mWm^-2) and the relatively cold forearc and foreland areas (40–75 mWm^-2). These variations in heat flow are primarily controlled by the thickness configuration of the uppermost layers (sediments and upper crystalline crust). The shallow cold slab affects the pattern of surface heat flow beneath the forearc

• (2)

  The modeled temperature configuration has implications for the rheology and, therefore, deformation patterns of the SCA. A cold, mafic, thick, and therefore potentially strong lithosphere beneath the broken foreland of the Sierras Pampeanas is prone to deformation processes that are controlled by inherited heterogeneities in the upper plate or by slab-pull, where the oceanic plate resumes steep subduction. In addition, such a lithospheric configuration may favor the coupling between the subducting and overriding plates, potentially contributing to a flat subduction setting

• (3)

  Sensitivity analysis of the v_s-to-T conversion shows that, at depths of <100 km, mantle composition has the strongest effect on the v_s-to-T conversion for v_s larger than 4.6 km s^-1, as is the case for the flat-slab segment. Nevertheless, compositional variations within the range of uncertainty do not significantly affect the main temperature and heat flow trends found in the tomographically studied mantle

The 3D thermal model presented in this publication can be accessed via GFZ Data Services (doi:10.5880/GFZ.4.5.2021.001).

The authors declare that they have no conflicts of interest.

This research was funded by the Deutsche Forschungsgemeinschaft (DFG) and the Federal State of Brandenburg under the auspices of the International Research Training Group IGK2018 “SuRfAce processes, TEctonics and Georesources: The Andean foreland basin of Argentina” (STRATEGy), DFG grant STR 373/34-1 to M. Strecker and M. Scheck-Wenderoth. We are grateful to Marcelo Assumpção and Mei Feng for providing the seismic tomography for the conversion to temperatures and to Corinna Kallich for the assistance on the figures. The color scales were taken from Crameri (2018), scientific color maps (doi:10.5281/zenodo.1243862). Data visualization was carried out with ParaView software (https://paraview.org). Figures were produced with Python.