Crustal structure of the Lazufre volcanic complex and the Southern Puna from 3-D inversion of magnetotelluric data: Implications for surface uplift and evidence for melt storage and hydrothermal fluids Open Access

Studies of subduction zones have shown that there are significant changes in the crustal and upper-mantle structure of these plate boundaries along strike. Some of the most important variations are the spatial and temporal changes in the angle of the subducting plate (Kay and Coira, 2009). These changes in slab angle have been shown to have a profound impact on the evolution of the continental crust above the subducting plate. These along-strike variations are obvious in the Central Volcanic Zone of the Andes, which has been formed by convergence between the Nazca and South American plates (Fig. 1). Of all the subduction zones around the world, this one is unique because the plate convergence has formed a thickened crust in the Altiplano and the Puna that locally exceeds a thickness of 70 km (Beck et al., 1996; Yuan et al., 2002; Heit et al., 2014; Ryan et al., 2016). This crustal thickening has occurred by a combination of processes that include magmatism, shortening, and perhaps mass redistribution by crustal flow (Husson and Sempere, 2003). The Central Andes are also characterized by changes in subduction angle in both space and time (e.g., Cahill and Isacks, 1992). Cessation of volcanism occurred as the slab flattened in the past, and episodes of intense magmatism (ignimbrite flare-ups) have occurred as the slab steepened (e.g., Coira et al., 1993; Kay and Coira, 2009). The distribution of partial melt in the crust is significant because it profoundly impacts the strength of the crust (Unsworth and Rondenay, 2012). Oblique convergence has caused along-strike migration of regions of flat-slab subduction (Allmendinger et al., 1997 and references therein), allowing a spatial study to investigate inferred temporal changes.

An enigmatic aspect of Andean magmatism was the discovery of regions that were undergoing persistent uplift at rates in excess of 10 mm per year extending over areas greater than 50 × 50 km (Pritchard and Simons, 2004). In several cases, these uplifts were not associated with active volcanic centers. It has been suggested that these locations could be places where pluton formation was occurring today, giving the opportunity to study this important geological process.

Detailed geological and geophysical studies of the Central Andes have contributed to our knowledge of the present-day structure of this region and provided important constraints on the temporal evolution. Previous geophysical studies have been focused on seismic exploration (Beck et al., 1996; Myers et al., 1998; Beck and Zandt, 2002; Bianchi et al., 2013; Heit et al., 2014). Geophysical methods that measure electrical resistivity provide an alternate view of structure and are especially sensitive to the presence of fluids, including partial melts. Although a number of magnetotelluric (MT) transects have been reported in the Central Andes (Fig. 1), the spatial coverage is far from uniform, and is much more limited than the seismic studies using surface waves (Ward et al., 2017).

In this paper, a new 200-km-long MT transect across the Southern Puna plateau in Argentina at a latitude of ~25°S is presented. The location is shown in Figures 1 and 2. The new data were collected as part of the PLUTONS project with the goal of addressing the two main questions: (1) to define the shallow crustal structure in a region that is undergoing a phase of increasing magmatism, perhaps in response to recent slab steepening and lithospheric delamination; and (2) to investigate the processes that cause surface uplift of relatively large regions that are not associated with centers of active volcanism. A list of participants in the PLUTONS project and details of funding can be found in Pritchard et al. (2018).

The Puna plateau is located at the southern limit of the Central Volcanic Zone between 22°–28°S and has higher elevations (>4000 m), rougher topography, and is narrower than the Altiplano to the north (Allmendinger et al., 1997). The southern limit of the Puna is defined by the southward change from normal subduction to the Pampean flat slab (Cahill and Isacks, 1992; Hayes et al., 2018). The Puna can be divided into distinct northern and southern regions. This study is focused on the Southern Puna in both Argentina and Chile.

The Southern Puna is characterized by several salars, which are depressions containing sediment with high salt content, and which may, or may not, be associated with the evaporation of a lake. Compared to the Northern Puna, the Southern Puna has a number of distinct features that include (1) greater relief (Isacks, 1988), (2) the pattern of Mio-Pliocene faulting with more normal faulting in the Southern Puna (Marrett et al., 1994), (3) presence of mafic lavas, and (4) less crustal shortening than in the Altiplano. It is also the location of the large ignimbrite deposits formed by the eruption and collapse of the Cerro Galan caldera, which is one of the largest on Earth with dimensions of 25 × 35 km (Sparks et al., 1985). Folkes et al. (2011) provide a summary of ages from 2.13 Ma to 2.06 Ma with the youngest reliable age reported from Ar/Ar dating by Kay et al. (2011). Based on these observations, Kay et al. (1994) proposed that piecemeal lithospheric delamination had occurred in the Pliocene beneath the Southern Puna. This argument can be used to explain the high elevation, associated with a thinner crust and shallow lithosphere-asthenosphere boundary (LAB) beneath the central part of the Southern Puna plateau. The delamination hypothesis was supported by additional data assembled by Heit et al. (2014).

The basement of the Southern Puna is mainly composed of crystalline rocks of primarily Precambrian and Ordovician age, which are overlain by younger sedimentary rocks of Cenozoic age deposited in a number of basins (Allmendinger et al., 1997; Gorustovich et al., 2011). These rocks have been deformed by compression that led to (1) more irregular and higher topography than the Altiplano; (2) crustal thickening with the thickest crust on both sides of the plateau and thinnest in the center; and (3) inversion of some of the sedimentary basins (Riller and Oncken, 2003). Beginning in the Miocene, the arc volcanism of the Central Volcanic Zone developed with several distinct pulses of activity that produced a range of extrusive igneous rocks, including both extensive lava flows and pyroclastic deposits. Compositions range from felsic ignimbrites to andesitic-dacitic volcanos and smaller mafic systems to the east (Coira et al., 1993). Whereas the modern volcanic arc strikes approximately north-south at this latitude, the crust of this region is characterized by several lineaments that trend primarily northwest-southeast (e.g., Salfity and Gorustovich, 1998). As suggested by many authors, these lineaments are believed to be long-lived, deep-seated features that appear to control the distribution of magmatism and hydrothermal ore deposits, but their depth extent and surface expressions are still under debate. In the study area of this paper, the Archibarca lineament extends southeast from Volcán Llullaillaco (near the Argentina-Chile border) to the eastern edge of the Puna where Cerro Galan is located (Richards and Villeneuve, 2002).

The crustal and upper-mantle structure of the Southern Puna has been investigated by a range of seismic studies. The most detailed of these was the 74-station Puna-delamination array recorded from 2007 to 2009, and the resulting data have been analyzed by a number of techniques. Teleseismic receiver function studies showed that the Southern Puna has a crustal thickness of 50–55 km, thickening to 60–75 km beneath the volcanic arc (Heit et al., 2014). These crustal thicknesses are somewhat less than those observed in the Northern Puna and Altiplano. This study also reported a v_p/v_s ratio less than 1.70, which showed that the crust was felsic and also a region with a high v_p/v_s ratio beneath Cerro Galan beginning at 10 km depth below sea level; this region was inferred to be a zone of partial melting. Using the same data set, Bianchi et al. (2013) described a large, low-velocity body beneath Galan and Antofalla volcanoes extending all over the Southern Puna and named it the Southern Puna Magmatic Body (SPMB). Attenuation tomography showed that the Southern Puna had a more heterogeneous structure, with stronger attenuation than the Northern Puna (Whitman et al., 1992). More recently, Liang et al. (2014) investigated seismic attenuation and reported a zone of low Q from 25°–27.5°S and high Q zones dipping from both western and eastern margins, and they interpreted this as evidence for delamination. The most recent analysis of this data set used local earthquake tomography and gave the most reliable image to date of the low-velocity feature beneath Cerro Galan, indicating the continued presence of partial melt in the crust from a depth of 10 km below sea level to the Moho. This study also confirmed the presence of a high-velocity mantle feature west of Cerro Galan that was interpreted as a block of delaminated lithosphere (Chen et al., 2020). The lithosphere below the Puna is estimated to be 80–100 km thick (McGlashan et al., 2008; Heit et al., 2014), giving a very thin lithospheric mantle, which is less than below the Altiplano where the lithosphere is more than 150 km thick. This thin lithospheric mantle provides further evidence for lithospheric detachment and thinning (Schurr et al., 2006), which numerical modeling has shown can produce surface uplift of several kilometers and significant crustal heating within a short amount of time on the order of 1–10 m.y. (Krystopowicz and Currie, 2013; Comeau et al., 2021).

The crustal structure of the Puna has been placed in a regional context with a large compilation of surface-wave dispersion data that now covers a significant proportion of the Central Andes (Delph et al., 2017; Ward et al., 2017). While this method has reduced depth resolution compared to the seismic tomography described above, it has the advantage of excellent spatial coverage in the crust. At a latitude of 25°–26°S in the Southern Puna, these studies show the presence of two distinct, low shear-velocity zones in the crust (1) east of Volcán Lastarria and (2) around Cerro Galan (Delph et al., 2017; Ward et al., 2017) as shown in Figure 2B. Although low, the shear velocities do not decrease to the very low values observed below the Altiplano where they are associated with the Altiplano-Puna Magma Body (APMB) at 22°S beneath Volcán Uturuncu. The pattern of crustal velocities described by Ward et al. (2017) supports ideas of a southward-propagating pulse of magmatism associated with a steepening of the Nazca plate that was first proposed by Kay and Coira (2009). Recently, a full waveform inversion tomography (Chen et al., 2020) detected similar low-velocity zones to those described by Ward et al. (2017).

The volcanic arc at this latitude is characterized by a 30-km-long chain of Quaternary volcanic centers near the Chile-Argentina border that strike N15°E and are often referred collectively as Lazufre, named from Volcán Lastarria and Cordon del Azufre (Pritchard and Simons, 2004; Ruch et al., 2008). The chain has produced andesitic to dacitic lavas and extends from Cerro Bayo in the south to the 5700 m Volcán Lastarria in the north. Present-day activity is characterized by persistent fumarolic activity and sulfur flows at Volcán Lastarria (Froger et al., 2007; Aguilera et al., 2016). A remarkable aspect of this area is the uplift with decadal time scales at a rate of 30 mm/year that was initially detected by INSAR data (Pritchard and Simons, 2004) and more recently investigated with surface geodetic data. Deformation studies placed the inflation source in the depth range 8–13 km below the surface (Perkins et al., 2016 and references therein). An analysis that considered uplift between 2003 and 2010 found no evidence for temporal changes in either the location or shape of the inflation source. The inflation source was best represented as a flat-topped body at a depth of 2–14 km below sea level (Remy et al., 2014). Since there have been no eruptions at this location in the historical record, it is inferred that the surface uplift could be due to the intrusion of magma at depth, leading to either a future eruption or the formation of a pluton (Pritchard et al., 2018).

The PLUTONS project was initiated to investigate the crustal structure at Volcán Lazufre and also at another major zone of surface deformation at Volcán Uturuncu at 22°S in Bolivia (Pritchard et al., 2018) with locations shown in Figure 1. The two locations have some similarities, but also some significant differences. At Uturuncu, the uplift is axisymmetric and surrounded by a zone of subsidence. Lazufre has an elliptical uplift pattern and no outer region of subsidence. At Uturuncu, the inferred inflation source is in the depth range 18−35 km below sea level and is located within a major crustal magma body—the Altiplano Puna Magma Body (APMB). The Lazufre uplift is localized in the middle or upper crust and not associated with a major magma body. Despite the uncertainties involved in data analysis, the inferred center of uplift at Uturuncu is significantly deeper than at Lazufre (Pritchard et al., 2018).

Seismic observations in the Puna have also provided evidence of the movement of magma in the crust. Mulcahy et al. (2014) reported seismic swarms located (1) under the Cerro Galan caldera and (2) south of the Peinado volcanic center, which is located south of Lazufre. These seismic swarms suggest the movement of magma in the crust.

Magnetotelluric (MT) exploration uses low-frequency electromagnetic waves to image the electrical resistivity of the Earth and provides complementary information about crustal and upper-mantle structure to that obtained from seismic studies. As shown in Figure 1, a growing number of MT surveys have been undertaken in the Central Andes and have shown that the resistivity structure varies significantly along strike, implying a variation in the melt and/or fluid content of the crust and mantle (Schwarz and Kruger, 1997; Lezaeta et al., 2000; Brasse et al., 2002; Brasse and Eydam, 2008; Díaz et al., 2012; Comeau et al., 2015, 2016; Kühn et al., 2017; Slęzak et al., 2021). The new profiles in southern Peru (P1 and P2) shown in Figure 1 were collected from 2017 to 2019 (Unsworth et al., 2020).

In the Puna region, a number of previous MT studies have been reported (Fig. 2A) around the zone of surface uplift at Lazufre on the border between Chile and Argentina. Budach et al. (2013) collected an array of 14 long-period MT stations around Volcán Lastarria and reported a deep-seated conductor that dipped beneath the Puna to the east. The size of the MT array limited the resolution of this feature beneath the Southern Puna in Argentina. A more detailed study around Lastarria used an array of broadband MT stations to image the upper-crustal magmatic system (Díaz et al., 2015). A 3-D inversion of these MT data revealed a shallow conductor that was interpreted as a magma body. The deeper conductor imaged by Budach et al. (2013) was also detected, but again, the array was not able to determine the crustal structure beneath the Southern Puna to the east.

A regional-scale MT study of the Puna at this latitude in the Andes has not yet been undertaken. In this paper, the first regional-scale MT study of the Southern Puna is presented. New broadband MT data were collected by the University of Alberta in 2013, and data were combined with the two previous MT studies to give the most complete image possible of the crustal structure beneath Lazufre and the Southern Puna.

The data used in this study were collected in three separate campaigns. The station distribution is shown in Figure 2A. An initial data set of 15 long-period stations were collected by the Free University of Berlin in 2010 (Budach et al., 2013). A more detailed broadband MT data set was collected around Volcán Lastarria by the University of Chile from 2013 to 2015 using Metronix ADU-07 instruments at 30 stations (Díaz et al., 2015). The coverage was extended to the east in Argentina by the University of Alberta using Phoenix Geophysics V5-2000 instruments to record broadband MT data at 38 stations in October and November 2013 on a transect that extended east from Volcán Lastarria to Cerro Galan across the Southern Puna. The recording time at each station was from 18–36 h resulting in data in the frequency range 1000–0.0003 Hz. Data from this broadband MT survey were of mixed quality at frequencies below 0.001 Hz, primarily because of the short recording times and low-resistivity material at the surface that resulted in weak electric fields. Some stations recorded just electric fields; therefore, processing of these stations used the magnetic field data from an adjacent station, typically located a few kilometers away. Each survey was limited by the logistical difficulties of travel in this region of the Andes. Each MT data set was processed using standard techniques associated with robust statistics and remote reference processing to obtain the highest quality estimates of the impedance tensor and magnetic field transfer functions.

Different instruments were used to collect the three MT data sets, thus resulting in three distinct frequency sets. The efficient 3-D inversion of the data required that they were combined into a single frequency set using two stages of editing. In the first stage, each station was edited using the original frequency set, with outliers being excluded. The data were then interpolated onto a common data set with ten frequencies per decade in the range 100–0.0001 Hz. The merged MT data were examined for problems arising from interpolation, and a second round of editing was undertaken.

The data are illustrated in Figure 3 as pseudosections and in Figure 4 as sounding curves at selected stations. The pseudosections show the apparent resistivity components (ρ_xy and ρ_yx) in geographic coordinates with the x-direction defined as north and the y-direction as east. Both components of apparent resistivity are heavily influenced by galvanic effects, which are frequency-independent perturbations of the surface electric field arising from near-surface resistivity structure. This results in the vertical striping seen in the apparent resistivity pseudosections. Despite this complication, a common pattern is that the high-frequency data are characterized by relatively low apparent resistivity. At lower frequency, a range of apparent resistivity values is observed. The phases show a clearer pattern with both components (φ_xy and φ_yx) showing high values in the frequency range 100–1 Hz, indicating a low-resistivity surface layer. From 1–0.01 Hz, the phase values are lower, indicating a deeper, higher resistivity layer. Below a frequency of 0.01 Hz, significant variations in the phases occur along the profile, with different patterns observed in φ_xy and φ_yx. Low values of φ_yx occur west of the volcanic arc and around Volcán Lastarria (L). High φ_yx values (>45°) are observed east of the volcanic arc and extend as far as the Salar de Antofalla (SA). Farther east, a region of high values of φ_yx is observed close to Cerro Galan (CG). Values of φ_yx decrease to 45° at the eastern end of the profile. The values of φ_xy show a different pattern at frequencies below 0.01 Hz. Values close to 45° are observed from the west end of the profile to the Salar de Antofalla. Values slightly greater than 45° are observed close to Cerro Galan.

In Figure 4A, a long-period MT sounding west of the volcanic arc is plotted and shows the moderate phase values described above. Figure 4B shows a broadband MT station close to Volcán Lastarria. Note that at the lowest frequencies, the apparent resistivity is not decreasing, and the phase values are ~45°, which indicates the absence of a shallow conductor directly beneath the volcano. Figure 4C shows a broadband MT station to the east and the decreasing apparent resistivity, and high phases indicate a major low-resistivity feature at the greatest depths measured by the data. Figure 4D shows a broadband MT station just east of the Salar de Antofalla, and apparent resistivity and phase both show that the deep structure is more resistive than to the west of the salar. Finally, Figure 4E shows a broadband MT station at the eastern end of the profile, close to Cerro Galan, and the decreasing apparent resistivity and high phases both indicate the presence of a strong conductor.

The vertical magnetic fields were used to compute the transfer functions that are also plotted in Figures 3 and 4. These quantities can also be plotted as induction arrows as shown in Figure 5. In the so-called Wiese convention, the real component of the arrows points away from regions of high conductivity (low resistivity). At frequencies in the range 1–0.1 Hz, the real induction arrows show a lot of scatter, reflecting heterogeneous near-surface resistivity structure (Figs. 5A and 5B). At a lower frequency of 0.01 Hz (Fig. 5C), a more consistent pattern emerges. West of the volcanic arc, the arrows point N50°W, indicating the presence of a conductor between the volcanic arc and the Salar de Antofalla. A second conductor is located east of the Salar de Antofalla and north of Cerro Galan by a reversal of the induction arrows at a frequency of 0.01 Hz. The quadrature, or imaginary, components are illustrated in Figure S1¹ and show a similarly complex behavior, giving further evidence for a 3-D resistivity structure beneath the study area.

The induction arrows give strong evidence for shallow 3-D resistivity structure, because of scatter of arrows at frequencies greater than 0.01 Hz. The arrows at frequencies below 0.01 Hz have a consistent orientation ~N60°W and perpendicular to the strike of the Salar de Antofalla.

The analysis of the MT data in the previous section gives evidence for a relatively complex subsurface resistivity structure beneath the Southern Puna. The collection of MT data along a profile was motivated by the plan to use a 2-D inversion to interpret the data. Before proceeding with a 2-D inversion, the dimensionality was investigated using two approaches.

In the first approach, the phase tensors were analyzed to provide a graphical way to visualize the dimensionality of the MT data (Caldwell et al., 2004). The orientation of the major axes of the phase tensor ellipse defines the regional strike, and the ellipticity defines azimuthal variations. A 1-D resistivity structure will result in a circle, while 2-D or 3-D resistivity structures will result in ellipses. The phase tensors are plotted in Figure S2. At frequencies greater than 1 Hz, the phase tensors are generally circles with small skew values, indicating a shallow, locally 1-D resistivity structure. From 1–0.01 Hz, a more complex pattern is observed with highly elliptical phase tensors and high skew magnitudes (>4°) showing 2-D and 3-D resistivity structures with no preferred strike direction. Below 0.01 Hz, a more consistent pattern emerges with the major and minor axes oriented parallel or perpendicular to N60°W. This pattern is consistent with the induction vectors described above. Note the 90° change in ellipse orientation at longitude 67°W indicating the location of a discrete conductor—in the same location as the induction arrow reversal.

In a second approach, tensor decomposition was applied to the MT data using the algorithm of McNeice and Jones (2001). This approach showed that the strike direction varied along the profile, making the application of a 2-D analysis challenging. The problems with a 2-D approach were confirmed by the inconclusive results of some 2-D inversions implemented with the algorithm of Rodi and Mackie (2001). Based on this analysis, it was decided to focus on a 3-D inversion approach for the MT data.

The MT data were inverted using the ModEM algorithm (Kelbert at al., 2014). In MT data analysis, inversion is a mathematical process that finds a 3-D resistivity model that fits the measured MT data to a specified statistical level. A key challenge of inverting geophysical data is that a range of Earth models can be found to fit a given data set. This non-uniqueness arises from (1) the physics of the geophysical methods used to study the Earth and (2) the presence of noise in the data. A common approach to address this non-uniqueness is to apply conditions such as spatial smoothness to the resistivity model. As a consequence, the 3-D inversion of the Puna MT data required a significant amount of experimentation to identify an inversion strategy and MT data subset that gave a stable inversion result. This section describes the preferred inversion model. The other 3-D inversions included a range of starting models, coordinate systems, error floors, and smoothing parameters. The effect of the MT frequency band inverted was also investigated.

All inversions included topography that was discretized using a 100 m vertical grid size. The average elevation across the measurement locations was more than 4200 m. The horizontal grid size was chosen to be 2 km in order to find a compromise between a fine grid and reasonable computation time. Thus, a number of stations around Volcán Lastarria were removed from the analysis. The complete MT data set at Lazufre was analyzed by Díaz et al. (2015) to give a model of the upper-crustal structure around the volcano. In contrast, this paper develops a resistivity model of the crust across the entire Puna. Outside of the core region, the horizontal grid size increased geometrically by a factor of 1.3. Beneath the layers used to model the topography, the grid size increased geometrically in the vertical direction by a factor of 1.1. The total grid contained 70 × 152 × 96 cells.

The starting model resistivity was found to have a significant effect on the final model. Models with a starting resistivity greater than 100 Ωm remained relatively resistive and were unable to fit the high values of φ_yx at low frequency observed east of the volcanic arc and north of Cerro Galan. Inspection of the apparent resistivity curves in the range 100–10 Hz gave a median value of 10–20 Ωm. Inversions that began from a half space with a resistivity in this range were able to fit the low-frequency MT data quite well and gave an acceptable overall misfit.

The effect of the low-resistivity ocean can be seen at a significant distance inland, at a period that increases with distance from the coast. The MT data used in this study were examined for the so-called coast effect, but no evidence was found that the data were influenced in this way by the Pacific Ocean.

The MT data were selected at 21 frequencies in the band 100–0.0001 Hz and included both the full impedance tensor and tipper. The impedance data were assigned an error floor of 5% on each of the diagonal and off-diagonal elements and 0.02 on the tipper. Many inversions were performed to determine a 3-D resistivity model that could be considered robust, and two of these models are illustrated in this paper. Two additional models are shown in the Supplemental Material (^{footnote 1}), along with the fit to the measured MT data.

The first inversion (puna_lazufre_60a) began from a 10 Ωm half space and used a covariance value of 0.3 in all three directions. The initial damping factor was set to λ = 10 and the inversion terminated when this parameter was reduced to λ = 10⁻⁸ after 207 iterations. A range of damping factors was investigated, and λ = 10 gave similar models to those that used higher values, while reducing the overall computation time. The root mean square (r.m.s.) misfit was reduced from 7.50 to 1.60. The model is shown in Figure 6 in both horizontal map view slices and vertical sections. The vertical section for the regional profile is taken as segments between sites on the profile (a fence section), so that only the best constrained parts of the resistivity model are displayed. Note that in all figures of the inversion models, the vertical axis shows depth below sea level. The surface elevation in the study area is typically in the range 4–5 km. The fit of the data in pseudosection format is illustrated in Figure S3b for this inversion. The misfit residuals are shown in Figure S4a.

A second inversion (puna_lazufre_61c) began from a starting model with a resistivity of 10 Ωm for the crust and mantle of the South American plate, and 1000 Ωm for the Nazca plate. The slab depth was determined from Hayes et al. (2018). The initial damping factor was set to λ = 10 and the inversion terminated when this parameter was reduced to λ = 10⁻⁸ after 218 iterations. The r.m.s. misfit was reduced from 7.54 to 1.63, and the model is shown in Figure 7. The fit of the data in pseudosection format is illustrated in Figure S3c for this inversion. The misfit residuals are shown in Figure S4b.

Above the Nazca plate, both 3-D resistivity models are quite similar. The MT data are relatively insensitive to the high resistivity expected for the subducting Nazca plate, as confirmed later in this paper by the sensitivity analysis. Note that the presence of the high-resistivity slab in Figure 7 forced the inversion to place a lower resistivity in the mantle wedge for the second inversion.

The features observed in both Figures 6 and 7 are listed in Table 1, and the cause of the low resistivity will be interpreted in later sections.

The 3-D resistivity models shown in Figures 6 and 7 were evaluated in a number of ways to determine the degree to which the main model features were required by the measured MT data.

Firstly, a range of inversions was performed to investigate the sensitivity of the final resistivity model to the inversion control parameters. A key aspect of this procedure was to compare the inversion model obtained by inversion of just the impedance data with the model obtained from the combined inversion of the impedance and vertical magnetic field (tipper) data. These models are compared in Figure S5 where it can be seen that the major conductors (C3, C4, and C5) are well imaged in the vertical sections for both inversions. Feature C3 was located directly beneath Lazufre in the impedance-only (Z) inversion but was imaged slightly to the east in the impedance plus tipper inversion (Z + T). The geometry of the dipping conductor C5 is largely similar in both models, which is important because this feature represents a major crustal feature and may represent the pathway of melt supply to the upper crust beneath the Lazufre inflation region. The resistivity models are different away from the area of the profile. This is as expected because the impedances are sensitive to structures below the profile, while tipper data are sensitive to structures that are off-profile. Other inversion parameters investigated included the following: error floors, starting model resistivity, azimuth of coordinate system, etc. Each of these inversions produced resistivity models that differed in some ways from the model presented in Figure 6. However, all of the main features remained. This exercise showed that the model in Figure 6 did not depend on a narrow choice of inversion parameters.

In the second stage, the 3-D resistivity model was edited, and selected regions of the model were made resistive. The forward response was computed, and the change in misfit was analyzed. Details are provided in Figures S6 and S7 (^{footnote 1}).

In a third stage, a synthetic MT inversion approach was used to investigate model resolution, and the results are illustrated in Figures S8 and S9.

Overall, these tests showed that the models in Figures 6 and 7 were resolved to depth of 50 km below sea level along most of the profile. Resolution was locally deeper in the Lazufre area where long-period MT data were available, and the crustal resistivity was higher than east of the volcanic arc beneath the Southern Puna. As expected, resolution is highest beneath the regions where MT stations are located. Feature C5 is located 10 km south of the main profile, but the data are sensitive to this major conductor. It should also be noted that this feature was detected by independent analyses of subsets of the MT data by Budach et al. (2013) and Díaz et al. (2015). The sensitivity test discussed previously shows that it affects the yx data component, since the east-west electric currents will be strongest in an east-west–oriented conductor.

The inversion models presented in Figures 6 and 7 and in the Supplemental Material show a number of common features that are independent of the choice of inversion control parameters or the data selected for inversion.

The near-surface structure has a variable resistivity that can be correlated with regions of volcanic rocks and internally drained salars.

The upper crust to a depth of ~20 km below sea level is generally resistive (>100 Ωm) and can be interpreted as dry sedimentary and igneous rocks in the upper crust. Within this layer, a number of discrete conductors are observed.

1. Beneath the Salar de Antofalla, a low-resistivity pathway (C1) extends through the upper crust. This feature is resolved most clearly in inversions that used both the impedance and tipper data.

2. Adjacent to Cerro Galan, the resistivity in this depth range is ~1 Ωm (C2).

3. Southeast of Lastarria, a low-resistivity feature (C5) can be observed dipping eastward from 5 km to 25 km below sea level in the horizontal slices shown in Figures 6 and 7.

The lower crust from 20 km below sea level to the Moho has a low resistivity in the range 3–30 Ωm. This layer has some variability with the lowest resistivity observed (1) around Cerro Galan in feature C4 and (2) in the region from the Salar de Antofalla to the Lazufre region as feature C3. West of the volcanic arc, the lower-crustal resistivity is more sensitive to the choice of data used in the inversion. Inversions that used just the impedance data extend feature C3 westward in the entire lower crust (Fig. S3c [^{footnote 1}]). In contrast, inversion of both impedance and tipper data shows that the low-resistivity layer dips to the Moho and may extend into the mantle wedge (Figs. 6 and 7).

The sensitivity tests in Figures S6, S7, and S8 showed that the MT data have limited sensitivity to upper-mantle structure. This is as expected for broadband MT data with the lowest frequency of 0.001 Hz in a relatively low-resistivity region. Collection of high-quality, long-period MT data would be required to image the upper mantle in this region.

The strongest crustal conductor detected in this study is the feature C2 at 67°W and located adjacent to the Cerro Galan caldera. The low resistivity of this feature most likely reflects a zone of partial melt and/or aqueous fluids. This feature (C2) is in the depth range 5–10 km below sea level with a resistivity of 1–3 Ωm. It should be noted that this feature is located to the north of Cerro Galan caldera in the model and that the MT data may not be able to determine if it extends beneath the caldera because the feature is some distance from the profile.

To distinguish between the various explanations of low resistivity, information about the melt or brine composition is needed. The Cerro Galan ignimbrites have been analyzed in significant detail by a number of petrological studies that have constrained the spatial and temporal evolution of the magma that caused this eruption (Sparks et al., 1985; Folkes et al., 2011). Kay et al. (2011) presented additional analyses of the ignimbrite composition and proposed that melt had evolved at three distinct levels within the crust. The deepest was in a crustal melt zone close to the Moho at a depth of 50 km. Melt then ascended to mid-crustal levels (20–30 km below sea level) where they formed a crystal mush zone. As eruptions took place, melt would have accumulated in transient upper-crustal magma bodies at 4–8 km depth below sea level.

The resistivity models in Figures 6 and 7 show evidence of each of these magma bodies. The deepest feature close to the Moho is not well resolved by the MT data presented in this paper, since the overlying conductors obscure it, a common problem in investigations of low-resistivity volcanic environments (Comeau et al., 2015).

Feature C4 occurs in the depth range 20–30 km below sea level and is regionally quite extensive. This feature can be identified with the intermediate depth magma body described by Kay et al. (2011). The Cerro Galan magmatic system was also investigated by Delph et al. (2021) with a combination of seismic and geochemical observations. Their synthesis was consistent with the petrological work of Kay et al. (2011). The resistivity is in the 3–10 Ωm range, and petrological data can be used to estimate the melt fraction required. In previous studies, the compilation of Pommier and Le Trong (2011) was used for this type of model interpretation. In this study, published studies for a specific composition are used to avoid problems arising from interpolation. Kay et al. (2011) report a silica content in the range 68%–71%, Na content of 3%, and an estimated water content of 2.5–4 wt%, although values greater than 5% are possible. The body is inferred to be at a temperature of 800–850 °C. With the surface at 5 km above sea level, the depth range 25–35 km below the surface corresponds to pressures in the range 625–875 MPa, with an average of 750 MPa. The equation of Guo et al. (2016) for felsic melts is valid for these conditions. At 850 °C, a water content of 3.25% and pressure of 750 MPa predict a melt resistivity of 1.5 Ωm (Fig. 8A). If water content varies from 2.5% to 4% and pressure ranges from 625 MPa to 875 MPa, this simple analysis predicts a range of melt resistivity from 2.39 Ωm to 0.94 Ωm. To account for the 3–10 Ωm observed in the MT model, a melt fraction in the range 15%–35% is the minimum required (Fig. 8D). Note that in Figure 8, the variable m is called the cementation exponent and controls the degree of interconnection of the partial melt, using the modified Archies’ Law of Glover et al. (2000). A value of m = 1 corresponds to well-connected melt, while m = 2 corresponds to melt in isolated inclusions. This range of melt fraction is consistent with the predictions from seismic studies that are described in Section 5.5.

In this region, a strong conductor (C3) is present at mid-crustal depths with resistivity values in the range 3–10 Ωm. In the east, close to Volcan Antofalla (elevation 6410 m), the lowest resistivity is observed at a depth of 30 km. The layer deepens to the west, and the lowest resistivity is observed at a depth of 50 km beneath the volcanic arc at Volcán Lastarria. As was the case for the mid-crustal structure beneath Cerro Galan, petrological data can be combined with the 3-D resistivity model to obtain estimates of the melt fraction to estimate the amount of melt that is present. This layer has been sampled by volcanic rocks sampled at both Volcán Lastarria (Stechern et al., 2017) and Volcán Antofalla (Richards et al., 2006). Since there are differences in the depth and composition of lavas at each location, separate calculations are presented.

At Volcán Antofalla, the earliest eruptions at 11 Ma were rhyolitic ignimbrites. From 9.1–1.6 Ma, the eruptions were primarily basaltic-andesite and andesitic-dacite in composition. Thus, to compute the present-day melt resistivity, an andesitic composition was assumed. A depth of 30 km below sea level corresponds to 35 km below the surface and an approximate pressure of 875 MPa. There are limited reports on the water content of the lavas; therefore, a range of 2%–3% was assumed, as reported elsewhere for lavas in the Central Andean subduction zone. Magma storage temperatures in the range 1000–1100 °C were reported from the Antofalla Volcanic Complex by Richards et al. (2006). Thus, at 1050 °C with a water content of 2.5%, a melt resistivity of 2 Ωm is predicted (Fig. 8B). With water content in the range 2%– 3% and pressure in the range 750–1000 MPa, an andesitic melt resistivity is predicted in the range 3.1–1.3 Ωm (Guo et al., 2016, 2017). To account for the 3–10 Ωm observed in the MT model, a melt fraction greater than 20% is needed (Fig. 8E).

At Volcán Lastarria, the deepest lava samples were identified as Group C by Stechern et al. (2017) and inferred to have erupted from pressures of 600–900 MPa and temperatures of 840−970 °C. Taking the mid-point of each range results in a pressure of 750 MPa and a temperature of 905 °C. Assuming an andesitic composition, a range of water contents 2%–3% predicts an andesitic melt resistivity of 4.6–10 Ωm (Guo et al., 2017). These relatively high melt resistivities would require very high melt fractions (50%–100%) and are inconsistent with the results beneath Volcán Antofalla.

It is important to recognize that geophysical studies image present-day structure, while petrology images the magmatic system at the time when lavas were erupted. The region of the Southern Puna on both sides of the Salar de Antofalla contains numerous mafic lavas of late Neogene age (Risse et al., 2013; Murray et al., 2015 and references therein). These lavas are believed to have been derived from mantle melting associated with lithospheric delamination and subsequent asthenospheric upwelling (Kay et al., 1994). It is possible that basaltic melt remains in the crust and is the cause of the low resistivity observed in the crust in this study. Data from laboratory studies are used to evaluate this hypothesis in the following sections.

The detailed study of Risse et al. (2013) showed that the mafic lavas erupted in the Southern Puna had undergone differentiation at two distinct depths within the crust. The first was at a depth of 60 km near the base of the crust with a temperature of 1375 °C and a pressure of 2 GPa. The second occurred in the mid-crust at depth of 20–35 km below sea level, which corresponds to the low-resistivity features C3 and C4 in Figures 6 and 7. Risse et al. (2013) showed that equilibration temperatures were in the range 1320–1140 °C with water contents low and less than 0.5 wt%. Laboratory studies are required to determine the electrical resistivity of these mafic melts. The study of Ni et al. (2011) was used, and the predictions are shown in Figure 8C where melts with water contents of 0.1, 0.3, and 0.5 wt% are predicted to have a resistivity of 0.3–0.7 Ωm at 1300 °C. It should be noted that 1300 °C is significantly below the temperature range used in the study of Ni et al. (2011). However, limited laboratory studies of mafic melt conductivity are available; thus, extrapolation of the results of Ni et al. (2011) is the most reliable option.

The bulk resistivity is computed with a similar approach to that used in sections 5.2 and 5.3, as illustrated in Figure 8F. To account for a bulk resistivity in the range 3–10 Ωm in features C4 and C5, at least 5%–15% interconnected basaltic melt is required. This value is certainly plausible, suggesting that basaltic melts may be responsible for the mid-crustal low resistivity beneath the Puna. Note the amount of basaltic melt required to produce a resistivity in the range 3–10 Ωm is significantly lower than corresponding amounts of rhyolitic and andesitic melt calculated in sections 5.2 and 5.3.

Directly beneath Volcán Lastarria, the upper crust in Figures 6 and 7 is characterized by relatively high resistivities, which agree with the 3-D inversions described by Díaz et al. (2015). The 3-D inversion of Díaz et al. (2015) identified a shallow conductor at a depth 3–5 km below sea level that was interpreted by Pritchard et al. (2018) as requiring a relatively high melt fraction if due to the presence of dacitic or andesitic melt.

Both prior studies by Díaz et al. (2015) and Budach et al. (2013) reported a deeper conductor to the east-southeast of Lastarria. In both studies, resolution of this feature was limited by the localized nature of the array of MT stations. This feature is imaged in the 3-D resistivity models presented in this paper and labeled C5. The combination of the three MT data sets gives better resolution than was possible in the individual inversions of Budach et al. (2013) and Díaz et al. (2015). However, the irregular grid of MT stations around Lastarria and just a profile east of Lastarria are far from ideal for resolving this feature. The feature begins at 5 km below sea level beneath the surface inflation center and then dips east where it merges with the mid-crustal conductor C3. The geometry of this feature in Figures 6 and 7 is similar to that reported by Budach et al. (2013). Feature C5 is also shown in a 3-D isosurface view in Figure 9A.

Models for uplift at Lazufre have computed inflation sources in the depth range 8–13 km below the surface (Pritchard and Simons, 2004; Perkins et al., 2016 and reference therein). These depths would place an inflation source close to the shallowest, more western part of C5. It is also significant that the upper limit of C5 is located close to the point where uplift occurs most rapidly, some distance south of Volcán Lastarria.

Geological and geodetic studies have shown that long-term uplift at Lazufre over the past 400,000 years has formed a topographic dome with horizontal spatial scale greater than 50 km (Perkins et al., 2016). It is proposed that the low-resistivity feature C5 represents the pathway by which melt has been supplied to the volcanic arc at Lazufre. The feature C5 in Figures 6 and 7 has a very low resistivity, which would require a high melt fraction (50%–100%). At the surface, it is characterized by fumarolic activity and numerous sulfur flows. It is also possible that the location of C5, on the edge of the array of MT stations, causes the inversion to assign this feature a lower resistivity than the actual value.

It is possible that aqueous fluids or hydrothermal alteration contribute to the low resistivity—i.e., the low resistivity is not solely due to partial melt, and feature C5 is due to the presence of saline fluids. These fluids can exhibit high conductivities with typical values reported in the range 5–30 S/m and possibly up to 100 S/m (Manning, 2018; Sakuma and Ichiki, 2016). Note that these conductivity values can be two orders of magnitude more conductive than partial melts (Comeau et al., 2020). The experimentally derived model of Sinmyo and Keppler (2017) and of Guo and Keppler (2019) can be used to estimate conductivities of H₂O-NaCl fluids for a range of depths and temperatures, for a range of possible salinities measured in wt% NaCl. Assuming saline fluids with conductivities of 10–30 S/m, corresponding to a salinity of 3.3 wt% to 5.6 wt%, values that are plausible based on geochemical data (Sinmyo and Keppler, 2017), an anomaly with a resistivity of 10 Ωm can be explained with a volume fraction less than 1%, whereas an anomaly of 1–3 Ω can be explained with a volume fraction of less than 5%. The presence of large amounts of saline fluids in the crust is significant as a number of studies have proposed that some of the fluids that formed the salars in this region of the Andes have a magmatic origin (Taussi et al., 2019; Godfrey and Alvarez-Amado, 2020).

The electrical resistivity model in Figures 6 and 7 shows a low-resistivity zone at mid- to lower-crustal depths in two main regions: (1) from Volcán Lastarria to Volcán Antofalla and (2) adjacent to Cerro Galan. This low resistivity can be interpreted as being due to the presence of partial melt, which will also lower the seismic velocity. Thus, a comparison of the 3-D resistivity model with a coincident seismic-velocity model allows both models to be critically evaluated.

Figure 10 shows coincident vertical sections from the 3-D resistivity model and two seismic models. The first is the P-wave tomography model of Chen et al. (2020). The second is the surface-wave tomography model of Ward et al. (2017). It should be noted that a shear-wave surface-wave tomography model for the southern Puna was initially published by Delph et al. (2017) and was expanded by Ward et al. (2017) to cover a greater area of the Central Andes. The figures in this paper use the surface-wave model of Ward et al. (2017). Figure 11 shows a comparison of horizontal sections of the tomography model of Ward et al. (2017).

In making a comparison between a 3-D resistivity model and a seismic-velocity model, it is important to recognize the different physics and data coverage of the various geophysical methods and surveys. MT impedance data are primarily sensitive to structures beneath a measurement station. Induction vectors are sensitive to structure laterally offset from the profile. Given the linear array used in the southern Puna for this study, resolution will be best close to the MT profile and will decrease to the north and south where there were no MT stations. The P-wave tomography model is also sensitive to structure beneath the array of seismic stations, which in the case of Chen et al. (2020), was oriented east-west. Away from the array of seismic stations, the sensitivity decreases. Surface-wave studies in the Central Andes, such as that of Ward et al. (2017) and Delph et al. (2017) used an array with a very large areal extent. This type of array gives generally good model resolution in the horizontal direction over the entire area shown in Figure 11. Vertical resolution is obtained from the frequency dependence of surface-wave velocity.

In Figures 10 and 11, the surface-wave model shows two prominent low-velocity zones with velocity <2.9 km/s and a decrease of up to 25% compared to the surroundings. These zones are labeled the Lazufre anomaly (LA) and the Cerro Galan anomaly (CGA). The features were observed at mid-crustal depths (15–30 km below sea level), but Delph et al. (2017) hypothesized that the two merged into a single feature in the lower crust at depths below 40 km. The low-velocity feature beneath Cerro Galan has an estimated volume of ~22,000 km³, compared to a cumulative erupted volume of ~1300 km³. These features were interpreted as being due to the presence of partial melt, with an estimated melt percentage of 22%. This melt fraction is similar to estimates derived from the MT data in this study. The low-velocity feature below Lazufre is smaller, whereas the low-velocity feature below Cerro Galan extends both north and south of the caldera location. Delph et al. (2017) proposed a model in which weak crust and crustal magma reservoirs are located at the edges of the Puna (the LA and CGA), and, after delamination of the lithosphere and upwelling of the mantle (which could have occurred in a short period of time on the order of 1–10 m.y.), basaltic magmas were intruded in the central region, which may correspond with the locations of feature C1, C3, and C4 in the electrical resistivity model.

A number of P-wave tomography studies have also shown evidence for low-velocity zones in the crust of the southern Puna. A low-velocity region was observed beneath Cerro Galan and the entire southern Puna by Bianchi et al. (2013) and was named Southern Puna Magma Body (SPMB). Chen et al. (2020) used local earthquake tomography to study crustal velocity and similarly reported a pronounced low velocity beneath Cerro Galan and referred to it as being part of the SPMB. The P-wave tomography model of Chen et al. (2020) is shown in Figures 10B and 10C. The low-velocity feature adjacent to Cerro Galan is coincident with features C2 and C4. The low-velocity feature beneath Lazufre is deeper than in the surface-wave tomography models and dips into the mantle wedge. The geometry of this low P-wave velocity feature is similar to that of a low-resistivity feature in the MT model (C6). It could be interpreted as a pathway of partial melt generated in the mantle wedge above the Nazca plate and rising through the mantle wedge and into the crust of the Southern Puna.

The low-resistivity zones detected in the mid and lower crust beneath the Puna require the presence of partial melt. The amount of partial melt is dependent on the composition, with 5%–10% mafic melt or 20%–40% andesitic or rhyolitic melt being required to account for the observed mid-crustal resistivity of 3–10 Ωm. The locations of zones with the lowest resistivity generally coincide with regions with the lowest S-wave velocities detected in seismic surveys. The consistency of the MT and seismic results gives strong evidence for melt distribution (1) close to the volcanic arc and (2) to the east where an elevated melt content likely reflects a melt influx to the crust that led to the eruption of Cerro Galan. The interpreted distribution of melt beneath the Puna is shown in Figures 9B and S11. The melt must originate in the underlying mantle wedge, although this region is not well imaged by MT data that do not include frequencies less than 0.001 Hz. Melt likely originates in this region due to an upwelling that is caused by downwelling of delaminated lithosphere. Melt may also originate close to the volcanic arc where fluids released from the subducting Nazca plate lower the melting point of the mantle wedge (C6).

Mid- and lower-crustal resistivity is significantly higher at 25°S than at 22°S on a profile at Volcán Uturuncu that crossed the Altiplano Puna Magma Body (APMB) (Comeau et al., 2015). At the broadest level, this implies a lower melt fraction across the southern Puna at 25°S than beneath the APMB. This observation is consistent with the overall tectonic framework, which proposes that changes in slab angle have migrated south (Kay and Coira, 2009). This interpretation suggests that at 22°S, the amount of melt in the crust has reached a peak with formation of the APMB. At 25°S, there is an ongoing influx of melts that has resulted in an increasing melt content and may have been associated with the eruption of Cerro Galan ca. 2 Ma (Ward et al., 2017).

The deep regional strike inferred from phase tensors and induction vectors is parallel to the Salar de Antofalla. The 3-D resistivity model is not extensive enough to see the exact structures causing this, but it is significant that this is not parallel to the volcanic arc.

The shallow magmatic system at Lazufre appears to be connected to zones of partial melt in the mid and lower crust to the east by the dipping feature C5 (Figs. 6, 7, and 9). The upper limit of C5 is within the depth range at which the inflation source is inferred to be located. Its spatial location is close to the center of the observed Lazufre surface uplift, which is offset from Volcán Lastarria. Lazufre provides another example of a magmatic system where melt transportation occurs in a direction that is non-vertical; it is becoming accepted that this is not uncommon (Lerner et al., 2020).

The resistivity structures at the two locations investigated by the PLUTONS project are significantly different within the crust, which thus implies significant differences in the distribution of partial melt within the Andean crust at these locations.

Volcán Uturuncu is underlain by the large conductor that was identified as the Altiplano Puna Magma Body (APMB). Uplift at this location has been determined to be short-lived and likely due to hydrothermal processes, rather than the formation of a diapir (Perkins et al., 2016). The MT models obtained at Volcán Uturuncu image both the deeper APMB and the shallow magmatic and/or hydrothermal system beneath the volcano but do not show a very strong conductor connecting the two.

At Lazufre, the grid of MT stations described in this study is more limited, mostly due to logistical issues than was the case at Volcán Uturuncu. However, the MT data collected to date show that there may be a connection between the shallow magmatic system and regions with a high melt content in the lower crust. This is an interesting result that could motivate a more detailed future study, one that could include a denser grid and broader coverage.

The 3-D electrical resistivity model presented in this paper provides new insights into the structure of the Southern Puna ~25°S in the Central Volcanic Zone of the Andes. The features observed in the 3-D electrical resistivity model were identified and interpreted. Together they suggest that the crustal melt distribution at the two locations studied during the PLUTONS project are significantly different. This result is consistent with the multi-disciplinary studies that have investigated the reasons for uplift, summarized in Pritchard et al. (2018). The main conclusions are as follows: (1) low-resistivity anomalies are identified in the mid-crust and can be realistically associated with the presence of the Southern Puna Magma Body (SPMB); (2) they are located near Lazufre and Cerro Galan; (3) they are consistent with the location of low shear-wave velocity anomalies; (4) below Cerro Galan, the low-resistivity anomaly can be explained with 22% partial melt; (5) below Lazufre, an eastward-dipping feature connects a shallow magmatic and/or hydrothermal system that requires fluids with a mid-lower crustal zone of partial melt; and (6) this feature is observed within the depth range of the inflation source below Lazufre, both of which are laterally offset from the volcanic cone.

Additional MT data are needed to fully image the resistivity structure of this important and dynamic region of the Central Andes, preferably with a regional-scale 3-D array (as allowed by the terrain), given the 3-D nature of the previously collected data.

Data collection was supported by the National Science Foundation (NSF) Continental Dynamics Program through award EAR-0908281 to Cornell University and Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant RGPIN-2020-04240 to Martyn Unsworth. The 3-D inversions used the ModEM algorithms, and we thank Gary Egbert, Anna Kelbert, and Naser Meqbel for allowing the use of their code. The inversions were performed on the clusters operated by Compute Canada. We thank Alan Jones and Gary McNeice for providing the STRIKE program. Many people are thanked for their assistance with the challenging task of fieldwork in the Puna; thanks to Victor Cuezzo, Dante Scatolon, Hermann Binder, and Pablo “Pichi” Maciel for their support. The original data and 3-D inversion models are available from the authors on request. This includes (1) the processed MT data in EDI format and (2) the 3-D inversion model in NetCDF format. We thank Suzanne Kay and an anonymous reviewer for comments that greatly improved this paper.