The Central Andean orogen formed as a result of the subduction of the oceanic Nazca plate beneath the continental South-American plate. In the southern segment of the Central Andes (SCA, 29°S-39°S), the oceanic plate subducts beneath the continental plate with distinct dip angles from north to south. Subduction geometry, tectonic deformation, and seismicity at this plate boundary are closely related to lithospheric temperature distribution in the upper plate. Previous studies provided insights into the present-day thermal field with focus on the surface heat flow distribution in the orogen or through modelling of the seismic velocity distribution in restricted regions of the SCA as indirect proxy of the deep thermal field. Despite these recent advances, the information on the temperature distribution at depth of the SCA lithosphere remains scarcely constrained.

To gain insight into the present-day thermal state of the lithosphere in the region, we derived the 3D lithospheric temperature distribution from inversion of S-wave velocity to temperature and calculations of the steady state thermal field. The configuration of the region – concerning both, the heterogeneity of the lithosphere and the slab dip – was accounted for by incorporating a 3D data-constrained structural and density model of the SCA into the workflow (Rodriguez Piceda et al. 2020a-b). The model consists on a continental plate with sediments, a two-layer crust and the lithospheric mantle being subducted by an oceanic plate. The model extension covers an area of 700 km x 1100 km, including the orogen (i.e. magmatic arc, main orogenic wedge), the forearc and the foreland, and it extents down to 200 km depth.

To predict the temperature distribution in the SCA, the model volume was subdivided into two domains: (1) a shallow domain, including the crust and uppermost mantle to a depth of ~50 km below mean sea level (bmsl), where the steady-state conductive thermal field was calculated using as input the 3D structural and density model of the area (Rodriguez Piceda et al., 2020a-b); (2) a deep domain between a depth of ~50 and 200 km bmsl, where temperatures were converted from S wave seismic velocities (Assumpção et al., 2013) using the approach by Goes et al. (2000) as implemented in the python tool VelocityConversion (Meeßen 2017).

The 3D model of Rodriguez Piceda et al. (2020) consists of the following layers: (1) water; (2) oceanic 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 sub-lithospheric mantle. For the computation of temperatures in the shallow domain, three main modifications were made to the 3D model of Rodriguez Piceda et al. (2020a-b). First, we removed the water layer thus considering the topography/bathymetry as the top of the model. Second, the horizontal resolution was increased to 5 km and, third, the layers were vertically refined by a factor of 3 to 32.

We assigned constant thermal properties (bulk conductivity λ and radiogenic heat production S) to each layer of the model according to each lithology (Alvarado et al. 2007, 2009; Ammirati et al. 2013, 2015, 2018; Araneda et al., 2003; Brocher, 2005; Čermák and Rybach, 1982; Contreras-Reyes et al., 2008; Christensen & Mooney, 1995; Gilbert et al., 2006; Hasterok & Chapman, 2011; He et al., 2008; Marot et al., 2014, Pesicek et al., 2012; Rodriguez Piceda et al., 2020; Scarfi & Barbieri, 2019; Vilà et al.,2010; Wagner et al., 2005; Xu et al., 2004).

The steady-state conductive thermal field in the shallow domain was calculated applying the Finite Element Method as implemented in the software GOLEM (Cacace & Jacquey, 2017; Jacquey & Cacace, 2017). For the computation, we assigned fixed temperatures along the top and base of the model as thermal boundary conditions. The upper boundary condition was set at the topography/bathymetry and it is the temperature distribution from the ERA-5 land data base (Muñoz Sabater, 2019). The lower boundary condition was set at a constant depth of 50 km bmsl for areas where the Moho is shallower than 50 km bmsl and at the Moho depth proper where this interface is deeper than the abovementioned threshold. The temperature distribution at this boundary condition was calculated from the conversion of S-wave velocities to temperatures (Assumpção et al., 2013).