Lithospheric density structure of the southern Central Andes constrained by 3D data-integrative gravity modelling

Abstract

The southern Central Andes (SCA) (between 27° S and 40° S) is bordered to the west by the convergent margin between the continental South American Plate and the oceanic Nazca Plate. The subduction angle along this margin is variable, as is the deformation of the upper plate. Between 33° S and 35° S, the subduction angle of the Nazca plate increases from sub-horizontal (< 5°) in the north to relatively steep (~ 30°) in the south. The SCA contain inherited lithological and structural heterogeneities within the crust that have been reactivated and overprinted since the onset of subduction and associated Cenozoic deformation within the Andean orogen. The distribution of the deformation within the SCA has often been attributed to the variations in the subduction angle and the reactivation of these inherited heterogeneities. However, the possible influence that the thickness and composition of the continental crust have had on both short-term and long-term deformation of the SCA is yet to be thoroughly investigated. For our investigations, we have derived density distributions and thicknesses for various layers that make up the lithosphere and evaluated their relationships with tectonic events that occurred over the history of the Andean orogeny and, in particular, investigated the short- and long-term nature of the present-day deformation processes. We established a 3D model of lithosphere beneath the orogen and its foreland (29° S–39° S) that is consistent with currently available geological and geophysical data, including the gravity data. The modelled crustal configuration and density distribution reveal spatial relationships with different tectonic domains: the crystalline crust in the orogen (the magmatic arc and the main orogenic wedge) is thicker (~ 55 km) and less dense (~ 2900 kg/m³) than in the forearc (~ 35 km, ~ 2975 kg/m³) and foreland (~ 30 km, ~ 3000 kg/m³). Crustal thickening in the orogen probably occurred as a result of stacking of low-density domains, while density and thickness variations beneath the forearc and foreland most likely reflect differences in the tectonic evolution of each area following crustal accretion. No clear spatial relationship exists between the density distribution within the lithosphere and previously proposed boundaries of crustal terranes accreted during the early Paleozoic. Areas with ongoing deformation show a spatial correlation with those areas that have the highest topographic gradients and where there are abrupt changes in the average crustal-density contrast. This suggests that the short-term deformation within the interior of the Andean orogen and its foreland is fundamentally influenced by the crustal composition and the relative thickness of different crustal layers. A thicker, denser, and potentially stronger lithosphere beneath the northern part of the SCA foreland is interpreted to have favoured a strong coupling between the Nazca and South American plates, facilitating the development of a sub-horizontal slab.

Introduction

The southern Central Andes (SCA) between approximately 27° S and 40° S formed as a result of subduction of the oceanic Nazca Plate beneath the continental South American Plate. This area is an ideal natural laboratory in which to study orogenic processes, since it is the type area for orogeny on a non-collisional, oceanic–continental convergent margin. One remarkable feature of the region is the along-strike variation in the subduction angle of the oceanic plate, which changes between 33° S and 34.5° S from sub-horizontal (< 5°) in the north to relatively steep (~ 30°) in the south (Fig. 1a; e.g., Cahill and Isacks 1992). The Andean orogen has a slightly curved shape that is convex towards the east, with major structures having a general N–S trend to the north of 33° S and an NNE–SSW trend to the south (Arriagada et al. 2013). This distinct change in structural trend is known as the Maipo Orocline (Farias et al. 2008).

Fig. 1

a Topographic elevation and bathymetric depths of the study area extracted from the ETOPO 1 global relief model (Amante and Eakins 2009), showing three subduction segments with increasing subduction angles from north to south. The extent of the forearc, the magmatic arc (red-shaded), the back-arc (grey-shaded), and the foreland are also depicted. The white rectangle shows the extent of the study area. The model covers a region measuring 700 km × 1000 km. Thick black lines indicate the main Quaternary tectonic faults (Sagripanti et al. 2017). The dashed black line shows the extent of the Maipo Orocline. On the right hand are shaded relief of the study area overlain with: b proposed sutures between terranes accreted to Gondwana during the Palaeozoic, modified from Ramos et al. 2010, and c morphotectonic domains with constituting morphotectonic provinces considered in this study: (i) forearc as white-shaded area (CC Coastal Cordillera, CD Central Depression), (ii) orogen as purple-shaded area (PC Principal Cordillera, FC Frontal Cordillera) and (iii) low-elevation back-arc and foreland regions as green-shaded area (CB Cuyo basin, EAB Extra-Andean basins, ESP Eastern Sierras Pampeanas, NB Neuquén basin, P Payenia, PrC Precordillera, TrB Triassic basins, WSP Western Sierras Pampeanas)

The continental crust of the Andean orogen (i.e., the magmatic arc and the main orogenic wedge) has been affected by both thick- and thin-skinned styles of deformation (Jordan et al. 1983). There are also systematic reductions from north to south in the amounts of crustal shortening and thickening, the topographic elevation, and the orogenic width (Kley and Monaldi 1998; Cristallini and Ramos 2000; Mescua et al. 2014; Giambiagi et al. 2015; Lossada et al. 2018). This configuration is the result of a multi-stage, tectono-magmatic evolution that has included both contractional and extensional cycles, lasting from the Neoproterozoic to the present day. While the crustal structure at shallow depths is relatively well documented, a number of unknowns remain regarding the characteristics and configuration of the deep crust and lithospheric mantle. Furthermore, the extent to which inherited variations in the composition and geometry of different layers in the upper plate may have influenced the location and style of the deformation remains poorly understood.

Seismic data from a number of different sources have been used to constrain the first-order discontinuities and seismic properties of the South American lithosphere, focusing primarily on the sub-horizontal slab region (e.g., from CHARGE, SIEMBRA, and ESP seismic experiments; Alvarado et al. 2005; Ammirati et al. 2013, 2018). These investigations were mainly used to generate data with which to construct 2D cross-sections or establish 3D models of sub-regions of the SCA. Their results have placed the Moho at an average depth of ~ 70 km below mean sea level (b.m.s.l) in the orogen and at depths ranging between 35 and 55 km b.m.s.l. in the foreland, with the depth decreasing eastwards (Alvarado et al. 2005, 2009; Heit et al. 2008; McGlashan et al. 2008; Porter et al. 2012; Perarnau et al. 2012; Ammirati et al. 2015, 2016, 2018; Pérez-Luján et al. 2015). West-to-east variations in the velocities of P and S waves have also been identified in the sub-horizontal slab region, with a relatively low-velocity zone beneath the orogen (between 69° W and 70° W) and the easternmost part of the foreland (between 65° W and 67° W), and a high-velocity zone in the back-arc and the westernmost part of the foreland (i.e., between ~ 67° W and 69° W; e.g. Ward et al. 2013; Marot et al. 2014; Ammirati et al. 2015, 2016, 2018). Most of these researchers have interpreted each of these velocity zones to reflect crustal terranes with different compositions that were accreted to the Gondwana margin during the early Paleozoic (e.g., Ramos et al. 2010). The most widely accepted of the various paleogeographic reconstructions addressing the nature of these terranes has been the Pampia–Cuyania–Chilenia model, which is based largely on surface and near-surface geological observations (e.g., Ramos et al. 1986, 2010; Ramos 2010). This model indicates the existence of two major terrane boundaries within the area that are thought to have had a direct impact on its tectonic evolution during the Cenozoic and consequently on its present-day characteristics (Fig. 1b).

In contrast to the sub-horizontal slab region, a small number of geophysical investigations to the south of 33° S have focused on the thickness of the crystalline crust (Gilbert et al. 2006; Heit et al. 2008). The results of these investigations suggest an eastward reduction in the Moho depth, from 50 km b.m.s.l. in the orogen to 40 km b.m.s.l. in the foreland. With regard to the crustal composition, possible large-scale lithological heterogeneities in the steep-slab region have yet to be investigated. Such investigations could furnish information on the extent of the crustal terranes that have been inferred to make up the overriding plate (e.g., Ramos et al. 2010). If these terranes do indeed exist and no significant overprint of their petrological compositions occurred subsequent to their accretion, compositional variations similar to those inferred for the sub-horizontal slab region leading also to variations in seismic velocity and density would also be expected. In view of these complexities and the remaining open questions, the generation of a 3D, lithospheric scale, data-integrative and gravity-constrained structural model of the SCA was considered likely to provide further insights into the deep lithospheric configuration beneath the Andean orogen. Models addressing this issue and aiming to decipher relevant lithospheric characteristics have been established for regions further north within the Central Andes (e.g., Prezzi et al. 2009; Meeßen et al. 2018; Ibarra et al. 2019), and for the oceanic Nazca Plate (Tassara et al. 2006; Tassara and Echaurren 2012). Despite these developments, a major gap remains in our understanding of the SCA.

We have therefore generated a 3D model of the lithosphere beneath the SCA and their foreland (29° S–39° S and 65°–72° W, Fig. 1a) by integrating geological and geophysical observations. This model describes the density and thickness variation of the main lithospheric layers. The modelled densities were derived from lithologic data and seismic velocities (e.g., Mescua et al. 2016; Schaeffer and Lebedev 2013). The main crustal interfaces (the Moho, top of the crystalline crust, and the base of the oceanic Nazca plate) are mostly based on gravity-independent data interpreted during earlier studies (e.g., Heine 2007; Assumpção et al. 2013). To further constrain the intracrustal structure, however, we relied on gravity forward and inversion techniques. We also assessed the sensitivity of the derived intracrustal structure to different density configurations in specific lithospheric units.

We discuss below the resulting thickness and density configurations for the crystalline crust in terms of their relationships with the spatial configuration of proposed tectonic domains and with respect to published models of accreted crustal terranes (e.g., Ramos et al. 2010). In our analysis, we first attempt to assess the extent to which the main discontinuities and compositional variations within the lithosphere relate to the tectonic events that have occurred in the region throughout its orogenic history. We also discuss the influence that the upper plate lithospheric structure has on the localization of short-term surface deformation, and on the subduction angle of the Nazca Plate.

Geologic setting

From west to east, the SCA comprise the forearc, the magmatic arc, the back-arc, and the foreland. Each of these morphotectonic provinces can be subdivided into different areas according to their Mesozoic–Cenozoic tectonic evolution (Fig. 1a; Jordan et al. 1983; Mpodozis and Ramos 1990). The forearc includes the Coastal Cordillera and the Central Depression; the magmatic arc is represented by the Eastern Principal Cordillera; the back-arc comprises the N–S oriented fold-and-thrust belts of Western Principal Cordillera, the thick-skinned structure of Frontal Cordillera, the Precordillera, and the Payenia volcanic province; the foreland includes the Sierras Pampeanas basement uplifts, the San Rafael block and the Cuyo and Neuquén basins (Fig. 1c). In this study, the magmatic arc and the main orogenic wedge (Principal Cordillera and Frontal Cordillera) are referred to as ‘Andean orogen’ (Fig. 1c).

The present-day configuration of the SCA results from a number of tectonic events that have occurred from the Neoproterozoic to the present day. Early stages in the tectonic evolution of the study area (Fig. 1) involved tectonic shortening and the accretion of crustal terranes to the proto-margin of Gondwana (e.g., Astini et al. 1995; Ramos et al. 1986, 1996, 2010; Rapalini 2005; Escayola et al. 2011). During the late Neoproterozoic-to-early Cambrian, the Pampean orogeny affected the eastern part of the study area between 27° S and 33° S, the area that is now known as the Sierras Pampeanas (e.g., Rapela et al. 1998). A further four subsequent major orogenies have been identified during the Paleozoic, these being (1) the early Ordovician Famatinian orogeny, (2) the middle–late Ordovician Ocloyic orogeny, (3) the late Devonian–early Carboniferous Chanic orogeny, and (4) the early Permian San Rafael orogeny (for reviews of these orogenies, see Heredia et al. 2018; Ramos et al. 2018; Rapela et al. 2018). The Famatinian orogeny resulted in the establishment of a magmatic arc along the western margin of Gondwana, vestiges of which still exist today in the central-western part of the Sierras Pampeanas (Ramos et al. 2010). The Ocloyic and Chanic orogenies have traditionally been associated with the accretion of the Cuyania and Chilenia exotic terranes, respectively (for reviews of these orogenies, see Ramos et al. 2010 and Giambiagi et al. 2015, and references therein). The San Rafael orogeny has been interpreted as representing a renewed onset of the subduction processes that generated a broad orogenic belt and crustal thickening along the western margin of the study area (e.g., Azcuy and Caminos 1987; Llambias and Sato 1990; Japas and Kleiman 2004).

In contrast to the compressional mountain building associated with these orogenies, the late Permian-to-early Jurassic was characterized by horizontal crustal extension associated with the break-up of Gondwana (e.g., Llambias et al. 1993). The Choiyoi volcanics were active during this time along most of South America's western margin, continuing until the early Triasssic (e.g., Ramos and Kay 1991). This magmatic episode was characterized by rhyolitic lavas and ignimbrites and volcaniclastic deposits (e.g., López-Gamundí 2006; Kleiman and Japas 2009; Spalletti and Limarino 2017). Regional horizontal crustal extension led to the formation of the Cuyo basin and other coeval basins during the Early-to-Middle Triassic, followed by the tectonic subsidence of the Neuquén Basin between Late Triassic and Early Jurassic (e.g., Uliana et al. 1989).

Subduction has been active along the Chile margin since at least the Jurassic, continuing to the present (Maloney et al. 2013 and references therein), while deformation and uplift of the SCA have been recorded since the Late Cretaceous (Fennell et al. 2015; Boyce et al. 2020). A major contractional phase occurred from the Miocene through the Quaternary (e.g., Jordan et al. 1997; Giambiagi and Ramos 2002; Giambiagi et al. 2003; Sagripanti et al. 2017; Horton 2018). From the Early Miocene to the Pliocene, the deformation migrated eastwards, forming the broken foreland of the Sierras Pampeanas (Jordan and Allmendinger 1986; Sobel and Strecker 1993; Ramos et al. 2002). The areal extent of this ongoing deformation episode coincides spatially with the sub-horizontal subduction segment of the Nazca Plate (Jordan et al. 1983), although such broken-foreland tectonism has also occurred in the thin-skinned sub-Andean foreland fold-and-thrust belt of the Santa Bárbara System (23° 30′ S–26° S, e.g., Mon and Salfity 1995). The Santa Bárbara system is still tectonically active and records the ongoing compressional inversion of Cretaceous normal faults (e.g., Mon and Salfity 1995; Kley and Monaldi 2002; Arnous et al. 2020). Finally, during the Late Miocene–Pliocene, the San Rafael block (e.g., Japas and Kleiman 2004; Folguera and Zárate 2009) and the Coastal Cordillera (Encinas et al. 2006) were subjected to horizontal crustal shortening.

Magmatic activity during the Mio-Pliocene between latitudes 33° S and 35° S was concentrated along the western slopes of the Principal Cordillera, but subsequently migrated to its current location along the border between Argentina and Chile (Fig. 1a, Stern and Skewes 1995). However, subsequent volcanism (Pliocene to Quaternary) also occurred in the present-day magmatic arc of the Principal Cordillera and in the retro-arc area of the Payenia volcanic province (Stern 2004; Folguera et al. 2009; Llambias et al. 2010), while the area between 27° S and 33° S has remained devoid of volcanic activity (Kay et al. 2006).

Methods

The workflow in this study involved two main steps: (1) data integration to set up an initial structural and density model (Section “Setup of the initial model”), and (2) gravity modelling that combined forward and inversion techniques (Section “Gravity modelling”). We first integrated geological and geophysical observations to derive the initial geometry and densities of the main lithospheric units of the SCA and their foreland. The interfaces between the modelled lithospheric units intend to represent abrupt variations in density in the lithosphere, delimiting bodies with homogenous physical properties. To further constrain our model, we tested the effect of its geometrical configuration on the gravity field. The differences between observed and predicted gravity fields reflect variations in the density or thickness that were ignored in the initial model. Since the density of the crystalline crust was considered to be homogeneous in the initial model despite indications of local variations in the physical properties (density and seismic velocity), we assumed that any major differences between observed and predicted gravity fields were indicative of density heterogeneities within the crystalline crust and, therefore, used them to obtain the horizontal density variations within this layer. For this purpose, the crust was split into two horizontal layers of different densities and the variations in the depth of the discontinuity between the two layers computed.

Setup of the initial model

The initial lithospheric model considered four major density interfaces: the ground-air/ground-water interfaces based on topography/bathymetry, the top of the crystalline crust (i.e., the base of the sediments), the Moho, and the base of the oceanic plate. To establish these interfaces, we compiled different datasets and interpolated them to a regular grid with a 25 km resolution, using the convergent interpolation algorithm of the Petrel software (Schlumberger 2019). The model covered an area measuring 700 km × 1100 km (white rectangle in Fig. 1), extending to a depth of 200 km.

Figure 2 shows the spatial distribution of the input data used to establish the initial 3D model, while Table 1 provides information on resolution and source. We used the ETOPO1 global relief model (Amante and Eakins 2009) for the topography and bathymetry. The sub-surface structures were defined by integrating information from a number of seismically constrained models, including the South American crustal thickness from the model by Assumpção et al. (2013), the sediment thickness from the CRUST1 crustal model (Laske et al. 2013), and the SLAB2 subduction zone geometry model (Hayes et al. 2018). We also included data from seismic reflection and refraction profiles across the Chilean margin (von Huene et al. 1997; Flueh et al. 1998; Araneda et al. 2003; Krawczyk et al. 2006; Sick et al. 2006; Contreras-Reyes et al. 2008, 2014, 2015; Moscoso et al. 2011), sediment-isopach maps from the ICONS intracontinental basin database (6 arc minute resolution, Heine 2007), and two oceanic sediment compilations: one along the axis of the southern part of the trench (0.05 degree resolution; Völker et al. 2013) and the other on a global-scale with 0.5 degree resolution (GlobSed, Straume et al. 2019).

Fig. 2

Distribution of input data. Sediment thickness derived from: the ICONS intracontinental basin database (Heine 2007); the CRUST1seismically constrained global model (Laske et al. 2013); the GlobSed global oceanic sediment thickness model (Straume et al. 2019) and a seismically constrained model of the sediment thickness along the southern part of the trench axis by Völker et al. (2013). Black-dashed lines indicate the location of 2D sections from seismic experiments used to constrain the top of the oceanic crystalline crust, the depth of the continental Moho, and the subduction interface close to the trench (von Huene et al. 1997; Flueh et al. 1998; Araneda et al. 2003; Krawczyk et al. 2006; Sick et al. 2006; Contreras-Reyes et al. 2008, 2014, 2015; Moscoso et al. 2011). Red circles indicate the data point constraints selected from seismic profiles that were used in the crustal thickness model of Assumpção et al. (2013) which was utilized to derive the depth of the continental Moho. The blue diamond shows the sample location for the mantle xenoliths (Conceicao et al. 2005; Jalowitzki et al. 2010). The black rectangle shows the extent of the study area. The thin black line indicates the South-American coastline

Table 1 Summary of input data used to constrain the initial 3D density and structural model of the SCA and their forelands

The uppermost layer for the oceanic part of the model is the seawater; its thickness was derived by subtracting the bathymetric depths obtained from the ETOPO1 model (Amante and Eakins 2009) from the mean sea level. The uppermost layer for the continental part of the model is the sedimentary layer; the sediment thickness was obtained from the ICONS (Heine 2007) and CRUST1 (Laske et al. 2013) data sets. The sediment thickness for the oceanic part of the model obtained from Völker et al. (2013) and GlobSed (Straume et al. 2019) and was cross-checked against available seismic profiles. The thickness of the continental crystalline crust was obtained from the model by Assumpção et al. (2013), with seismic profiles used for additional constraints close to the trench.

The thickness of the subducting section of the oceanic plate (h_s) was obtained from the SLAB2 model (Hayes et al. 2018). The thickness of the non-subducting section (h_ns) was approximated by the thickness of the thermal boundary layer, as defined by the characteristic thermal diffusion distance (e.g., Turcotte and Schubert 2002):

$${h}_{ns}=2.32\sqrt{k.t},$$

(1)

where k is the thermal diffusivity of the oceanic plate, which in this case we assigned a constant value of 8 × 10^–7 m²/s (Hasterok 2013), and t is the age of the oceanic plate obtained from a global model by Müller et al. (2008). The oceanic plate was subsequently subdivided into three units. The upper two units correspond to the basaltic and eclogitic layers of the oceanic crystalline crust, with the transformation from basalt to eclogite occurring at a depth of about 35 km b.m.s.l. (Faccenda and Dal Zilio 2017). The lower unit represents the mantle component of the oceanic slab. For the oceanic crustal layers, we assumed a constant thickness of 7 km, which is consistent with local estimates based on seismic data (e.g., von Huene et al. 1997). The thickness of the mantle component of the slab was computed by subtracting the 7 km from the total thickness of the oceanic plate (h_s and h_ns for the subducting and non-subducting sections of the oceanic plate, respectively). To calculate the true vertical thickness of each individual layer (h_v), we took into account the subduction angle α extracted from the SLAB2 model (Hayes et al. 2018) using the following equation:

$${h}_{v}=\mathrm{h}/{{\mathrm{cos}}{\alpha} }.$$

(2)

The depths to the model interfaces were computed by downward stacking of the modelled layer thicknesses. Depths to the top of the crystalline crust for both continental and the oceanic areas were obtained by adding the thicknesses of continental or oceanic sediments, respectively, to the topographic elevation or bathymetric depth. The depth to the Moho was calculated by adding the thickness of the crystalline crust to the top of the crystalline crust.

The depth of the interface between the subducting oceanic plate and the upper plate (subduction interface) was defined by interpolating the depth to the top of the oceanic crystalline crust with the depth to the top of the slab from the SLAB2 model. Seismic section interpretations were used (where available) to sort out any inconsistencies between different data sources. The depth of the oceanic Moho was computed by adding the true vertical thickness of the oceanic crust (h_v) to the depth of the subduction interface. To obtain the depth to the base of the oceanic plate, we added the thickness of the mantle component of the slab to the depth of the oceanic Moho.

Once the main lithospheric interfaces had been defined, we assigned constant densities to each of the units in the initial model, except for the lithospheric mantle. Table 1 summarizes the model layers and their assigned densities. The continental sediments were modelled with a density of 2400 kg/m³ as this has previously been estimated to be the mean density of the sediments in the Neuquén and Cuyo basins (Fig. 1c; Sigismondi 2012; Mescua et al. 2016 and references therein). The initial density assigned to the continental crystalline crust was 2800 kg/m³, which is equivalent to a felsic composition (Christensen and Mooney 1995). It also reflects the average composition of the exposed basement rocks (in the Cuyania and Chilenia crustal terranes and the Choiyoi magmatic province; see references in Mescua et al. 2016).

For the oceanic plate, we assigned a density of 2900 kg/m³ to the basaltic crustal layer, which is typical for a crust that has experienced hydrothermal alteration by seawater circulation (Stern 2002). The eclogitic crustal layer was assumed to have a density of 3200 kg/m³ (Hacker et al. 2003; Faccenda and Da Zillio 2017). For the lower unit of the slab, we assigned a density of 3340 kg/m³, which corresponds to a peridotitic mantle (Hacker et al. 2003).

The thermal field and its effect on the density distribution within the mantle can be assessed through the analysis of S-wave velocities due to their temperature sensitivity. We computed mantle densities using the SL2013sv global upper mantle tomographic model (0.5° resolution; Schaeffer and Lebedev 2013). This model is based on the inversion of surface and S waveforms for depths between 25 and 700 km b.m.s.l. At shallow depths (< 75 km b.m.s.l), the modelled velocities reflect both crustal and mantle components as a result of sampling over a depth range of 50 km. We therefore only used velocities from depths > 75 km b.m.s.l.

The S-wave velocities from the SL2013sv model were converted into densities using the Python Velocity Conversion tool (Meeßen 2017; see Supplementary Material), following a modified approach developed by Goes et al. (2000). To compute mantle densities, this method required a mantle composition, which we derived from the petrology and geochemistry of mantle xenoliths found in the Agua Poca volcanic area (Conceição et al. 2005; Jalowitzki et al. 2010; Fig. 2). These xenoliths consist mainly of anhydrous spinel lherzolite, with minor quantities of harzburgite and banded pyroxenite (Bertotto et al. 2013). A summary of the relevant petrological and geochemical analyses is provided in the Supplementary Material. The final step involved interpolation of the obtained mantle density distribution over a 3D grid (voxel) and subsequent integration of this distribution into the initial model.

Gravity modelling

We used gravity modelling to derive additional constraints on the lithospheric structure of the study area. This involved checking the gravity field calculated from the model (referred to below as the ‘calculated gravity’) against the measured gravity (referred to below as the ‘observed gravity’). For our study, we compared the calculated gravity to the free-air anomaly (observed gravity) from the EIGEN 6-C4 global gravity model (Förste et al. 2014; Ince et al. 2019) at a fixed elevation of 6 km above the datum, which had a horizontal resolution of 10 km (Fig. 3). EIGEN 6-C4 model integrates data from both terrestrial and satellite sources, and although some of the short wavelength response that a full terrestrial data set would possess main remain unresolved, it was considered appropriate for the scale of the modelling.

Fig. 3

Free-air gravity anomaly from the EIGEN 6-C4 global gravity field model (‘Observed gravity’; Förste et al. 2014; Ince et al. 2019) at 6 km above the datum used for the gravity modelling of our study. The positive values correspond spatially to the topographic highs (e.g. in the orogen, the Payenia volcanic province and the Sierras Pampeanas) and the negative values in the northern part of the study region, with the oceanic domain and the foreland basins

One of the methods used to compare calculated gravity with observed gravity relies on forward modelling. This method involves iteratively modifying the free parameters (i.e., the density or geometry of the layers) of an initial model until the calculated gravity matches the observed gravity. To forward model the free-air anomaly associated with our initial model, we used the IGMAS + software (Götze and Lahmayer 1988; Schmidt et al. 2010), which is an interactive tool for processing and integrating gravity data. The input surfaces are interfaces between different units, each with its own constant density. The software triangulates between previously defined parallel vertical sections (working sections) and the input interfaces, thus generating polyhedra that represent the three-dimensional units of the model. The algorithm then calculates the volumetric contribution of each of the units to total the gravity (the calculated gravity), by transforming the volume integral of each polyhedron to a surface integral (Götze 1984; Götze and Lahmayer 1988). The IGMAS + software allows density variations to be introduced within the various units by integrating gridded data (voxels) with the voxelization tool, which in our case was used to integrate the density values in the mantle derived from seismic velocities. The final calculated gravity is then the sum of the different contributions from the homogeneous units and the voxels. The IGMAS + software can also be used to compute the difference between the observed gravity and the calculated gravity (gravity residual), so that the two fields can be automatically compared. The user can iteratively modify the geometry or the density of the modelled units and recompute the calculated gravity, and thus the residual gravity, after each modification.

For our research, the initial model was subdivided into 69 working sections running from west to east, with a spacing of 25 km equivalent to the horizontal resolution of the model. To avoid any edge effects, we used an enlarged version of the initial model by extending the interfaces by 300 km in the N–S direction and by 500 km E–W. The voxelization tool was used to take into account any density variations within the mantle, using the previously described voxel grid. The gravity field was continued upwards to 6 km above the datum, thus allowing comparisons to be made between the calculated gravity and the observed gravity.

Because different lithospheric density models can yield similar calculated gravity fields, we aimed to minimize the number of parameters to be modified by comparing the geometries and densities of different layers with those in the previously described data sets. Since most of the parameters were substantially constrained by gravity-independent observations, we considered the residual between the observed gravity and the calculated gravity to be due to density variations within the unit for which the least information is available, i.e., the continental crystalline crust. To derive a density discontinuity within this layer, we inverted the gravity residual obtained from the forward modelling. For this purpose, the continental crystalline crust was assumed to consist of two layers with different average densities: an upper crystalline crust of 2800 kg/m³ and a lower crystalline crust of 3100 kg/m³. These values were based on the results obtained by Ammirati et al. (2015) in the sub-horizontal slab area of the SCA. These authors converted the P-wave velocities at different depths in a regional model into densities, using the empirical relationship between P-wave velocity and density reported by Brocher (2005). Inverse modelling of the gravity residual was carried out using the "Fatiando a Terra" Python tool (Uieda et al. 2013) and an approach similar to that followed by Meeßen et al. (2018). A schematic diagram of the gravity inversion procedure is provided in Fig. 4. The algorithm sets initial density perturbations (seeds) at the base of a prism mesh and propagates them until the residual gravity reaches a prescribed minimum threshold. In this study, we fixed the initial seed distribution at 500 m above the Moho and the density perturbations within the continental crystalline crust were introduced by thickening the lower crustal unit after each propagation step. The thickness was increased by 500 m between each step. To avoid non-unique solutions, the seeds were only allowed to grow upwards and were constrained to the continental crystalline crust. The resulting interface separating the upper and lower continental crystalline crustal units was subsequently integrated into the initial lithospheric configuration, to derive the final model.

Fig. 4

Schematic plots illustrating the use of gravity inversion to determine the thickness of the lower crust. Exemplary cross-sections showing the thickness of the main crustal layers down to 60 km b.m.s.l. a before and b after inversion of the initial model gravity residual. CS continental sediments, CC continental crust, UC upper continental crystalline crust, LC lower continental crystalline crust, and CM continental mantle. The observed and calculated gravity responses throughout the cross-sections are shown in solid red and dashed blue lines, respectively

Results

Initial structural model

The main characteristics of the initial 3D structural model derived from the integration of different data sets are illustrated in Fig. 5. Figure 5a shows a 3D representation of the initial model, together with the topographic elevation and the bathymetric depths extracted from the ETOPO1 model (Amante and Eakins 2009). The highest topographic elevation (~ 4.8 km a.m.s.l.) occurs in the central-northern part of the SCA. The elevation of the orogen then decreases southward down to ~ 2 km a.m.s.l. The highest topographies in the eastern back-arc and the foreland occur in the Precordillera, the Payenia volcanic province and the Sierras Pampeanas, where elevations range between ~ 2 and ~ 3.5 km a.m.s.l., with maximum elevations in the Precordillera and western Sierras Pampeanas. The forearc and most of the foreland have a maximum elevation of ~ 1.5 km a.m.s.l. The oceanic domain has a maximum depth of 4.5 km b.m.s.l., close to the trench.

Fig. 5

Setup of initial model. a 3D representation of the initial model. Abbreviations used for the main lithospheric layers are: CS continental sediments, CC continental crust, CM continental mantle, OCe eclogitic oceanic crust, MS mantle slab, OM oceanic mantle. b Depth to the top of the crystalline crust obtained from the difference between the topographic elevation and sediment thickness. c Thickness of oceanic and continental sediments obtained from the thickness maps compilations in the ICONS database (Heine 2007), the GlobSed model (Straume et al. 2019), and the seismically constrained CRUST1 model (Laske et al. 2013); SC–T Santa Clara–Tupungato depocentre in the Cuyo Basin, DC Dorso de los Chihuidos depocentre in the Neuquén basin. d Depth to the Moho from Model A (Assumpção et al. 2013) with additional seismic constraints (e.g., von Huene et al. 1997). e Thickness of the crystalline crust obtained from the difference between the depth to the top of the crystalline crust and the depth to the Moho. f Depth of the interface between the subducting oceanic plate and the upper plate coinciding with the depth to the top oceanic crystalline crust to the west of the trench and with the depth to the top surface of the subduction zone model SLAB2 (Hayes et al. 2018) to the east of the trench with additional seismic constraints (e.g., von Huene et al. 1997). The location of seismic profiles used to derive d, f are shown in Fig. 3. Boundaries of the main morphotectonic provinces are also marked: for abbreviations, see Fig. 1c

Figure 5b shows the depth to the top of the crystalline crust, which increases beneath the Neuquén and Cuyo basins to 4.5 and 5 km b.m.s.l., respectively. The maximum sediment thickness reaches ~ 5 km in the main depocentres of both of these basins (Fig. 5c), these being the Dorso de los Chihuidos in the Neuquén basin and the Santa Clara–Tupungato depression in the Cuyo basin. There is also an abrupt increase in the depth to the top of the crystalline crust from the orogen at 32° S to the Cuyo Basin at 33° S. Other basins in the foreland area of minor extent and the depths to the top of the crystalline crust do not exceed 3 km b.m.s.l.; the sediment thickness in these areas is between 2 and 3 km. In the oceanic domain, the depth to the top of the crystalline crust increases towards the trench to ~ 4.5 km b.m.s.l. and the sediment thickness reduces to less than 500 m.

The depth of the Moho is shown in Fig. 5d. A strong E–W variation is evident in the northern part of the study area where the Moho is less than 10 km b.m.s.l. in the oceanic domain, but reaches approximately 68 km b.m.s.l. within the orogen; it is about 38 km b.m.s.l. in the eastern part of the foreland. In contrast, variations in the Moho depth are less pronounced in the central and southern parts of the study area. In the central part of the area, a shallow Moho (~ 30 km b.m.s.l.) has been documented in the forearc region (71° S–72° S) and a deeper Moho (~ 50 km b.m.s.l.) within the orogen. The Moho in the southern part of the area is generally shallow (~ 30–35 km b.m.s.l.) in the Payenia volcanic province and along the western margin of the Neuquén basin, and deeper (~ 42–48 km b.m.s.l.) in the Andean orogen. The foreland in the rest of the study area has a more or less constant Moho depth of ~ 38 km b.m.s.l.

Figure 5e illustrates the thickness of the crystalline crust, obtained by subtracting the depth to the top of the crystalline crust (Fig. 5b) from the Moho depth (Fig. 5d). The thickness of the crystalline crust shows similar trends to those described above for the depth of the Moho. In the northern part of the study area, we identified a marked variation in thickness, with a thin (~ 7–10 km) crystalline crust in the oceanic domain, becoming thicker (up to ~ 70 km) beneath the orogen and then thinner again (between ~ 35 and 40 km) towards the eastern foreland. The thickness variations are less pronounced in the central part of the study area, where the thickness of the crystalline crust is between ~ 20 km and ~ 40 km beneath the forearc and foreland and ~ 50 km beneath the orogen. In the southern part of the study area, the greatest thickness for the crystalline crust (~ 50 km) occurs beneath the orogen and the least thickness (~ 20 km) beneath the Neuquén Basin. In the rest of the foreland, the thickness of the crystalline crust ranges between ~ 30 km and ~ 40 km.

Figure 5f shows the depth to the top of the oceanic crust, which is also the top of the subducting oceanic plate. The oceanic plate dips towards the east, with the subduction angle varying from north to south. In the northern part of the study area (29° S–33° S), the slab subducts sub-horizontally beneath the orogen and part of the eastern Sierras Pampeanas, and then steepens further to the east. The sub-horizontal subduction segment lies at depths of between 90 and 95 km b.m.s.l. The width of the sub-horizontal segment narrows progressively southward from 33° S to 35° S, after which the entire slab has a steeper dip. The change in subduction angle follows the westward migration of the trench from north to south (Fig. 1). The general depth of the subducting slab reaches ~ 300 km b.m.s.l. in the north-eastern part of the study area and ~ 500 km b.s.m.l in the south-eastern part.

Density distribution within the mantle

The density distribution within the mantle is dependent on both its temperature and its composition; we first assessed the influence of temperature. Figure 6a, b shows the temperature distribution at 100 km and 200 km b.m.s.l., respectively, obtained by converting S-wave velocities from the SL2013sv tomographic model (Schaeffer and Lebedev 2013). At 100 km b.m.s.l. (Fig. 6a), there is a small regional trend of increasing temperatures from north to south, ranging from ~ 1000 °C in the northernmost part of the study area to ~ 1050 °C in the southernmost part. The tomography also indicates a low-temperature zone (~ 850 °C) in the western part of the study area between 33°S and 35°S that reaches a minimum beneath the orogen. The temperature follows a similar distribution at 200 km b.m.s.l., but differs in the absolute values, which range from ~ 1200 °C in the north to ~ 1300 °C in the south (Fig. 6b). A low-temperature zone (~ 1000 °C) is again present in the west at roughly the same location as in 100 km depth distribution.

Fig. 6

Results from the conversion of shear-wave velocities of the SL2013 tomographic model (Schaeffer and Lebedev 2013) to temperatures and densities in the mantle. Temperature distribution shown for the depths of: a 100 km and b 200 km (base of the model), c Average mantle density, calculated pointwise vertically over a depth range of 75–km from the 3D voxel grid. Boundaries of the main morphotectonic provinces are also marked: for, see Fig. 1c

Figure 6c illustrates the distribution of vertically averaged density within the mantle (between 75 and 200 km b.m.s.l.), obtained by converting the SL2013sv S-wave velocities. The values range from 3340 kg/m³ in the north-western and southern parts of the study area to 3360 kg/m³ in the centre (0.5% variation). There is a northeast–southwest trend in relatively high values (~ 3356 kg/m³), which is similar to the trend observed in the temperature distributions. It is interesting to note that this low-temperature and high-density feature coincides spatially with the area in which the subduction angle of the oceanic slab changes from sub-horizontal in the north to a steeper gradient in the south (Fig. 5f); it also corresponds with the areal extent of the Maipo Orocline (Fig. 1a). However, this low-temperature and high-density feature does not extend beyond south of 38°S, where the dip of the slab is also steep, but the mantle has the lowest densities in the study area. These low densities occur beneath an area of thin crust in the vicinity of the Neuquén basin (Fig. 5e).

Forward gravity modelling

The gravity and gravity residual calculated from the initial configuration with the densities presented in Table 1 are shown in Fig. 7a, b, respectively. Figure 7c is a histogram of gravity residuals. The calculated gravity indicates positive values of up to 360 mGal between 70 and 72° W (Fig. 7a). To the east, the Precordillera, the Triassic basins, and parts of the Sierras Pampeanas all exhibit a negative anomaly down to − 180 mGal. In contrast, the Payenia province and the easternmost Sierras Pampeanas are characterized by a positive anomaly of up to 180 mGal. The calculated gravity exhibits a number of discrepancies with respect to the observed gravity field. First, there is a long-wavelength positive residual ranging from 80 to 240 mGal in the northern part of the study area, which indicates a substantial mass deficit in the model (Fig. 7b). Second, there is an isolated positive, short-wave-length residual at ~ 30.5° S–70.5° W that correlates spatially with the Elquí-Limari batholith (Nasi et al. 2010). Thirdly, negative residuals (up to − 280 mGal) are common on the western part of the study area within the oceanic domain and within the south-central parts of the forearc and the orogen (69.5° W–72° W). The back-arc area of the Payenia volcanic province also exhibits negative residuals of about − 80 mGal. These features indicate a mass excess in the initial model. In the remaining parts of the study area, the gravity residual ranges between − 30 and 30 mGal, suggesting that the initial model yields a good approximation of the observed gravity in those areas.

Fig. 7

a Calculated gravity and b gravity residual—as the difference between the observed gravity and the calculated gravity—from the initial model. Positive residuals indicate a mass excess in the model, and while negative residuals indicate a mass deficit. c Histogram and statistics (average, standard deviation σ, and minimum and maximum values) for the gravity residual from the initial model. All values are expressed in mGal. Boundaries of the main morphotectonic provinces are also marked in (a, b): for abbreviations, see Fig. 1c

Inverse modelling

We computed the depth to the top of the lower continental crust (Fig. 8a) by inverting the gravity residual of the initial model and assuming that any changes to the density configuration required for the calculated gravity to match the observed gravity would occur within the crystalline crust. The modelled depth of the intracrustal discontinuity (separating an upper crust, with a density of 2800 kg/m³, from a lower crust with a density of 3100 kg/m³) increases to ~ 40 km b.m.s.l beneath the orogen (69.5° W–71° W). In contrast, beneath the forearc (71° W–72° W), the back-arc area of the Precordillera, the northern part of the foreland (the Sierras Pampeanas), and the Neuquén basin, the depth to this intracrustal density discontinuity decreases down to 5 km b.m.s.l.

Fig. 8

a Depth to the top of the continental lower crystalline crust computed from the inversion of the gravity residual of the initial model. b Thickness of the continental lower crystalline crust calculated by adding the depth to the top of the lower continental crystalline crust to the depth to the Moho. c Gravity residual following inversion, derived from the final model. d Histogram statistics (average, standard deviation σ, and minimum and maximum values) for the gravity residual from the final model. All values are expressed in mGal. Root-mean-square error (RMSE) = 20 mGal. Boundaries of the main morphotectonic provinces are also marked in (a–c): for abbreviations, see Fig. 1c

Figure 8b shows the thickness of the lower continental crystalline crust layer, obtained by subtracting the depth to the top of the lower continental crystalline crust from the Moho depth. The modelled lower crystalline crust in the northern part of the study area is relatively thick, ranging between 30 and 50 km beneath the Precordillera and the western Sierras Pampeanas, and parts of the northern forearc and the Neuquén basin also have a thick lower crystalline crust (~ 25 km). In contrast, a thin lower crystalline crust was modelled for the south-central part of the orogen (~ 17 km), the volcanic province of Payenia (~ 15 km), and the oceanic domain (~ 5 km). In the remaining foreland area, the lower crystalline crust has a thickness of between ~ 20 and 25 km.

Through the integrated forward and inverse modelling, we derive the final model. The residual gravity trends of this final model are shown in Fig. 8c,d. The final model accurately reproduces the observed free-air anomaly over most of the study area (RMSE = 20 mGal). Nevertheless, short-wave-length (< 20 km) residuals remain in the north-central part of the orogen, possibly due to intracrustal heterogeneities at a scale that the model is unable to resolve.

We also computed the density distribution for the entire continental crystalline crust by averaging the densities of the upper and lower layers, taking into account the spatial variations in their thicknesses. Figure 9 presents the resulting average-density distribution, overlain with the boundaries of the morphotectonic provinces (as in Fig. 1c). Areas with high average densities (> 3000 kg/m³) in the continental crystalline crust correlate spatially with the north-central part of the back-arc Precordillera fold-and-thrust belt, the basement uplifts of the western Sierras Pampeanas, the north-eastern Sierras Pampeanas, and the foreland basins of Cuyo and Neuquén. In contrast, areas with lower average crystalline crustal densities (~ 2900 kg/m³) correlate with the oceanic domain, the Principal Cordillera, the Frontal Cordillera, the southern part of the Precordillera, and the Payenia volcanic province. Intermediate average crustal densities (between 2950 and 3000 kg/m³) occur in the narrow part of the southern Coastal Cordillera, the forearc basin of the Central Depression, and the remaining foreland areas (the south-eastern Sierras Pampeanas, the San Rafael Block, and the sedimentary basins away from Andes).

Fig. 9

Average crustal density of the total (upper and lower) crystalline crust as derived from the final model. Boundaries of the main morphotectonic provinces are also marked: for abbreviations, see Fig. 1c

Figure 10a, b depicts the gradient of (i) topographic elevation and (ii) average crustal density, also showing those crustal seismic events with a moment magnitude of 4 or more (International Seismological Centre 2020) and the geodetically derived orientations of contractional and extensional strain rates (Drewes and Sánchez 2017). Although the geodetically derived data cover our entire study area, only those sites with the highest horizontal strain rates (> 5.10^–8 y⁻¹) are shown. A clear spatial correlation is evident between high horizontal strain rates (Drewes and Sanchez 2017) and the steepest topographic gradients in the study area (Fig. 10a). A spatial correlation is also observed between the horizontal strain rate and areas with abrupt variations in the average crustal density (Fig. 10b). Most of the large magnitude (≥ 4) seismic events also occur in areas with marked variations in the average crustal density.

Fig. 10

Gradient of a topography and b average crustal density overlain with: seismic events of moment magnitude ≥ 4 (small magenta squares; International Seismological Centre 2020) and geodetically derived orientations of contractional and extensional strain rate (Drewes and Sanchez 2017). The lengths of the red and blue lines indicate the magnitude of the maximum horizontal strain rate. Only the sites with the highest maximum horizontal strain rate (> 5e⁻⁸ y⁻¹) are shown. Blue lines indicate compression. Red lines indicate extension. Boundaries of the main morphotectonic provinces are also marked: for abbreviations, see Fig. 1c

Sensitivity analysis

The chosen densities for the final model are between specific error ranges (e.g., Chirstensen and Mooney 1995). We, therefore, analysed the sensitivity of the gravity field to different density configurations for the continental sediments, the continental crystalline crust, and the mantle. We first established alternative configurations to those of the initial model by imposing a density perturbation on the continental sediments, the continental crystalline crust, and the mantle, one at a time. As a second step in the sensitivity analysis, we carried out a gravity inversion (as in Section “Inverse modelling”) for each configuration to compute the individual effect of each density perturbation on the resulting thickness of the lower crystalline crust. We also calculated the residual gravity following the inversion for each alternative configuration, to determine the most robust configuration, i.e., the configuration that most closely matched the observed gravity field. Since the gravity field is dependent on the density, depth, and mass of the rock units (e.g., Turcotte and Schubbert 2002), this type of analysis allows the effect that density perturbations have on the resulting lower crustal thickness to be quantified. We were thus able to determine whether or not the gravity residuals did indeed represent intracrustal heterogeneities or whether they simply resulted from an erroneous attribution of densities to the other lithospheric units.

First, for the continental sediments, we selected densities of 2300 and 2600 kg/m³, based on the maximum and minimum densities encountered for the Cuyo and the Neuquén basins (Mescua et al. 2016). For a density of 2300 kg/m³, the thickness of the lower continental crystalline crust increases by up to 1.5 km (Fig. 11a) compared to the final model. In contrast, for a sediment density of 2600 kg/m³, the lower continental crystalline crustal unit is up to about 3 km thinner than in the final model (Fig. 11b). These variations occur mainly beneath the Neuquén and Cuyo basins (Fig. 11a, b). For both of these densities, the inversion enabled us to match the observed gravity field within an RMSE of 21 mGal (Fig. 12a, b). Nevertheless, negative residuals (< − 30 mGal) remained along the north-western border of the study area, possibly indicating the presence of intracrustal heterogeneities that cannot be resolved at the scale of the model.

Fig. 11

Variation of the thickness of the continental lower crystalline crust with respect to the final model induced by varying the density in selected modelled units: a continental sediments = 2300 kg/m³; b continental sediments = 2600 kg/m³; c continental upper crystalline crust = 2750 kg/m³; d continental upper crystalline crust = 2850 kg/m³; e continental lower crystalline crust = 3000 kg/m³; f continental lower crystalline crust = 3200 kg/m³; g homogeneous mantle = 3360 kg/m³; h density distribution within the mantle obtained through conversion of S-wave velocities from the SEMUM2 tomographic model (French et al. 2013). Positive values indicate a thinner lower crust than in the final model. Boundaries of the main morphotectonic provinces are also marked: for abbreviations, see Fig. 1c

Fig. 12

Gravity residual (difference between observed and calculated gravity) from models with density variations in selected modelled units after inversion: a continental sediments = 2300 kg/m³; b continental sediments = 2600 kg/m³; c continental upper crystalline crust = 2750 kg/m³; d continental upper crystalline crust = 2850 kg/m³; e continental lower crystalline crust = 3000 kg/m³; f continental lower crystalline crust = 3200 kg/m³; g homogeneous mantle = 3360 kg/m³; h density distribution within the mantle obtained through conversion of S-wave velocities from the SEMum2 tomographic model (French et al. 2013). RMSE (root-mean-square error) in mGal. Boundaries of the main morphotectonic provinces are also marked: for abbreviations, see Fig. 1c

Second, another important predefined feature in the modelling process was the density of the continental crystalline crust. We investigated the effect of varying the density of the upper layer by ± 50 kg/m³ (1.7%). Modelling the lowest density in that range for the upper continental crust yielded resulted in a maximum increase of 9 km in the thickness of the lower continental crystalline crust in the northern part of the orogen (Fig. 11c). In contrast, there was a reduction in thickness of the lower continental crystalline crust of 1.5 km in the oceanic domain (Fig. 11c). Modelling the highest density in the range for the upper continental crust reduced the thickness of the lower continental crystalline crust to 8 km, as opposed to an increase in thickness of the same layer of 1.5 km in the oceanic domain (Fig. 11d). The models with reduced or increased upper continental crustal densities were able to accurately reproduce the gravity field following inversion, with RMSEs of 19 mGal and 22 mGal, respectively (Fig. 12c, d). However, residual anomalies in the northern part of the study area were predominantly > 30 mGal, indicating that the alternative models are not able to accurately reproduce the observed gravity in this part of the study area (Fig. 12c, d).

Third, we investigated the effects of varying the density of the lower continental crystalline crust by ± 100 kg/m³. Reducing the density of this lower crystalline unit resulted in an increase in its thickness of up to 17 km in the northern part of the orogen (Fig. 11e). In contrast, increasing the density by the same amount resulted in a lower crust that was up to 12.5 km thinner (Fig. 11f). With regard to the gravity residuals of these last two configurations, the inversion yielded a less accurate match with the observed gravity than the final model and the other alternatives described above in this section (RMSE = 33 mGal for the less dense lower continental crystalline crust and 25 mGal for the denser lower unit; Fig. 11e, f). The final model (Sect. 4.3) was therefore preferred over these alternative configurations.

The fourth step in the sensitivity analysis involved investigating the effect of varying the density distribution in the mantle. We first considered a constant mantle density of 3360 kg/m³. Following inversion, the lower crust was up to 3.5 km thicker along the western part of the study area than in the final model (Fig. 11g). This area has the largest density differences between the preferred model and the model with a constant mantle density. In addition, the resulting lower crust was up to 2 km thinner in the southern part of the study area. With regard to the gravity residuals, the inversion yielded an accurate match with the observed gravity field (RMSE = 21 mGal; Fig. 11g), but negative residuals (< − 30 mGal) remained in the north-western part of the study area.

Finally, we modified the mantle density configuration by considering a different tomographic model than the SL2013sv tomography model used for the final configuration, in which S-wave velocities were converted to densities (the SEMum2 model, French et al. 2013). This model has a coarser horizontal resolution (2°) than the SL2013sv tomography (0.5° resolution), but it can be used to investigate the effect that the mantle density has on the gravity field, and therefore also on the resulting intracrustal structure. The resulting lower crust was up to 2 km thicker in the western part of the study area and up to 1 km thinner in the southern part (Fig. 11h). This alternative configuration also yielded an accurate match with the observed gravity, with similar values to those in the final model (RMSE = 21 mGal; Fig. 11h).

Table 2 summarizes the results of the sensitivity analysis for different densities in the various units. To compare the different outcomes, we calculated a factor f by dividing the induced variations in lower continental crustal thickness by the percentage that parameters were altered by to induce these variations. This factor measures how the lower continental crustal thickness changes in response to variations in the predefined densities of the modelled units. This approach indicates that the intracrustal structure is more sensitive to density perturbations in the mantle (f = 5.8 km/%) than to those in the crystalline crust (f =  ~ 5.3 km/%). The areas in which the largest variations in thickness occur generally have a thick crust, e.g., beneath the northern part of the orogen and the adjacent foreland (Fig. 11e, f). This would imply that perturbations in the average crustal density are relatively small (< 18 kg/m³) compared to the average density for the entire crustal column. The major trends in the resulting density distribution would therefore persist, even in scenarios with different upper and lower crustal densities. In contrast, the densities of the sediments do not have any significant effect on the thickness of the lower crust, or therefore on the resulting crustal-density distribution (f = 0.37 km/%).

Table 2 Summary of results from sensitivity analysis

On the basis of these effects resulting from variations in the input data, we conclude that the main spatial trends in the depth to the top of the lower crust would persist despite changes to the input data, within the range of natural variability.

Discussion

Our results indicate that the crystalline crust in the orogen is thick (Fig. 5e, 40–70 km) and has a low density (Fig. 9, ~ 2900 kg/m³), which suggests a felsic composition. Although data on the crustal composition beneath the SCA are scarce, seismic surveys to the north of 36° S have indicated low seismic velocities beneath the orogen compared to beneath the adjacent foreland (e.g., Ward et al. 2013; Marot et al. 2014; Ammirati et al. 2016). Our model agrees with these observations in that the average density of the crystalline crust beneath the orogen is lower than beneath most of the foreland (Fig. 9, ~ 3000 kg/m³), which suggests compositional differences between the two areas.

The mechanisms that produced an overall less dense crust beneath the orogen are likely to have been related to plutonic additions with felsic-to-intermediate compositions, associated with the various Andean-type magmatic arcs that were first established along the western margin of the study area at 580 Ma and that are currently represented by the Cenozoic Andean volcanic chain (e.g., Stern and Skewes 1995; Kay et al. 2006). Mescua et al. (2016) have also suggested that the upper levels of the crystalline crust in the orogenic interior have relatively low densities. This area of low crustal densities correlates spatially with the vestiges of the extensive extensional Choiyoi magmatic province (Fig. 13a; Mescua et al. 2016; Spalletti and Limarino 2017), which was responsible for the extrusion of felsic magmas and the intrusion of granitic material during the Late Permian–Triassic (e.g., Mpodozis and Kay 1990; Llambias and Sato 1990; Cristallini and Ramos, 2000; Giambiagi and Ramos 2002). Moreover, alternating cycles of extension and contraction could have led to a depletion in mafic components as a result of crustal delamination (e.g., Kay et al. 1994; Allmendinger et al. 1997). Other driving mechanisms associated with the Cenozoic Andean orogeny could also have been responsible for the lower average densities in the crust beneath the orogen. These may have involved the stacking of upper crustal rocks by thrusting (e.g., Ramos et al. 2002) and/or the underthrusting of cratonic material, possibly introduced from the east and bringing less dense material into the deep crust (Armijo et al. 2015). This last mechanism is, however, deemed less likely to have occurred in the SCA, since there is no geophysical evidence for westward-directed thrusting of the cratonic material beneath the Andes between 29° S and 39° S (e.g., Alvarado et al. 2005, 2009; Heit et al. 2008; McGlashan et al. 2008; Porter et al. 2012; Perarnau et al. 2012; Ammirati et al. 2015, 2016, 2018; Pérez-Luján et al. 2015). The low average crustal density beneath the orogen may be characteristic of ocean-continent subduction scenarios, as has previously been suggested in other gravity modelling studies of the Central Andes (e.g., Tassara et al. 2006; Ibarra et al, 2019). In contrast, collisional orogens on continent–continent margins usually possess a dense mafic lower crust, as observed in the Urals, the Himalayas, and the European Alps (e.g., Stadtlander et al. 1999; Kumar et al. 2019; Spooner et al. 2019).

Fig. 13

a Crystalline crustal-density distribution and morphotectonic units overlain with the extent of the Permo-Triassic Choiyoi magmatic rocks (modified from Mescua et al. 2016 and Spalletti and Limarino 2017). b Average shear-wave velocity from ambient-noise tomography in the Central Andes (Ward et al. 2013) overlain with morphotectonic units boundaries. The white areas are areas with no data coverage (e.g., oceanic domain, Neuquén basin, and Extra-Andean basins). c Crystalline crustal density distribution overlain with proposed lines of suture between different terranes accreted to Gondwana during the Palaeozoic (modified from Ramos et al. 2010). Thick grey lines show the extent of the ophiolitic belt to the west (OB) and the Valle Fértil Lineament to the east (VFL). For abbreviations of morphotectonic province names, see Fig. 1c

To the east of the Andean orogen, the remainings of the back-arc and foreland exhibit marked variations in the thickness (Fig. 5e) and average density (Fig. 9) of the crystalline crust. Three distinctive domains stand out: (1) a northern domain (29° S–34.5° S) with a thick crust (40–60 km) and high densities (~ 3000–3050 kg/m³); (2) a southern domain (37° S–39° S) with a thinner crust (~ 20 km) and high densities (~ 3050 kg/m³); and (3) a central domain (34.5° S–37° S) with intermediate crustal thicknesses (35–45 km) and low-to-intermediate densities (~ 2900–2950 kg/m³).

The high average crustal densities in the northern part of the foreland (including most of the Sierras Pampeanas) and the north-central sector of the Precordillera (29° S–34.5° S) suggest a mafic crustal composition adjacent to the orogen. This is compatible with the results from a number of field investigations that have concluded that the basement rocks are mainly gabbros and leuco-gabbros (e.g., Vujovich and Kay 1998; Ramos 2010; Pérez-Luján et al. 2015). Our results suggest that the middle and lower crustal levels also have mafic compositions. The higher crustal densities in this part of the foreland compared to the orogen interior could be associated with partial eclogitization of the lower crust, as has been previously suggested on the basis of geophysical studies, in particular in view of the increase in V_p and the V_p/V_s ratio beneath the foreland (Gilbert et al. 2006; Alvarado et al. 2007, 2009; Ammirati et al. 2013, 2015, 2018; Marot et al. 2014). The factors causing partial eclogitization in the lower crust remain unclear (e.g., Ammirati et al. 2015), although some studies suggest that this process could be facilitated by small amounts of water released from the sub-horizontal section of the subducted oceanic crust (Leech 2001; Yang et al. 2008).

In the eastern part of our study area, high crustal densities (3000–3100 kg/m³; Fig. 9) predominate in parts of the Sierras Pampeanas. Isostatic processes and the overall higher density crust in this morphotectonic province could explain why the elevation is relatively low in this particular area (Allmendinger et al. 1990; Bellahsen et al. 2016) despite its large crustal thickness (~ 50–60 km, Fig. 5e). However, the style of deformation, which involves isolated reverse-fault bounded uplifts of crystalline blocks, may also provide an explanation for the lack of crustal stacking and topographic growth. It has been proposed on the basis of the surface geology that the eastern Sierras Pampeanas are underlain by a felsic crust (e.g., Ramos et al. 2010 and references therein). However, our own results suggest that, apart from the shallower crustal layers, the average crustal composition in this part of the Sierras Pampeanas is likely to be mafic. Nevertheless, mafic rocks are not the only materials in the lower crust that can attain densities above 3000 kg/m³. For example, according to Hacker et al. (2011), rocks with high proportions of felsic components such as pelites and wackestones can attain high densities when metamorphosed under extreme P/T conditions (such as 1050 °C at 3 GPa). These authors referred to the addition of such highly metamorphosed sediments to the base of the crust as ‘crustal relamination’. Possible lithologies with this potential are those of the metasedimentary rocks belonging to the Puncoviscana Formation of north-western Argentina (Turner et al. 1960), as proposed by Meeßen et al. (2018). These rocks have been interpreted to have been deposited in the Puncoviscana ocean along the western Gondwana margin between 540 and 524 Ma (Escayola et al. 2011). Although outcrops of the Puncoviscana Formation are only found to the north of 27° S (e.g., Turner 1960), it has been proposed that correlative sedimentary rocks extended further south (Escayola et al. 2011). Indirect evidence for the existence of these sedimentary rocks within our study area may possibly exist in the Ordovician igneous rocks of the Sierras Pampeanas (Otamendi et al. 2009; Ducea et al. 2015; Alasino et al. 2016). The results of petrological and isotopic analyses suggest that the main petrogenetic processes operating at that time involved metasedimentary materials. Nevertheless, no unambiguous evidence has been found to confirm that these metasedimentary rocks did indeed originated in the passive-margin basins of Gondwana. The corroboration of such a hypothesis would require isotopic evidence from rocks that are equivalent to the Puncoviscana Formation; unfortunately, no such data are currently available from within the study area.

In contrast to most of the Sierras Pampeanas, the south-eastern part of the Sierras Pampeanas has a low average crustal density (2960 kg/m³, Fig. 9). This is in agreement with the seismic results and outcrop observations from the area, which indicate that the underlying crystalline crust has a felsic composition (e.g., Perarnau et al. 2012; Alvarado et al. 2009; Ammirati et al. 2018). These felsic rocks may have originated from the intense magmatic activity that shaped the eastern calc–alkaline magmatic belt of the Pampean orogen during Neoproterozoic–Cambrian times (Ramos 2010 and references therein).

The crystalline crust in the southern part of the foreland (37° S–39° S, which correlates spatially with the Neuquén basin) is thin (~ 20 km, Fig. 5e) and has a high average density (~ 3050 kg/m³, Fig. 9). This area is likely to have been affected by magmatism associated with the Choiyoi extensional event (Fig. 13a), leading to crustal thinning and underplating of mafic material during basin initiation (Kay et al. 1989; Llambias et al. 1993), as has been interpreted for the magmatic Kenya Rift (KRISP Working Party 1991). These past processes could account for the present-day crustal configuration of the Neuquén basin. Despite the Neuquén and Cuyo basins having similar sediment thicknesses, our model suggests a substantial difference in the configuration of the crystalline crust beneath these basins (Fig. 5c). In contrast to the thin, dense crust beneath the Neuquén basin, the crust beneath the Cuyo basin is relatively thick (50–40 km, Fig. 5e) and less dense (~ 2950–3000 kg/m³, Fig. 9). The modelled difference in the configuration of the crystalline crust could relate to differences in the main process that drove the evolution of the two basins, even though the evolution of these depocentres was in both cases initiated by rifting and thermal subsidence. This was followed by a partial compressional reversal of the normal faults and flexural subsidence due to the load of the adjacent Andean orogen (e.g., Giambiagi et al. 2008; Barredo 2012). The thin, dense crystalline crust beneath the Neuquén basin may indicate that the main processes involved in the generation of accommodation space were the earlier rifting processes, rather than the subsequent flexure in response to shortening and thrust-stack loading during the Andean orogeny. In contrast, the relatively thick and light crystalline crust beneath the Cuyo basin may indicate that the flexural processes had a larger influence on the generation of basin subsidence than the previous extensional processes.

There is a distinct area of low average crustal density (~ 2900 kg/m³, Fig. 9) evident within the central part of the back-arc region adjacent to the orogen (34.5° S–37° S) that corresponds to the Payenia volcanic province. The low crustal densities compared to surrounding areas are likely to be related to the felsic effusive products from the Choiyoi volcanics (Fig. 13a). The density distribution may also have been affected by regional magmatic and thermal processes, since intraplate volcanism has been active in the region for the last 20 Ma (Kay et al. 2006; Folguera et al. 2006; Llambias et al. 2010).

Our results for the area to the west of the SCA suggest that the crystalline crust beneath the forearc comprises two domains: a northern domain of high average density (~ 3000 kg/m³; Fig. 9) and a southern domain of slightly lower average density (2950 kg /m³, Fig. 9). The presence of intermediate-to-mafic rocks from Jurassic–Cretaceous magmatism (Parada 1990) could explain the high crustal densities in the northern domain. These crustal densities are also in agreement with the results from regional geophysical investigations where P-wave crustal velocities of up to 6.5 km/s have been interpreted to indicate crystalline crustal rocks (Contreras-Reyes et al. 2014, 2015). In the southern part of the forearc, the spatial extent of the relatively low modelled crustal densities correlates well with the distribution of exposed late Paleozoic andesitic and medium-grade metamorphic rocks (Muñoz et al. 2000; Hervé et al. 2013), and it is consistent with the P-wave velocities of ≤ 6 km/s modelled for the entire crustal column of the area (Krawczyk et al. 2006; Contreras-Reyes et al. 2008). Toward the oceanic domain, the reduction in average crustal density modelled along the entire forearc may relate to eroded or highly tectonized volcanics (Contreras-Reyes et al. 2014, 2015).

On a regional scale, there is a rough spatial correlation between the modelled crustal density distribution (Fig. 9) and average S-wave velocities derived from ambient-noise tomography of the Central Andes (Fig. 13b; Ward et al. 2013). Areas of lower crustal densities generally have lower seismic velocities, while areas with higher crustal densities crust have higher velocities. However, this pattern does not apply to the Payenia volcanic province. A careful interpretation of the seismic velocities beneath this area is therefore required, specially in view of the fact that the small number of seismic recording stations in the area means that they are only able to provide limited coverage and relatively few seismic records, which compromises the reliability of the seismic results (see Fig. 1 of Ward et al. 2013).

Our study has revealed that the orogen, the southern part of the forearc, the Payenia volcanic province, and the southern part of the Sierras Pampeanas are areas with low crustal densities and seismic velocities, which suggests that they may all be underlain by crust of mainly felsic composition. In contrast, zones with high densities and seismic velocities are likely to have a more mafic crust, as is the case for the northern part of forearc and most of the foreland. The low crustal density and low seismic velocity areas are therefore likely to have higher concentrations of radiogenic elements than their mafic counterparts (e.g., Vilà et al. 2010) and these felsic areas may therefore contribute significantly to the heat budget of the lithosphere, making it hotter and therefore probably weaker and more susceptible to deformation processes. In contrast, a more mafic crust is likely to be cooler and to be able to withstand higher stress levels before reaching the P–T window for brittle or viscous deformation (e.g., Burov 2007).

Our results show that compositional variations in the overriding plate, such as those we derived in our study, need to be taken into account when investigating the rheological behaviour of the lithosphere and associated deformation patterns in tectonically active areas adjacent to a convergent plate margin. In this context, our analysis has shown that the geodetically derived distribution of high horizontal strain rates and high-magnitude seismic events correlates with areas that have high topographic gradients (Fig. 10a) and marked variations in gravity-constrained average crustal densities (Fig. 10b). This suggests that short-term (i.e., earthquake-related) deformation is likely to be localized in these areas. These spatial correlations are probably related to an increase in gravitational potential energy (GPE) due to the high topographic gradients and abrupt variations in crustal density. Part of the increased GPE would be transferred into the elastic strain field, ultimately facilitating seismogenic processes (Barrows and Langer 1981; Barrows and Barrows 2010).

Another question regarding density contrasts and deformation processes is whether the relationship between deep-seated density contrasts and strain distributions also applies to the long-term deformation features of the SCA. The spatial correlations between elevation, crustal thickness, and long-term deformation in a sector of the Central Andes north of the area of our study (21–29° S) have been extensively discussed by Ibarra et al. (2019) in terms of variations in lithospheric strength and GPE. They argued that the high elevation and large crustal thickness of the orogen results in high a GPE that prevents internal deformation within the orogen, and that the deformation front is therefore propagated to areas with lower GPE, i.e., cold, strong areas such as the foreland. The first-order characteristics of the lithosphere in the sector of the Central Andes investigated by Ibarra et al. (2019) are similar to those in our study area, with felsic compositions, high elevations, and a thick crust beneath the orogenic part of the area and with more mafic compositions, lower elevations, and a thinner crust beneath the northern part of the foreland and the Neuquén basin. This suggests that similar mechanisms to those proposed for the Central Andes between 21° and 29° S are likely to have been active in the SCA, and that pre-strained areas with different densities and rheological characteristics may focus deformation under changing geodynamic boundary conditions. This has also been suggested, for example, for the early tectonic evolution of the Eastern Cordillera of Argentina north of latitude 26°, and for regions farther to the west that are now part of the Andean Plateau (the Puna; Montero Lopez et al. 2020). The thicker, more rigid basement of the Puna region was less affected by the initial Paleogene Andean deformation processes, while the sedimentary rocks exposed in the Cordillera Oriental experienced pervasive deformation. The zone between these two morphotectonic provinces and the areas of early uplift and deformation coincides with the north–south oriented Ordovician Eastern Puna Eruptive Belt (Bahlburg 1990; Bahlburg and Hervé 1997), where basement heterogeneities such as metamorphic fabrics and shear zones were preferentially reactivated during Andean mountain building (Hongn and Riller 2007; Hongn et al. 2010; Montero-López et al. 2020). However, additional rheological and structural studies focusing on the localization of strain over geological time-scales will be required to further test the possibility that these principles also apply to the SCA and other mountain belts.

Geodynamic numerical modelling of convergent non-collisional margins has suggested that the strength of the overriding plate influences the angle of subduction of the oceanic plate (e.g., Manea et al. 2012; Sharples et al. 2014). Upper plates in a subduction system that are thick, mechanically strong, and moving towards the trench therefore favour sub-horizontal subduction due to their strong plate coupling (Sharples et al. 2014). This could be the case for the northern part of the SCA (27° S–33° S), where a gently dipping slab has been identified beneath the thick, dense, and potentially strong foreland of the South American Plate. In view of these considerations, the configuration of the overriding plate could indeed contribute to the maintenance of a shallow subduction angle in the northern part of the SCA.

We also investigated the thermal field and composition of the lithospheric mantle. Temperatures and densities converted from S-wave velocities show a heterogeneous distribution, with a northeast–southwest low temperature and high-density trend at all depths (Fig. 6a–c). We interpret these lower temperatures and higher densities as being due to sampling of the slab, which is situated farther east in the northern part of the foreland and farther west in the central part of the foreland because of the N–S variation in the subduction angle (Fig. 5f). It remains unclear, however, why this anomaly is not visible in the southern sector of the foreland.

Finally, our results clearly indicate that the density distributions in the crystalline crust and the lithospheric mantle are not spatially correlated with terrane boundaries inferred mainly from surface observations, such as those of the Chilenia–Cuyania–Pampia terranes (Fig. 13c; Ramos et al. 2010 and references therein). Although we observe heterogeneities in the average crustal densities, the boundaries between the different density domains do not generally coincide with those proposed for shallow crustal depths. The only exception to this observation is to the north of the Cuyania Terrane (29° S–33° S), where a sharp transition between areas of different crustal density coincides with the proposed terrane boundaries. These horizontal density contrasts correlate spatially with the loci of major reverse faults such as that responsible for the Valle Fértil Lineament in the east (e.g., Ramos et al. 1986), and with an ophiolitic belt in the west (Haller and Ramos 1984), interpreted as the Cuyania–Pampia and Cuyania–Chilenia boundaries, respectively (Fig. 13c). The near-surface expression of these features terminates at about 33 or 34° S (Alvarez et al. 2012; Boedo et al. 2016) and their southward continuation remains poorly constrained. If, however, these boundaries continue farther south, the individual crustal compositions of each terrane are likely to have been overprinted by tectonic processes that occurred subsequent to their accretion to Gondwana. An alternative scenario is that the crustal heterogeneities of each terrane stem from tectonic processes that occurred prior to their accretion to Gondwana. In line with this latter interpretation, previous investigations (e.g., Ramos et al. 1998) have argued that the crustal composition of the Cuyania terrane is heterogeneous on the basis of differences in basement composition between the Precordillera and the westernmost Sierras Pampeanas, which would have to be older than the proposed terrane collision during the Middle-to-Late Ordovician. It is therefore also plausible that other pre-collisional crustal heterogeneities might extend farther south. Nevertheless, from the relationship between crustal density domains and tectonic events that occurred after the terrane accretion discussed above, we favour the hypothesis that the crustal heterogeneities developed after terrane accretion. Indeed, continuous crustal domains sutured in deformation phases that occurred prior to the Andean orogeny appear to act as coherent assemblages. A possible explanation for this observation is that the dominant process controlling the present-day configuration of the SCA lithosphere is the subduction of oceanic crust at the convergent margin. Nevertheless, further investigations are required before a robust conclusion can be reached, including, for instance, high-resolution seismological experiments to constrain P- and S-wave velocity variations within the crust and the use of additional constraints from gravity gradients. Such detailed studies would help to clearly define the terrane boundaries if present.

Conclusions

By integrating different geological and geophysical datasets, including seismic, gravity, and petrological data, we have derived a 3D structural and density model describing the first-order characteristics of the lithosphere beneath the southern Central Andes (SCA), their forearc, and their foreland. Our findings suggest that:

• 1. The crystalline crust is thicker (~ 55 km) and less dense (~ 2900 kg/m³) beneath the orogen (magmatic arc and main orogenic wedge) than beneath the forearc (~ 35 km, ~ 2975 kg/m³) and foreland (~ 30 km, ~ 3000 kg/m³). The foreland comprises three principal crustal domains: (1) a thick, dense northern domain (40–60 km, ~ 3000–3050 kg/m³); (2) a thin, dense southern domain (~ 20 km, ~ 3050 kg/m³); and (3) a central domain with intermediate crustal thickness (35–45 km) and low-to-intermediate crustal densities (~ 2900–2950 kg/m³). These variations are likely to reflect differences in the long-term tectono-magmatic evolution of each domain.

• 2. The average-density configuration within the crystalline crust correlates spatially with the location of major seismic events and the present-day horizontal surface strain rates. The highest horizontal surface strain rates occur in the areas with the highest topographic gradients within the orogen, and across regions with abrupt variations in crustal density. This suggests that crustal compositions and associated rheological variations exert an important influence on short-term deformation patterns.

• 3. North of 33° S, the thick crust of the orogen and the thick and potentially strong crust of the foreland may enhance the coupling between the subducting and overriding plates, which could facilitate the establishment of a sub-horizontal subduction setting.

• 4. There is generally no spatial relationship evident between the boundaries proposed for the Pampia–Cuyania–Chilenia terrane model and boundaries based on either the lithospheric density distribution or the surface-deformation patterns. If these boundaries were originally present, they must have been subsequently obliterated by tectono-magmatic processes following the accretion of these terranes to the margin of Gondwana. These results also suggest that, in the SCA, the terrane provenance has little influence on the present-day deformation.

• 5. The northeast–southwest trending low-temperature, high-density feature modelled for the lithospheric mantle reflects the perturbation of the thermal field by the cool subducting slab of the Nazca Plate.