Magnetotelluric (MT) data were used to create a three-dimensional electrical resistivity model of the Altiplano-Puna magma body (APMB) in the area surrounding Volcán Uturuncu in southern Bolivia. This volcano is at the center of a zone of surface deformation with a diameter of 150 km and persistent inflation of ∼10 mm/yr. Low electrical resistivities (<3 Ωm) at a depth of 14 + 1/–3 km below sea level (16–20 km below surface) are interpreted as being due to the presence of andesite melts in the APMB, and require a minimum melt fraction of 15%. The upper crustal resistivity structure is characterized by finite-length, dike-shaped conductors, oriented approximately east-west near sea level. A combination of dacite partial melts and aqueous fluids is required to explain the observed low-resistivity values. Geodetic data do not require any deformation in these shallow regions. The geometry of the upper surface of the APMB beneath Volcán Uturuncu is consistent with that predicted by geodynamic models that suggest that the APMB bulges upward directly beneath Volcán Uturuncu, near the measured inflation center (∼3 km west of Volcán Uturuncu). Viscosity estimates from the MT-derived resistivity model gives a maximum value of 1016 Pa·s and is consistent with models that propose diapir-like ascent of magma above the APMB. Resistivity models are compared and quantitatively correlated to seismic velocity models, showing good agreement on the spatial extent and depth of the APMB. A forward modeling study shows that the small differences in the depth to the top of the APMB between the different geophysical methods could be explained by variations in the composition of the magma body.
The Andean subduction zone is a convergent plate boundary, where the oceanic Nazca plate subducts beneath the continental South American plate (Cahill and Isacks, 1992). The Central Andes region contains the Altiplano and Puna, which together compose one of the largest high plateaus on Earth (Isacks, 1988; Allmendinger et al., 1997). High plateaus typically form in continent-continent collision zones, such as the Tibetan Plateau or the East Anatolian–Iranian plateau (Royden et al., 2008; Unsworth, 2010), and rarely form above a subduction zone. The Altiplano-Puna plateau is the only modern example found in the 55,000 km of subduction zones present on the Earth today (Stern, 2002). The Altiplano-Puna has an average elevation of 4000 m and an average crustal thickness of 65 km that reaches a maximum of 75 km (Beck et al., 2014). Crustal thickening has occurred by both magmatic addition and thrusting (Allmendinger et al., 1997). This region is currently associated with subduction at a relatively steep angle of 30°, although flat-slab subduction has previously occurred (e.g., Cahill and Isacks, 1992).
Volcanism in this region is characterized by a well-developed volcanic arc and a number of backarc volcanic centers, including Volcán Uturuncu, located in southern Bolivia (22.265°S, 67.187°W). A steepening of the subducting slab caused a major ignimbrite flare-up during the past 10 m.y. that formed the Altiplano-Puna volcanic complex (APVC) (de Silva, 1989). Geophysical studies have shown that the Altiplano is underlain by (1) anomalously low seismic velocities (ANCORP Working Group, 2003; Chmielowski et al., 1999; Beck and Zandt, 2002; Zandt et al., 2003), (2) high seismic attenuation (Haberland et al., 2003; Schurr et al., 2003), (3) low electrical resistivities (e.g., Schwarz and Kruger, 1997; Brasse et al., 2002; Schilling et al., 2006), (4) anomalously high heat flow values (Henry and Pollack, 1988; Springer and Förster, 1998; Hamza et al., 2005), and (5) a negative Bouguer anomaly (e.g., del Potro et al., 2013). Together these anomalies give evidence for the presence of a major mid-crustal magma body, the Altiplano-Puna magma body (APMB), and make this region a key location for understanding silicic magmatism in continental crust. Furthermore, satellite-based geodetic data have detected a spatially broad (∼150 km diameter) and temporally continuous (decades long) deformation pattern coincident with the APMB and centered on Volcán Uturuncu, which includes a central ring of uplift (10–15 mm/year) and an outer ring of subsidence (2 mm/year) (Fialko and Pearse, 2012; Henderson and Pritchard, 2013). Geodetic models place the inflation source in the depth range 15–20 km below sea level (BSL) (Pritchard and Simons, 2004), located within the APMB, as determined from geophysical data (Zandt et al., 2003). The shape of the body causing the inflation is not well constrained, and several models have been proposed, which are all consistent with the available geodetic data. The available geophysical data near Volcán Uturuncu (seismic and gravity) are not able to resolve crustal features smaller than a few kilometers above the APMB. These features may be important because they may define the location of melt ascending from the APMB to the surface. Understanding magmatism, and the potential for future eruptions in this region, requires that these features are resolved. Therefore, additional geophysical data are needed.
The electrical resistivity of a rock is sensitive to the quantity, geometry, and composition of melt in a partially molten region. Magnetotelluric (MT) data can measure the electrical resistivity at depth, and these data can constrain subsurface magma distribution (Chave and Jones, 2012). An array of 180 broadband MT sites was collected from 2011 to 2013 to investigate the crustal magma distribution below the APVC; site locations are listed in the KML file (Supplemental File 11). An initial analysis of these MT data was presented in Comeau et al. (2015); the MT data were used to generate an ∼300-km-long regional two-dimensional (2-D) resistivity model extending across the APVC. The dominant feature of the model was a low-resistivity layer at a depth of 15–20 km BSL that was identified as the APMB. The upper surface of this layer bulged upward below Volcán Uturuncu. Comeau et al. (2015) also presented a 3-D resistivity model around Volcán Uturuncu using a subset of the MT data, which showed both a regional conductor at a depth of ∼15 km BSL (∼20 km below surface) and a shallow conductor at a depth near sea level (∼5 km below surface). Note that the 2-D inversion placed the APMB conductor at a greater depth than the 3-D inversion, likely due to the effect of static shifts (Comeau et al., 2015). The depth obtained from the 3-D MT inversion gives better agreement with the seismic data than the 2-D inversion, and is thus considered more reliable.
This paper presents a more complete analysis of this MT data set, including a 3-D inversion that uses additional MT sites and a quantitative comparison with seismic velocity models, and attempts to reconcile the fact that different geophysical methods image the APMB at slightly different depths. In addition, the resistivity values derived from the MT data inversion are used to estimate melt fractions in the magma bodies and to estimate their effective viscosities in order to better constrain geodynamic models.
3-D MT INVERSION AND MODEL RESOLUTION TESTS
3-D MT Inversion
In Comeau et al. (2015), the 3-D inversion of a set of 73 MT sites (grid S73) was described using the algorithm of Siripunvaraporn et al. (2005). This MT inversion algorithm finds a resistivity model that fits the MT data to within a specified statistical tolerance. Because many models can be found to fit a given data set, the resistivity model is also required to satisfy conditions such as being spatially smooth, i.e., it is regularized. The inversion algorithm is iterative and starts from a user-specified resistivity model. Additional regularization can be implemented through the use of a second (a priori) resistivity model. The computational cost of an inversion increases as the number of MT sites is increased, so the entire data set (180 sites) could not be inverted. With careful data selection it was possible to invert the 96 MT sites shown in Figure 1 (grid S96) with 18 frequencies from 115 to 0.00083 Hz.
The 3-D inversion of S96 also used the algorithm of Siripunvaraporn et al. (2005), and a two-stage approach was used. For the first stage, both the prior model and the starting model were a uniform half-space (10 Ωm) and the inversion algorithm searched for a spatially smooth resistivity model close to the prior model that was able to fit the measured data. After 8 iterations the root mean square (r.m.s.) misfit was reduced from 8.74 to 1.47, using a 10% error floor on all impedance components (diagonal and off diagonal). A second inversion run was then begun, using the final model from the first stage as the prior model. After 3 iterations the r.m.s. misfit was reduced to 0.95; this low misfit indicates an acceptable fit to the measured MT data. The preferred 3-D resistivity model was obtained from the joint inversion of all elements of the impedance tensor and is shown in Figure 2. The resistivity model is also available in digital form online. Inversions without the diagonal components of the impedance tensor (i.e., Zxx, Zyy) were also run, but it is preferred to include the full impedance tensor to properly resolve the small-scale 3-D structures (Tietze and Ritter, 2013; Kiyan et al., 2014; Newman et al., 2015).
Tipper data were measured in the field and were used in the 2-D inversion of Comeau et al. (2015) to locate the edges of the large, deep conductor interpreted to be the APMB. In the 3-D inversion, models from inversions with and without the tipper data gave similar results. This is because the inversion domain around Volcán Uturuncu did not extend as far as the edges of the APMB. The inversion without tipper was preferred because it allowed a larger model to be used.
Topographic effects influence the higher frequency MT data, and this problem can be addressed by including topography in the inversion. However, this increases the number of model parameters, and significantly slows the inversion. Topography was included in the 2-D model of Comeau et al. (2015), but was not included herein because the focus was on intermediate depth (∼5 km below surface) structure and the highest frequency was 115 Hz; therefore, topographical effects may not be significant. The surface was defined as 5000 m above sea level. Earthquake hypocenters from Jay et al. (2012) (Fig. 2) appear clustered near sea level at the edges of conductors. The data fit is illustrated in Figure 3, which shows the sounding curves at select MT sites. Figure 4A shows the model as a set of resistivity-depth profiles beneath each MT site.
Many combinations of inversion parameters were investigated and the main features of the resistivity models did not depend on any specific parameters choices, implying that the features are robust and are required by the MT data. The various inversion parameters tested are summarized in Appendix 1 (Table A1).
Resistivity Model Features
The resistivity model presented in Figure 2 is similar to that in Comeau et al. (2015). Three distinct layers are present: (1) a shallow heterogeneous surface conductor ∼1 km thick (C1); (2) a deep, spatially uniform conductor that begins at a depth of 10–15 km BSL, below which there is no resolution (C2); (3) an intermediate high-resistivity layer with a prominent set of east-west–oriented dike-shaped conductors located at approximately sea level (C3).
Shallow Features (<1 km below Surface)
The layer C1 is confined to the upper 1 km and has highly variable resistivity, with regions as low as 0.1 Ωm and as high as 10,000 Ωm. The low-resistivity regions are likely caused by a shallow aquifer and a layer of hydrothermal alteration within and below the regionally extensive ignimbrite near-surface layer. Similar features were observed beneath ignimbrites in the Taupo Volcanic Zone in New Zealand and it appears that the low-resistivity layer is composed of clay minerals formed by low-temperature hydrothermal alteration over a time scale of ∼1 m.y. (Bibby et al., 2005). The ignimbrites in the MT study area were erupted in four major pulses (Vilama, 8.5 Ma; Guacha, 5.5 Ma; Atana, 4 Ma; Pastos Grandes, 3 Ma), each eruption moving successively farther west by a total of ∼100 km (Salisbury et al., 2011). Given their age, it is reasonable to assume that there has been sufficient time for alteration to occur and produce the observed low resistivity.
Intermediate Depth Features (–3–5 km BSL)
At intermediate depths (3 km above sea level to 5 km BSL; 2–10 km below surface), series of east-west–striking, dike-shaped conductors are observed beneath Volcán Uturuncu (C3), the Volcán Uturuncu high-conductivity zone. These conductive features are located between C2 and the surface and are significant because they may represent present-day zones of magma and aqueous fluids or hydrothermal alteration caused by pathways of past fluid motion. The conductors below Volcán Uturuncu may also be interpreted as shallow magma chambers, as suggested by other studies based on geochemical evidence (Muir et al., 2014), seismicity (Jay et al., 2012), and a low-velocity zone identified by a local ambient noise tomography (ANT) network (Jay et al., 2012). The 3-D resistivity model shows a resistivity of ∼3 Ωm for these features. The approximate volume of each conductor is ∼500 km3. These dike-shaped conductors extend vertically downward toward the APMB and seismicity is confined to their edges. These conductive dikes appear to be feeding the volcanic system from below, rising from the upward bulge in the APMB below the inflation center near Volcán Uturuncu. Low resistivity in a volcanic environment can be caused by aqueous fluids, partial melt, or hydrothermal alteration, so interpretation must use additional information to distinguish between these possibilities (Evans, 2012). Similar long and narrow conductors were detected by the Earthscope (http://www.earthscope.org/) MT data in the northwestern U.S. (Meqbel et al., 2014) and were interpreted to be the response of the isotropic 3-D inversion to large-scale electrical anisotropy. This is unlikely to be the case for this data set because the shallow conductors below Volcán Uturuncu are much larger than the MT site spacing, and are relatively well resolved.
Deeper Features (>10 km BSL)
The major resistivity feature at depth is the prominent conductor C2. Resistivity depth profiles show this feature clearly (Fig. 4A). Note that there is a high degree of variability in the upper 10 km of the resistivity model, where different MT sites are detecting intermediate conductors of different geometries near sea level, such as C3. However, at depths >10 km BSL, a similar resistivity-depth profile is observed at many sites. The model resistivity decreases smoothly to a minimum resistivity of 0.3 Ωm at a depth of ∼20 km BSL, and then smoothly increases. As with all inversions, the nonuniqueness of the model should be considered when interpreting resistivity values and model features.
Defining the depth of the top of C2 must be undertaken with care, because the MT inversion approach used in this study seeks the spatially smoothest resistivity model consistent with the MT data. Thus even if the true resistivity model was a layer with sharp edges, it would appear in the MT inversion as a smooth feature. This phenomenon is illustrated in Figures 4B and 4C using synthetic MT data. Two 1-D synthetic resistivity models are shown representing a low-resistivity half-space and a low-resistivity layer. Synthetic MT data were generated for these models and 5% Gaussian noise was added to the impedance data. The MT data were then inverted using the same approach as for the field MT data. Figure 4 shows that the inflection point of the resistivity-depth curve corresponds to the true depth of the top of the low-resistivity region, for both the half-space case and the layer case, and that the bottom of the layer is not well resolved. Thus when interpreting the resistivity model of the APMB (Fig. 2), interpretation should focus on the inflection point that occurs at a depth of 13–14 km BSL (18–19 km below surface). Defining the depth from a specific resistivity contour is subjective. For example, if the top was defined as the 3 Ωm contour in the 3-D model, then the conductor would be correctly located at a depth of ∼14 km BSL, but other choices would give incorrect results.
The exploration depth of MT is controlled by the electromagnetic skin depth phenomena, and is dependent on both the period of the signal and the ground resistivity. MT data can give a reliable estimate of the top of a conductor, but it is not always possible to record long enough periods (low frequencies) to detect the base of a conductor. This is the case with C2, and the hachured pattern on the resistivity model (Fig. 2) indicates where there is minimal resolution (calculated as 1 skin depth at a period of 1000 s). Synthetic inversions also confirm that model resolution is lost within C2 (Comeau et al., 2015).
Based on the horizontal extent, the low-resistivity values observed, and the laterally uniform resistivity, C2 can be identified as the APMB. The low resistivity can be interpreted as a zone of partial melt, and this interpretation is justified and quantified later in this study. The thickness of the APMB is not constrained because the very low frequencies needed to pass through the APMB were not measured. However, the MT data show that the APMB has a very high conductance that exceeds 20,000 S (conductance is the product of a layer’s conductivity and thickness). If the APMB is assumed to have a minimum resistivity of 0.3–1 Ωm this gives a minimum thickness of 6–20 km for the APMB from the MT data. In stable continental regions the average crustal conductance is generally no more than 1000 S, although a conductance of as much as 20,000 S has been reported in Tibet (Wei et al., 2001). Previous MT studies of the Altiplano found a conductance >20,000 S (Schwarz and Kruger, 1997; Brasse et al., 2002).
Alternatively, it is possible to estimate the thickness of the APMB from its lateral extent, as detected with MT data, and its volume (estimated from the plutonic to volcanic ratio; intrusive:extrusive ratio). A plutonic to volcanic ratio of 2:1–5:1 is commonly used in the Andes (Francis and Hawkesworth, 1994; de Silva and Gosnold, 2007), although higher values have been proposed and include values of 10:1 (Kay et al., 2010) and 30:1 (Ward et al., 2014). The APVC erupted ignimbrite volume is estimated to be 15,000–30,000 km3 (de Silva et al., 2006; Schilling et al., 2006). The area of the APMB estimated from the regional resistivity models (Comeau et al., 2015; Comeau, 2015) is ∼14,000 km2. This gives an estimated thickness of 2–11 km for the APMB, as shown in Comeau (2015). Earlier seismic receiver function studies imaged the APMB as a thin sill with a thickness of ∼1 km (Chmielowski et al., 1999; Zandt et al., 2003; Leidig and Zandt, 2003). The combination of ambient-noise tomography with receiver function data has greatly improved the reliability of the seismic results, as shown by Ward et al. (2014), who imaged the APMB as a low-velocity zone ∼10–15 km thick.
Resolution Tests of the 3-D MT Inversion Model
The 3-D resistivity model in Figure 2 was investigated using a set of synthetic MT inversions. These allow the model resolution and sensitivity to the measured MT data to be investigated. The synthetic resistivity model (model A in Fig. 5A) was kept simple and included the main resistivity features seen in the data inversion model in Figure 2. It included three main layers, as described in the following.
A 1-km-thick surface layer with resistivity of 3 Ωm represents the hydrothermal alteration within and beneath the ignimbrite layer.
An upper crustal resistive layer (100 Ωm) contains 3 parallel dike-like conductors, each with a resistivity of 3 Ωm, and with tops at a depth of –3 km BSL. These dike-like conductors are 15 km long, 4.5 km wide, and 6 km thick.
The model terminated with a low-resistivity half-space (3 Ωm) at a depth of 15 km BSL, representing C2.
Synthetic MT data were computed in the period band 0.01–1500 s, similarly to the real MT data inversion, and 5% Gaussian noise was added. The actual MT site distribution was used in the synthetic inversions. The data were inverted with the algorithm of Siripunvaraporn et al. (2005). For model A, the r.m.s. misfit was reduced from 6.38 to 1.44, after 5 iterations. The synthetic forward model and the inversion results are compared in Figure 5. The true model and the synthetic inversion model are shown as depth-resistivity profiles in Figure 6. It can be seen that the depth of the deep conductor (C2) is recovered reliably, and the intermediate dikes (C3) are also imaged in the expected positions with the correct resistivities. The geometry of the dikes was not fully recovered, owing to the irregular spatial distribution of MT sites.
Model B was created by extending the dikes to infinity in an east-west direction, in order to investigate whether the MT data are sensitive to the east-west extent (Fig. 5C). For model B, the synthetic inversion reduced the r.m.s. misfit from 6.08 to 0.91, after 5 iterations. Model features C2 and C3 are imaged in the expected positions (Fig. 5D). However, owing to the irregular distribution of MT sites the structure is not as uniform as in the true model. Despite this, it is clear that the 3-D inversion can distinguish between dikes with a finite and infinite extent. This gives confidence that the features in the 3-D inversion model in Figure 2 can properly define the limits of the dikes.
Sensitivity of Model Features to Measured MT Data
An additional test was conducted to determine if the features present in the 3-D resistivity model were required by the MT data. This involved editing the resistivity model so as to remove the intermediate-depth dikes (C3), as shown in Figure 7. This analysis used a subset of the MT sites, and used the S73 3-D inversion model of Comeau et al. (2015), because this model focused on the central area of interest. The resistivity of the model in the depth range of 2–12 km below surface was set to a moderately resistive value (30 Ωm). The forward response of the edited model in Figure 7B was computed and compared to the response of the inversion model and the measured MT data. It was found that this increased the overall r.m.s. misfit from 0.95 to 2.27. The response of the forward modeled MT data as apparent resistivity and phase curves is shown in Figure 8. The edited model does not fit the measured MT data at periods in the range 3–300 s, the approximate period range of the intermediate features that were removed. Furthermore, it is clear that the yx impedance component is poorly fit, compared to the xy component. This is as expected because it is the yx impedance data, derived from the east-west electric fields, and analogous to the transverse electric mode in 2-D MT, that is sensitive to the east-west–oriented conductors. This implies that the dikes are required in the resistivity model, and are required to be oriented approximately east-west. Overall, by removing the conductive dikes the apparent resistivity would be expected to increase, and the phase to decrease. This is seen in the forward model responses of Figure 8, which shows that the MT responses with and without the dikes are quite different, and can be readily distinguished with the measured MT data.
The next step in the analysis was to start an inversion from the edited resistivity model, to determine if the inversion algorithm would reintroduce the dikes in order to fit the MT data, or whether it would be possible to fit the MT data with a different model. The inversion was restarted with the prior model equal to the edited model in Figure 7B, and the starting model a half-space of 10 Ωm. Inversion parameters were the same as in the original 3-D inversion. The resulting resistivity model in Figure 7C shows resistivity features similar to those of the original model in Figure 7A. Figures 7D and 7E show horizontal slices of the original and rerun model. This test shows that feature C3 of the inversion model is required by the MT data. The original inversion model r.m.s. misfit was 0.95, whereas the inversion model of the edited model has an r.m.s. misfit of 1.11, and the main features are returned. When the starting model was changed to the edited model the inversion result also returned the intermediate conductors, as expected. The diagonal components are noisy for periods <10 s and are therefore difficult to interpret. However, examining their fit (Figs. 8D, 8E) can provide physical insight. Once the model was edited, the magnitude of the diagonal components at MT sites located at the edges of C3 (i.e., site 4) were reduced by a factor of 100. This was as expected because 3-D features were removed and the resulting 2-D model has small diagonal components. The response of the inversion model shows that the magnitudes of the diagonal components increased as the 3-D dikes were replaced. For MT sites located away from the 3-D dikes (i.e., sites 1 and 3) little change was observed in the diagonal components when the model was edited, or the inversion restarted because the model is relatively 2-D in all cases. For MT sites located directly above the 3-D dike (i.e., site 2), a significant decrease in the magnitude of the diagonal components was observed for the yy component but not the xx component. When the model was reinverted the diagonal components are similar to the original inversion because the 3-D dike was replaced.
INFERRING MELT PROPERTIES FROM MT RESISTIVITY
Predicted Values of Melt Resistivity
The resistivity values derived from MT exploration can be used to determine the composition of a magma body. However, this is a nonunique process because the bulk resistivity depends on several parameters, including the melt composition, melt quantity, and melt geometry (see review by Pommier, 2014). Melt resistivity is dominated by both the water content and the concentration of sodium ions (Pommier and Le Trong, 2011; Nesbitt, 1993). The SIGMELTS relation of Pommier and Le Trong (2011) allows the electrical conductivity of silicate melts to be computed as a function of composition, temperature, and pressure. A geochemical analysis of Volcán Uturuncu lavas provides compositional values for rocks in the study area and allows the melt resistivity to be predicted (Muir et al., 2014; Sparks et al., 2008), as shown in Figure 9.
Dacite lava samples with silica contents of ∼65% were inferred to be derived from a shallow magma chamber (Muir et al., 2014; Sparks et al., 2008), corresponding to feature C3 in Figure 2. This shallow chamber is at a depth of 5 km below the surface (approximately sea level), corresponding to a pressure of ∼100 MPa (Muir et al., 2014). The dacite samples had water contents of 2.5–4 wt%, and a sodium content of 2–3 wt% (Muir et al., 2014). Thermobarometry by Sparks et al. (2008) predicts an average temperature of 838 °C for these dacite samples. Based on these values a pure melt resistivity of 5–15 Ωm is predicted by SIGMELTS (Fig. 9A). These melt resistivity values match the laboratory studies of Laumonier et al. (2014), who measured a melt resistivity of 2.5–6.5 Ωm at 870 °C and 130 MPa on a dacite sample. These resistivity values are relatively high for arc magmas, and likely reflect the relatively low sodium content of the Volcán Uturuncu lavas, compared to other systems (e.g., 3.5–5 wt% for Mount St. Helens, Pinatubo, Unzen; Muir et al., 2014, and references therein).
The analysis was repeated for andesite samples with a silica content of ∼60%, which were inferred to have been erupted from the deeper magma body beneath Volcán Uturuncu, corresponding to the APMB (C2 in Fig. 2). This body is at a depth of ∼15 km below surface, with a pressure of ∼375 MPa (Sparks et al., 2008). Based on melt inclusion data, the andesite magmas had a temperature of 980 °C, a sodium content of 2–3 wt%, and were water-rich, possibly with water contents in excess of 7 wt% (Sparks et al., 2008). These values predict a melt resistivity of 0.1–2 Ωm (Fig. 9B). These melt resistivity values match the laboratory studies of Laumonier et al. (2014), who predicted a melt resistivity of 0.2–2 Ωm at 870 °C and 1000 MPa. The combination of a higher temperature and an elevated water content gives the andesite melt a much lower resistivity than the dacite melt.
Melt Fraction Estimates
A suitable equation to combine the resistivity of the two phases (melt and rock matrix) is required. Archie’s Law ignores electrical conduction through the rock matrix, and this assumption is not valid at the high temperatures of the APMB where conduction through the rock matrix can be significant. Thus in this study, the modified Archie’s Law (MAL) of Glover et al. (2000) was used. The Hashin-Shtrikman upper and lower bounds were also used (Hashin and Shtrikman, 1963) because they give reliable bounds on the resistivity. Deriving the melt fraction from the observed resistivity represents an inverse problem, and using the minimum values of melt resistivity will result in the determination of the minimum melt fraction required to explain the resistivity of the 3-D model.
For the dacite melt a resistivity of 5 Ωm was determined and used for further calculations (Fig. 9A). The resistivity of the rock matrix was computed as 780 Ωm using the relation of Hashim et al. (2013) at a temperature of 838 °C. The bulk resistivity was then computed for a range of melt fractions using MAL with cementation exponents in the range m = 1–2. The Hashin-Shtrikman bounds were also used because crustal melts are generally interconnected (ten Grotenhuis et al., 2005; Unsworth and Rondenay, 2012), making the lower bound value the most realistic. Figure 9C shows the melt fraction required to explain the observed resistivity values of the shallow dacite body. The dike-like conductor beneath Volcán Uturuncu (C3 in Fig. 2) has an observed resistivity of 3–7 Ωm in the inversion resistivity models, and dacite melt cannot explain the observed resistivity without having an unrealistically high melt percentage (>80%) (Fig. 9C). Therefore, another conducting phase is needed to explain the observed resistivity values. Saline aqueous fluids, in addition to magma or alone, could explain the observed resistivity values. Saline fluids would be more conductive, helping to explain the observed low-resistivity value. The fluids could be exsolved from magma crystallizing in C3, or perhaps the APMB. Aranovich et al. (2013) suggested that such fluids can move vertically over significant distances and accumulate in the shallow crust.
The analysis was repeated for andesite lava samples, which were inferred to have been erupted from the deeper magma body beneath Volcán Uturuncu, corresponding to the APMB. The bulk resistivity was computed for andesite melt having a resistivity of 0.1 Ωm, as described here. The resistivity of the rock matrix was 325 Ωm for the deeper body at a temperature of 980 °C. Figure 9D shows the melt fraction required to explain the observed resistivity values of the deeper andesite body. The deeper conductor beneath Volcán Uturuncu (C2) has an observed resistivity of <1 Ωm in the inversion models. The Hashin-Shtrikman lower bound predicts that 15%–45% andesitic melt is required to explain resistivity values of 0.3–1 Ωm (Fig. 9D). Therefore, a minimum melt fraction of 15% is required if andesite melt is the cause of this conductive anomaly. These estimates are consistent with other studies that have suggested a melt fraction of 14%–27% for the APMB (Schilling et al., 2006). If higher values of pure melt resistivity are used then a higher melt fraction is required. Therefore these values represent the minimum melt fractions required by the MT data.
COMPARISON OF ELECTRICAL RESISTIVITY AND SEISMIC VELOCITY MODELS
It is difficult to interpret physical rock properties with a single geophysical method and multiple methods are commonly used to reduce the uncertainty. Both electrical resistivity and seismic velocity are sensitive to the composition of a rock, especially the fluid content (Bedrosian et al., 2007). Although there is no universal relationship between electrical resistivity and seismic velocity, they can be linked in certain lithologies through the amount of fluid present (Bedrosian et al., 2007). As the amount of fluid increases, both electrical resistivity and velocity generally decrease. This trend has been observed for rocks that contain both aqueous fluids and partial melts. However, some compositional changes will influence only the electrical resistivity (Unsworth and Rondenay, 2012). For example, increasing the sodium content of the melt will decrease the melt resistivity but will not significantly change the velocity. A joint inversion of velocity and resistivity data is possible, but can be of limited use if the two parameters are not controlled by the same physical process (Roux et al., 2011; Moorkamp et al., 2007). In this study a quantitative correlation is used to investigate the velocity and resistivity models of the APMB below Volcán Uturuncu in southern Bolivia.
Comparison and Correlation with Velocity Model of Ward et al. (2014)
A number of seismic studies have been undertaken to investigate the structure of the APMB in southern Bolivia and adjacent regions. Chmielowski et al. (1999) and Zandt et al. (2003) used receiver functions (RF), while Ward et al. (2013) used ambient noise tomography (ANT). Ward et al. (2014) combined ambient noise tomography and receiver function data (ANT + RF) to give the most complete image of the crustal velocity structure, and the best spatial coverage in the vicinity of Volcán Uturuncu. This study was used for comparison with the MT resistivity models in the following.
Figure 10 shows a comparison of the S-wave velocity model of Ward et al. (2014) and the 3-D electrical resistivity model interpolated onto a common spatial grid (that of the velocity model), at lat 22.3°S. The common grid is spatially coarse, with a horizontal cell size of 11 km, in comparison to the original resistivity model (Fig. 2), which has a horizontal cell size of 0.75 km. The resistivity model and seismic velocity model clearly show a correlation. The low-resistivity region below a depth of ∼10 km BSL is coincident with a broad low-velocity region, which is interpreted as the APMB (red in Figs. 10A, 10B). The higher velocity lid above the APMB is coincident with the higher resistivity region above ∼10 km depth BSL (blue in Figs. 10A, 10B). Because of the resolution differences between the two models, the resistivity model shows much finer scale spatial structures above the APMB that are not resolved in the seismic velocity model.
A quantitative approach to correlating two independent geophysical data sets was developed by Bedrosian et al. (2004). By correlating both parameters, distinct domains can be defined with certain ranges of resistivity and velocity. In order to do this type of quantitative correlation the finer resolution resistivity model needs to be interpolated so that it shares a common set of grid points with the velocity model. Next the poorly resolved areas of the models must be omitted to avoid inaccurate correlation. In the analysis shown here, the edges of the 3-D resistivity model where it is not constrained by MT sites and the region below the large conductor where there is minimal resolution have been removed to ensure proper correlation. The correlation results are then presented as a two-parameter histogram that shows the number of occurrences of each combination of resistivity and seismic velocity values (Fig. 10C). Distinct domains can be identified and mapped back onto the spatial model to better understand which zones show significant correlation (Fig. 10D).
The correlation results show a general trend of increasing resistivity with increasing velocity (Fig. 10C). Zone 1 can be spatially identified with the APMB and with the near-surface layer, with a correlation of low velocity (<2.75 km/s) and low resistivity (<20 Ωm). Zone 2 is identified as the lid above the APMB, and shows a correlation of high velocities (>2.75 km/s) and high resistivities (>10 Ωm). Zone 3 is an anomalous zone that does not follow the expected trend. It shows a correlation of low resistivity (<10 Ωm) with high velocity (>2.75 km/s). This zone is caused by the small-scale low-resistivity features above the APMB seen in the resistivity model but not resolved in the velocity model. Therefore the explanation of this anomalous zone may be due to the resolution difference between the methods, as discussed herein. Mapping the zones spatially back onto the model confirms the identification of the zones (Fig. 10D).
Reconciling Seismic and MT Estimates of the Depth to the APMB
An objective comparison of the models obtained from different types of geophysical data requires that the characteristics of each exploration method are understood. It was shown here (Fig. 4) that the regularized MT inversion will image a sharp step in resistivity as a gradational transition over several kilometers. It was shown that the inflection point (location of maximum vertical resistivity gradient) gives a good approximation to the depth of the top of the conductor. This approach gives a more robust method to locating the depth of a conductor than selecting an arbitrary resistivity contour value. In the 3-D MT inversion shown in Figure 11, the inflection point occurs at a depth of 13–14 km BSL, and corresponds to a resistivity value of 3 Ωm. This depth can be considered as the average depth to the top of the APMB over the region sampled in the 3-D MT inversion. Note the depth of the APMB as defined with MT is not constant and varies by several kilometers across the extent of the APMB. This depth can be compared to depth estimates from seismology. The early receiver function results of Zandt et al. (2003) gave a depth of 13–15 km BSL for the top of the APMB, as shown in Figure 11C, and matches the MT results quite well (Fig. 11A). However, the joint ANT + RF results (Fig. 11B) have an improved resolution and show a low-velocity zone in the depth range 6–20 km BSL beneath station LCOL (48 km west of Volcán Uturuncu), and 10–23 km BSL beneath station PTLP (48 km northeast of Volcán Uturuncu), using a velocity of 2.9 km/s to define the top of the APMB (Ward et al., 2014). The lowest velocity occurs at a depth of 13–15 km BSL (Ward et al., 2014). Therefore, if the top of the APMB is defined by the top of the 2.9 km/s velocity contour, then there is a discrepancy of 3–8 km between the two methods. However, if instead the APMB is located where the slowest velocities are recorded (i.e., <2.5 km/s), then all 3 methods agree on a depth of ∼14 km BSL to the top of the APMB and there is no discrepancy.
One issue with the ANT results is that this method uses a broad station spacing and is not sensitive to small-scale vertical structures. It only images a single smooth and continuous low-velocity zone (Ward et al., 2014) that only shows the long-wavelength structure of the APMB, and likely smears the small-scale features seen in the MT models in the horizontal direction. This could explain why the ANT + RF model gives a shallower depth to the low-velocity zone of the APMB, compared to other RF seismic models and the MT models. A synthetic study of how the ANT + RF joint inversion will image and smooth a sharp boundary is required to better explain these differences.
INTEGRATED PETROLOGICAL GEOPHYSICAL CONSTRAINTS ON THE STRUCTURE OF THE APMB
It was shown herein that both the 3-D MT resistivity model and the seismic velocity model detect a low-resistivity–low-velocity layer at a comparable depth that can be identified as the APMB. The agreement in depth for the two anomalies is relatively good, and the differences may be within the data uncertainty, or due to resolution differences between the two methods. However, it is possible that the properties of the APMB vary with depth and the two geophysical methods are sensing different layers of the APMB. If the upper part of the APMB was characterized by a low-velocity zone, but not a low-resistivity zone, then a difference in depths would be observed. This could be achieved by varying the composition of the upper portion of the APMB in such a way that the resistivity was affected but the velocity was not. In the following, a model of the APMB is developed and a range of properties investigated to determine if a realistic layered model could explain the differences between the ANT and MT models of the APMB.
Seismic Velocity Structure of a Layered Magma Body
The shear velocity is primarily sensitive to the amount of fluid or melt, so if the melt fraction is constant, the velocity will be constant throughout the magma body. The shear wave velocity as a function of melt fraction can be computed using the relations of Dvorkin (2008) and Takei (2000) (for details, see Appendix 2). A layer with 25% melt fraction results in a velocity of ∼1.5 km/s from the Dvorkin relation, which matches the ANT + RF results of Ward et al. (2014). A 15% melt fraction predicts a shear velocity of ∼2.2 km/s, still within the values for the low-velocity zone as defined by Ward et al. (2014).
Electrical Resistivity Structure of a Layered Magma Body
A model is developed to allow the resistivity-depth profiles to be computed for a layered magma body. The resistivity of the melt was calculated from the SIGMELTS model of Pommier and Le Trong (2011). Melt resistivity decreases as the concentration of water and/or cations increases. Sodium ions are the most effective charge carrier owing to their high mobility (Pommier, 2014). In contrast, the melt resistivity increases with silica content owing to an increase in the melt viscosity. The resistivity of the solid grains is controlled by ionic diffusion through the solid and is thermally activated with an exponential temperature dependence (Nover, 2005; Hashim et al., 2013). To calculate the resistivity of the partial melt, the resistivities of the two phases, the melt fraction, and the melt geometry are required. This calculation uses a two-phase mixing law because at the high temperatures found in a magma body, some electrical conduction will occur through the rock matrix. The MAL of Glover et al. (2000) was used for this analysis.
In the following, the effect of varying silica content, water content, melt connectivity, melt fraction, and sodium content are presented in a concise form. These models represent layered magma bodies with vertical compositional gradients (e.g., Hildreth, 1981; Bachmann and Bergantz, 2008). Representative models and their predicted MT responses are shown in Figure 12. Further details are provided in Appendix 2, including a summary of the various magma body models (see Table A2). Note that all resistivity models include a shallow conductive layer from the surface to 1 km depth below the surface that simulates the shallow low-resistivity layer (C1 in Fig. 2). This layer reduces the resolution of the underlying resistivity structure in MT exploration, and was included to make the modeling exercise more realistic.
Varying the Melt Composition, Melt Interconnection, and Melt Fraction
The modeling approach is shown in Figure 12. The temperature in Figure 12A increases linearly with a geothermal gradient of 65 °C/km until it reaches the magma body at a temperature of 980 °C (Sparks et al., 2008), where it becomes isothermal due to convection. The rock matrix resistivity is temperature dependent and was calculated with the equation of Hashim et al. (2013) (see Equation 1). The bulk resistivity is constant in the melt-free upper portion of the model because it equals the resistivity of the host rock (1000 Ωm), but begins to decreases as the increasing temperature makes conduction through the rock matrix possible. The pressure is calculated from the depth, assuming lithostatic conditions and a rock density of 2700 kg/m3.
The reference model (model 1) is shown in Figure 12 and is a uniform magma body from 15 to 25 km depth with a melt fraction of 25%, a silica content of 60 wt%, a water content of 6 wt%, a sodium content of 3 wt%, and interconnected melt with a cementation exponent of m = 1. All subsequent models involve a layered body with an upper portion from 15–20 km depth below the surface, and a lower portion from 20–25 km depth below the surface. The depths are fixed and the properties of each layer are varied. Variations in structure are considered in the following.
Variable silica content. Model 3 attempts to make the upper layer more resistive by making it more felsic (70% silica content) than the lower layer (60%). Both models 1 and 3 have a uniform water content of 6 wt%, and a uniform sodium content of 3 wt%. The bulk resistivity is shown in Figure 12G and it can be seen that the more felsic upper layer makes the top of the magma body slightly more resistive (increases from 2.7 to 4.6 Ωm, a factor <2). Variations in the silica content are petrologically possible (Hildreth, 1981), but changing just the silica content results in relatively minor changes in the resistivity of the magma body.
Water content. Model 5 attempts to make the upper part of the magma body more resistive by lowering the water content in the upper portion of the magma body (4 wt%). This produces a relatively small change in the resistivity from 2.7 to 6.9 Ωm, a factor of ∼2.5. The partitioning of water in this way is likely unrealistic because crystal settling should create a higher water concentration in the top portion (Bachmann and Bergantz, 2008).
Melt connectivity. Model 7 uses a combination of elevated silica content and decreased melt connectivity to raise the resistivity of the upper layer. The resistivity of the upper layer increases from 2.7 to 17.5 Ωm, producing a change of a factor of ∼6.5. A change in melt connectivity near the top of a magma body is expected from petrological studies (e.g., Hildreth, 1981; Bachmann and Bergantz, 2008), and this variation causes one of the largest resistivity changes found in the modeling exercise.
Sodium content. Model 10 attempts to increase the resistivity of the upper layer by decreasing the sodium content of the melt from 3 wt% to 1 wt%. This is geologically unlikely because a value of 3 wt% is already quite low (Muir et al., 2014). For model 10, the resistivity of the upper layer increases from 2.7 to 11.6 Ωm, a factor of ∼4.
The apparent resistivity and phases curves predicted for the resistivity models all exhibit the same shapes, but detectable variations are seen (Figs. 12H, 12I). Model 7, which combines physically realistic changes in both silica content and melt connectivity, exhibits the largest differences compared to model 1. As seen in Figures 12J and 12K, the apparent resistivity difference is larger than a factor of 2, and the phase difference is ∼10°, at periods >100 s. These differences are greater than noise levels found in typical MT data. For example, a noise level of 10% in apparent resistivity corresponds to a noise level of 3° in phase. These values are both below the differences in model responses shown in Figure 12.
The scenarios here show that a layered magma body with higher resistivities in the upper part and lower resistivities at depth could be produced by varying several properties. Varying the silica content from 70% to 60% with depth would produce a small change in surface MT data. Varying the melt connectivity (e.g., from interconnected, m = 1, to isolated, m = 2) produces a larger and detectable effect, as does changing the melt fraction. Changes in sodium content and water content may be petrologically unrealistic, but would produce similar observable changes. A combination of these parameters can produce larger changes (see Table A2), to a maximum of a factor of 16. A higher resistivity for the upper portion of the magma body could be detected by MT data, causing a greater depth to the top of a conductive layer to be inferred. However, it remains to be determined if these differences are resolvable within the uncertainties associated with MT data and inversion algorithms.
Seismic velocities for the layered models are expected to all be very similar, because they are relatively unaffected by such a layered composition. Therefore, if the layered magma body is detected as the top of a low-velocity zone, it may appear at a shallower depth in seismic data compared to MT data, which is more sensitive to the magma body composition.
MT Inversion Model for a Layered and Nonlayered Magma Body
The forward modeling exercise showed that compositional differences within a layered silicic magma body may produce detectable variations in the MT data that would be measured at the surface. A synthetic inversion is needed to test whether the resistivity values and geometry of the different models can be recovered. Random Gaussian noise (5%) was added to the forward modeled MT data and then MT inversions were run for both model 1 (nonlayered body) and model 7 (layered body). The synthetic profile had a vertical cell spacing of 500 m in the depth of interest. Figure 13A shows the synthetic resistivity models and Figure 13B shows the 1-D MT resistivity-depth profiles produced by the synthetic inversion. The depth to the top of the uniform magma body (i.e., the inflection point) is identified at a shallower depth than the layered magma body. Figures 13C–13F show the true models and the results of the synthetic inversion as vertical sections. The surface conductor masks the signal, causing the lower layers to be less well resolved. The results show for the nonlayered case, the depth to the top of the magma body is well resolved. However, for the layered case, the model shows the depth to the lower layer is well resolved because it is a strong conductor, but the depth to the top of the upper layer is not interpreted as readily, because it acts as a relative resistor.
Can a Layered Magma Body Explain the Depth Differences between Seismic and MT Measurements?
Both the forward modeling of MT data (Fig. 12) and the synthetic inversions (Fig. 13) clearly show that MT data can distinguish between a uniform magma body (model 1) and a layered magma body (e.g., model 7). The forward modeled MT data for several models showed significant differences, and these differences would typically be above the error of MT data. The synthetic inversion models show a depth difference to the interpretable top of the magma body between the layered magma body (model 7) and the uniform magma body (model 1). Previously the shallower depth to the top of the APMB determined from the ANT + RF method as compared to the MT method was explained as being due to resolution differences between the different geophysical methods.
The results give another potential explanation for the difference in depths to the top of the APMB from seismic ANT + RF and MT results; they can be explained by a compositionally layered magma body. If the APMB has a silica-rich and noninterconnected upper layer this can drastically affect the resistivity and therefore the MT data, showing a greater depth to the top of a conductive layer. However, the seismic velocities will be relatively unaffected by this layering, showing a shallower depth to the top of a low-velocity zone, resulting in a difference in depths between the two methodologies.
CONSTRAINTS ON VISCOSITY
Viscosity is an important physical property in volcanology because it controls the speed at which magma moves through the Earth and predicts the explosivity and potential hazards of volcanoes (e.g., Pommier et al., 2013; Dingwell, 2006). Determining the viscosity of a melt in the laboratory is relatively straightforward (Dingwell and Hess, 1998); however, when considering the structure at depth in the crust, geophysical constraints must be used. Ryder et al. (2007) used postseismic crustal deformation to estimate viscosity, and Rippe and Unsworth (2010) were able to link MT resistivity measurements and crustal viscosity by estimating crustal flow due to a pressure gradient from surface topography. Resistivity and viscosity are both related to the amount of fluid in a rock, and measurements of resistivity can be used to constrain viscosity (Rippe and Unsworth 2010; Pommier et al., 2013). Both electrical conductivity and viscosity are parameters that depend on the chemical composition of the sample, but viscosity has a much higher sensitivity, with small changes in any parameter creating changes of several logarithmic units for viscosity.
While the electrical conductivity of melt is controlled by the presence of charge carriers such as sodium ions, it is the silica framework of a melt that controls the viscosity (Pommier et al., 2013). An increase in silica content raises the viscosity of melt, because the structure becomes more rigid and its ability to flow is reduced. The viscosity of a melt has an inverse relation with the temperature, because the melt begins to crystallize as it cools, limiting its ability to flow (Dingwell, 2006). The viscosity also has an inverse relationship with the water content. Water acts to make the silica more mobile and decreases the viscosity (Dingwell, 2006). In comparison, electrical resistivity depends on silica content, water content, and temperature in a similar way as viscosity; the resistivity decreases as water content or temperature is increased, and the resistivity increases as silica content is increased.
A two-stage calculation is needed to compute the viscosity from geophysical measurements. In the first stage, the viscosity of the liquid melt is calculated. Composition values for andesite and dacite samples from Volcán Uturuncu from Sparks et al. (2008) and Muir et al. (2014) were used in this analysis. Two approaches are used to calculate the viscosity of a pure melt. Method A uses the relation of Giordano et al. (2008), which depends on the temperature, composition (primarily silica and sodium), and water content. Method B uses the relation of Pommier et al. (2013), which also includes the electrical conductivity of the melt as a parameter. Note that the conductivity of the melt is dependent on the temperature, the composition, and the water content (Pommier and Le Trong, 2011).
The second stage of the calculation is to consider the overall viscosity of a partial melt, i.e., a mixture of pure melt and solid particles, such as crystal rock grains. Semiempirical relations derived from laboratory studies show that the viscosity increases as the crystal fraction is increased (Caricchi et al., 2007; Costa, 2005; Costa et al., 2009). Experimental results can be plotted as relative viscosity values by comparing the apparent viscosity of the sample (measured) with the viscosity of the pure melt (with no crystals), which allows different experimental setups and different sample compositions to be compared. A crystal fraction (melt fraction = 1 – crystal fraction) of 35%–50% can increase the viscosity of the melt by a factor of 10; a crystal fraction of 60% (melt fraction of 40%) can increase the viscosity of the melt by a factor of 102–104; a crystal fraction of 80% (melt fraction of 20%) can increase the viscosity of the melt by a factor of 106–108 (Caricchi et al., 2007).
Method A—Relation of Giordano et al. (2008)
Using the model of Giordano et al. (2008) and the melt composition values from Sparks et al. (2008), the viscosity of the pure andesite melt in the deeper magma reservoir of the APMB is predicted to be ∼103 Pa·s (Fig. 14A). A viscosity of 105.5 Pa·s is predicted for the melt in the shallow magma reservoir. The bulk viscosity of the partial melt in both the APMB and the shallow upper magma reservoir can then be calculated using the Caricchi et al. (2007) model, which shows that if the magma reservoir has a melt fraction of 20% (crystal fraction of 80%) then the partial melt viscosity will increase by a factor of 106–108. This gives a final effective viscosity for the partial melt in the shallow magma reservoir of 1011.5–1013.5 Pa·s for 20% melt fraction. Correspondingly, the partial melt in the deeper magma reservoir will have an effective viscosity of 109–1011 Pa·s for 20% melt fraction.
Method B—Relation of Pommier et al. (2013)
Viscosity and conductivity are theoretically related through the modified Stokes-Einstein equation (Pommier et al., 2013). However, a relation that clearly expresses the dependence on measurable quantities from the rock composition is useful for geophysical interpretations. Pommier et al. (2013) proposed a relation between the viscosity of melt and electrical conductivity, for a given melt composition, water content, and temperature. This relation is based on empirical evidence from laboratory studies of silicate melt conductivity and viscosity models (Pommier and Le Trong, 2011). The composition is quantified through the optical basicity (OB) parameter, which is derived from estimates of the oxide ion activities (Pommier et al., 2013; Duffy and Ingram, 1975; Zhang and Chou, 2010).
Figure 14B shows the predicted viscosity for various resistivity values of a melt that has a variable composition (water content = 3–6 wt%; temperature = 900–980 °C; OB = 0.59–0.64, silica content of ∼60%–65%). The viscosity can vary by many factors of 10 (e.g., 1010) if the resistivity is varied by only a factor of ∼2. Using a melt resistivity value of 0.2 Ωm, as determined for the APMB (see Fig. 11; cf. values from Laumonier et al., 2014), the relation gives a viscosity of 106 Pa·s for the APMB. For a temperature of 900 °C (the lowest acceptable value for the relation) and a resistivity of ∼1 Ωm to represent the shallower body (Laumonier et al., 2014), a viscosity value of 108 Pa·s is calculated. These values are higher (by ∼2 orders of magnitude) than those discussed here using the model of Giordano et al. (2008).
The viscosity of the partial melt in both the APMB and the shallow upper magma reservoir can be calculated using the model of Caricchi et al. (2007). This gives an effective viscosity for the partial melt in the shallow magma reservoir of 1014–1016 Pa·s for a 20% melt fraction. Correspondingly, the partial melt in the deeper magma reservoir will have an effective viscosity of 1012–1014 Pa·s for a 20% melt fraction.
How Do the Viscosity Values Compare?
Assuming a minimum melt fraction of 20%, the effective viscosity of the APMB has been calculated to be a maximum of 1014 Pa·s, while the shallow magma reservoir has a maximum effective viscosity of 1016 Pa·s. These values are in good agreement with the Fialko and Pearse (2012) geodynamic model, which estimated a viscosity of 1016–1018 Pa·s for the APMB and the magma in the crust above it. The addition of rheological information is useful when interpreting the conductivity anomalies from MT data, because it can provide constraints on magma motion in the crust above the APMB. These viscosity values will later be used to compare with proposed geodynamic models.
IMPLICATIONS FOR DEFORMATION MODELS OF THE APMB
A range of geodetic models have been proposed to explain the pattern of surface deformation around Volcán Uturuncu observed with InSAR (interferometric synthetic aperture radar). Figure 15 shows several of these models overlain with the 3-D resistivity model obtained from the MT data. Details of the geodetic models are as follows: model 1 is a spherical inflation source centered at 17 km BSL (Pritchard and Simons, 2004), model 2 is a horizontal ellipsoid source centered at 19 km BSL (Pritchard and Simons, 2004), model 3 is a vertical ellipsoid source centered at 16 km BSL (Hickey et al., 2013), model 4 is a diapir extending vertically from the APMB at a depth of 15 km BSL (Fialko and Pearse, 2012), and model 5 is a flat-topped body rising from the APMB (Walter and Motagh, 2014). The inflation sources are located within the APMB, or slightly below it, as determined from seismic models (e.g., Zandt et al., 2003; Ward et al., 2014) and MT resistivity models.
Can the MT resistivity model distinguish between the proposed inflation models? The MT resistivity model may not allow us to distinguish between the geometry of model 1 and model 2, but it can be seen that both these inflation sources are inconsistent with the resistivity model. Model 3 proposes that the surface deformation is caused by a vertically elongated inflation source driven by the lateral motion of magma (Hickey et al., 2013). This model predicts a shallowing in the upper surface of the APMB that is consistent with the 3-D resistivity model (see Fig. 2). Volcán Uturuncu is located above the highest point of a long-wavelength upwelling of the APMB (see Comeau et al., 2015). This upward bulge may be caused by the inflation source. The 3-D resistivity model shows that, at a local scale, the MT results are consistent with a diapir geometry (such as model 4), and fit the dimensions proposed by Fialko and Pearse (2012). Furthermore, viscosity calculations show that the MT derived resistivity model is consistent with model 4. The MT model shows that the diapir is more resistive than the APMB, which may imply that it has a lower melt fraction, or a different melt composition than the APMB. A flat-topped magma reservoir extending above the top of the APMB does not match the MT data at the scale proposed by Walter and Motagh (2014). It is possible that this geometry would be difficult to resolve from model 4 at smaller scales, so the MT model does not preclude this geometry.
The deformation source may be radially symmetric and located within the APMB. However, above a depth of 10 km BSL the magma distribution pattern does not appear to be radially symmetric. Instead the MT model shows that magma and fluids are moving upward in a series of discrete dike-like zones striking approximately east-west, as opposed to the broad homogeneous upwelling expected from a diapir.
EVIDENCE FOR A SHALLOW MAGMA CHAMBER BELOW VOLCÁN UTURUNCU
InSAR data imply a deep inflation source that is believed to be related to the APMB (e.g., Henderson and Pritchard, 2013; Fialko and Pearse, 2012), but a shallow magma body does not appear to contribute to the deformation pattern. The deformation signal of the shallow magma chamber could easily be hidden by the larger signal of the APMB (Sparks et al., 2008). The mineralogy observed in dacite samples from Volcán Uturuncu gives evidence for the existence of a shallow magma chamber (de Silva and Gosnold, 2007). Petrological evidence for magma storage at shallow depths (4–10 km) has been seen at other volcanoes (Blundy and Cashman, 2005; Bachmann and Bergantz, 2008). Sparks et al. (2008) proposed that dacite magmas derived from the APMB ascend periodically to form shallow magma bodies at depths of 4–8 km below the surface throughout the APVC, and that these shallow magma bodies supply magma during eruptions. The brittle-ductile transition zone is shallow below the APVC (5 km depth below surface, near sea level; Jay et al., 2012) due to the thermal heating of the upper crust, which explains why magma ascends and ponds at these shallow levels. Muir et al. (2014) identified a preeruptive magma storage location below Volcán Uturuncu at a depth of 2–6 km below surface from the geochemistry of erupted dacites (S1 in Fig. 15).
To date, the only geophysical evidence of a shallow magma chamber below Volcán Uturuncu has come from a local seismometer network that imaged a low-velocity zone at 2–6 km depth below surface (Jay et al., 2012) (S2 in Fig. 15). The small low-velocity zone was also spatially correlated with earthquake hypocenters (see Fig. 15). However, the low-velocity zone is one continuous feature, and separate dikes cannot be identified, as seen with MT data. The 3-D MT resistivity model shows a low-resistivity feature at a depth of 3–7 km below surface (C3 in Figs. 15 and 2), which correlates well with the position of the assumed shallow magma chamber.
Correlation of Electrical Resistivity and Earthquake Locations
The resistivity models derived from MT data show that the volcano-tectonic earthquake hypocenters below Volcán Uturuncu are clustered at the perimeters of the conductive dike-like feature (C3; see Fig. 2). Volcano-tectonic earthquakes are thought to occur because of stresses from the migration of magma through the crust, for example as dike propagation, or from the injection of new magma into the crust from a deeper source (Roman and Cashman, 2006). Both horizontal slices and vertical sections through the 3-D resistivity model show that earthquakes surround the shallow conductive body below Volcán Uturuncu that is thought to represent a magma storage region (C3).
A number of other MT studies have observed that seismicity occurs at the perimeters of conductors that are interpreted as fluid-bearing zones (e.g., Wannamaker et al., 2009; Ogawa et al., 2001) or conduits for fluid migration (Bertrand et al., 2012). Turkoglu et al. (2008) noted that two mechanisms are possible to link electrical resistivity due to fluids and seismic activity: either fluids are controlling crustal deformation, causing seismicity by weakening the crust, or deformation produces a network of interconnected cracks enhancing the effect of fluids on electrical resistivity. Below Volcán Uturuncu it is likely the first mechanism that is dominant.
Summary of Depths to the APMB from Different Methods
Table 1 provides a summary of published depths to the APMB. As discussed here, various geophysical methods have detected the APMB. Several seismic methods detected a low-velocity zone below the Altiplano that was interpreted as a region of partial melt. A high-conductivity zone was detected with MT data, and also interpreted as a zone of partial melt. InSAR data show a ground deformation pattern around Volcán Uturuncu that can be modeled with source depths and geometries that match other evidence. Table 2 provides a summary of depths to the shallow preeruptive magma chamber. Seismic ANT results, earthquake hypocenters, geochemical data, and MT data all agree that a shallow chamber is located below Volcán Uturuncu at a depth near sea level, or just above it. Geodetic modeling and InSAR methods have not detected this chamber, probably because any signal would likely be dominated by the larger APMB.
The inversion model derived from the MT data gives new constraints on the distribution of magma beneath Volcán Uturuncu. The MT model shows that the APMB is a horizontally continuous regional-scale magma body, containing >15% andesitic melt. The top of the APMB appears to shallow beneath Volcán Uturuncu, with pathways of low resistivity connecting the APMB to the surface near Volcán Uturuncu. These pathways indicate an active magmatic system, as magma and aqueous fluids are required to explain the low resistivity.
The deep area of high conductivity (C2) also closely coincides with the area of very high attenuation reported by Haberland et al. (2003). High temperatures alone can explain the attenuation results, but cannot explain the low-resistivity values observed, while electrical conductors such as graphite can explain the low-resistivity values, but would not affect the seismic attenuation. Therefore, these correlated results (high conductivity, low velocity, high attenuation) jointly support the interpretation of interconnected partial melts, which would significantly increase the conductivity, lower the velocity, and increase the attenuation, giving stronger evidence of a large crustal magma body.
The MT data presented in this paper represent the first large-scale, high-resolution, broadband MT study of the Altiplano-Puna plateau and the only MT study of Volcán Uturuncu to date. Conclusions derived from this study include the following.
The depth to the top of APMB as determined from the 3-D electrical resistivity model is 14 +1/–3 km BSL (16–20 km below surface). This depth is in reasonable agreement with seismic velocity models. The APMB has an estimated thickness of 6–11 km.
The existence of a shallow (2–7 km below surface) magma storage region below Volcán Uturuncu, as inferred by petrological studies (Sparks et al., 2008; Muir et al., 2014), is supported by the MT resistivity models that show a conductor (∼5 Ωm) at that location, the Volcán Uturuncu high-conductivity zone. Three conductive dikes, oriented approximately east-west, are imaged with the MT data; each has a volume in excess of 500 km3.
The geometry of the upper surface of the APMB is characterized by an upward bulge beneath Volcán Uturuncu where it is shallowest. Pathways of low resistivity connecting the APMB to the surface below Volcán Uturuncu are imaged with the MT data.
The resistivity of the APMB requires that it contains a minimum of 15% interconnected andesitic melt. The shallow conductor below Volcán Uturuncu cannot be explained with dacite melt alone, but requires the addition of saline aqueous fluids to explain the observed resistivity values, hinting at a possible hydrothermal origin, as proposed by Jay et al. (2012).
Correlation of the geometry of conductors and volcano-tectonic earthquakes below Volcán Uturuncu supports the interpretation of fluid or magma migration through the crust (e.g., as dike propagation) and the injection of new magma from a deeper source below.
A comparison and correlation of the electrical resistivity models and seismic velocity models revealed that areas of low-velocity and low-resistivity were coincident with the APMB location.
Synthetic modeling showed that MT data can distinguish between a uniform and a layered magma body. A layered magma body that has a more resistive, silica-rich, and noninterconnected top layer would cause the MT method to image the top of a high-conductivity zone (e.g., the APMB) at a greater depth than the low-velocity zone imaged with seismic methods.
The viscosity of the APMB for 20% partial melt is calculated to be <1016 Pa·s. This is consistent with the geodynamic model of Fialko and Pearse (2012), who proposed that a diapir could be rising from the APMB.
This research was supported by a Natural Sciences and Engineering Research Council of Canada Discovery grant to Martyn Unsworth and National Science Foundation grant EAR-0908281 to Cornell University (Ithaca, New York). We thank Weerachai Siripunvaraporn for the three-dimensional (3-D) inversion algorithm, and Compute Canada for access to the WestGrid computer clusters that made the 3-D inversions possible. We thank the following for assistance in the field: Djamil Al Halbouni, Grover Arismendi, Herman Binder, Hernan Barcelona, Enrique Calvo, Victor Cuezzo, Daniel Diaz, Monica Doorenbos, Matthew Drew, Alicia Favetto, Benjamin Heit, Sergio Leon, Pablo Maciel, Gisel Peri, Cristina Pomposiello, Dauth Porcil, Santiago Ramos, Dante Scatolon, Mayel Sunagua, Eusebio Ticona, Faustino Ticona, Pablo Vargas, Fernando Zamudio, and Marcos Zeballos. This work has benefited from discussions with many colleagues at the University of Alberta, including Benjamin Lee and Robert Ferner.
APPENDIX 1. MAGNETOTELLURIC INVERSION AND MODEL RESOLUTION TESTS
Many combinations of inversion parameters were investigated and the main features of the resistivity models did not depend on any specific parameters choices, implying that the features are robust and are required by the magnetotelluric (MT) data. The inversion parameters tested are summarized in Table A1. Figure A1 shows maps of the impedance phase at representative periods, illustrating the data fit. The model fits the observed data well. Note that at a period of 1.2 s both the xy and yx impedance components appear similar because the electrical resistivity structure is relatively one dimensional. Differences at a period of 150 s arise due to the presence of the dike-like east-west conductors.
APPENDIX 2. INTEGRATED PETROLOGICAL-GEOPHYSICAL CONTRAINTS ON THE STRUCTURE OF THE ALTIPLANO-PUNA MAGMA BODY
We showed that both the three-dimensional (3-D) magnetotelluric (MT) resistivity model and the seismic velocity model detect a low-resistivity–low-velocity layer at a comparable depth, which can be identified as the Altiplano-Puna magma body (APMB). The agreement in depth for the two geophysical anomalies is relatively good, and the differences may be within the data uncertainty, or due to resolution differences between the two methods. However, it is possible that the properties of the APMB vary with depth and the two geophysical methods are sensing different layers of the APMB. If the upper part of the APMB was characterized by a low-velocity zone, but not a low-resistivity zone, then a difference in depth would be expected. This could be achieved by varying the composition of the upper portion of the APMB in such a way that the resistivity was increased but the seismic velocity was unchanged. In the following, a model of the APMB is developed and a range of properties investigated to determine if a realistic layered model could explain the differences between the ANT (ambient noise tomography) and MT models of the APMB. We present here and a more detailed explanation of the summary of the analysis in the text.
Seismic Velocity Structure
Note that the velocity is primarily sensitive to the amount of fluid or melt. Takei (2000) showed that there was some dependence on the melt geometry, through the dihedral angle. As the dihedral angle becomes smaller, the melt becomes more interconnected and this reduces the shear velocity through a reduction in the rigidity of the rock matrix. Shear velocity is quite insensitive to the composition of the melt (Dvorkin, 2008). Bottinga and Weill (1970) gave equations that allow the calculation of magma density based on its chemical composition and temperature. Using the composition and temperature values from Sparks et al. (2008), the density of andesitic magma is ∼2400 kg/m3. Assuming a rock matrix density of 2700 kg/m3 (del Potro et al., 2013) and a rock matrix velocity of 3.5 km/s, the shear wave velocity can be computed for any melt fraction. Using the relation of relation of Dvorkin (2008), the shear wave velocity structure of model 1 (section B of Fig. A2) is shown in section C of Figure A2. The velocity is constant throughout the magma body, so long as the melt fraction is constant. A layer with 25% melt fraction results in a velocity of ∼1.5 km/s from the Dvorkin relation, which matches the ANT + RF (receiver function) results of Ward et al. (2014). A 15% melt fraction gives a velocity of ∼2.2 km/s, still within the values for the low-velocity zone as defined by Ward et al. (2014).
Electrical Resistivity Structure of a Layered Magma Body
Model 1 is the reference model with a melt fraction of 25% in a uniform magma body extending from 15 to 25 km depth, a silica content of 60 wt%, water content of 6 wt%, a sodium content of 3 wt%, and interconnected melt with a cementation exponent of m = 1. All other models involve a layered magma body with upper and lower portions in the depth ranges 15–20 km and 20–25 km below surface. The models used in the forward modeling study are listed in Table A2. To compute the bulk resistivity, the modified Archie’s Law (Glover et al., 2000) was used with 25% partial melt and the melt was assumed to be well connected with a cementation exponent of m = 1. The bulk resistivity is shown in Figure 12G, and decreases significantly within the magma body. Variations in structure are considered in the following.
Variable Silica Content
Figure A3 shows several models that explore the effect of variations in the silica content of the melt. Three scenarios are considered: model 1 (the reference model) represents a uniform magma body with a silica content of 60%; model 2 represents a layered magma body where the upper portion is more felsic, having 65% silica, and the lower portion has 60% silica content; model 3 is also layered, with 70% silica content in the upper portion and 60% silica in the lower portion. These three models have uniform water content of 6 wt% and a uniform sodium content of 3 wt%. Model 1 shows a constant resistivity within the magma body, while models 2 and 3 predict that resistivity will decrease downward through the magma body, because the upper portion with a higher silica content is more resistive. The upper portion of the magma body in model 3 is nearly twice as resistive as in model 1 (4.6 Ωm compared to 2.7 Ωm). This shows that the upper portion of the magma body could be more resistive than the lower portion if there is a vertical change in the silica content of the melt (e.g., 60% silica versus 70% silica).
The next stage is to determine if these resistivity models produce detectable differences in the MT data that would be measured at the surface. Forward modeling was used to calculate the predicted 1-D MT responses for these resistivity models, as shown in Figures A3, sections F and G. The apparent resistivity curves predicted for the three models are all quite similar. Figure A3, section H shows the ratio of the resistivity values for each model compared to model 1, and Figure A3, section I shows the phase difference for each model and model 1. At long periods (>1000 s) model 3 shows a resistivity difference to as much as a factor of 1.4 larger (10 Ωm versus 14 Ωm) from model 1 and from 100 to 1000 s a phase difference of 4° is seen. These differences are above the noise levels of typical MT data (e.g., 10% error in apparent resistivity for 30 Ωm = 3 Ωm; an absolute phase error of 3°), although they are small. In summary, Figure A3 shows that depth variations in silica content produce modest variations in the resistivity of the magma body.
Variable Water Content
Because water content can have a large effect on melt resistivity, partitioning of water within the magma body was investigated. The models considered in Figure A4 have water contents that vary with depth. Model 4 has 9 wt% water in the upper portion and 6 wt% water in the lower portion, which produces a lower resistivity in the upper portion. Note that this trend is the opposite of that required to make the upper portion of the magma body more resistive. Making the upper part of the magma body more resistive requires a lower water content, as illustrated in model 5, with 4 wt% water in the upper portion and 6 wt% water in the lower portion.
The predicted MT responses are shown in Figure A4, sections F and G, and there are observable differences between models 1, 4, and 5. The apparent resistivity curves for model 1 and model 4 differ by a factor of 5, while model 1 and model 5 differ by a factor of 1.7 (Fig. A4, section H). The phase curves for models 1 and 4 have a difference of as much as 22°, while models 1 and 5 have a difference of as much as 6.5° (Fig. A4, section I). Note that water will preferentially partition into the upper layers of the magma body (Hildreth, 1981). Thus water partitioning is not an obvious explanation for a magma body that has higher resistivity at the top and lower resistivity at the base, as is required to have the MT-imaged APMB appear deeper than the seismically defined APMB. However, it may be possible to have discrete stacked magma bodies with differing water content, rather than a single body with a gradient in water concentration (Annen et al., 2006). In this way the upper body may have a lower water content, increasing the resistivity of the upper portion of the magma body.
Effect of Melt Interconnection
Another factor that could produce a magma body with a resistivity that decreases with depth is the degree of melt connection, and this is investigated in Figure A5. If the melt is not well connected, the bulk electrical resistivity will be higher. Model 6 is identical to model 1, but with melt in the upper portion of the magma body in isolated inclusions, as represented by a cementation exponent of m = 2 in the modified Archie’s Law. The resistivity depth profile in Figure A5, section E shows a higher resistivity in the upper portion of the magma body.
Model 7 is more realistic and includes variations in both silica content and melt interconnection, with 70% silica in the upper portion and isolated melt (m = 2). Model 7 shows a large difference in apparent resistivity compared to model 1 of a factor of 2.2, and a phase difference of as much as 9° as seen in Figures A5 and A5, sections H and I.
Variation in Melt Fraction
Variations in melt fraction as a function of depth are also considered, and are shown in Figure A5. These are certainly expected in magma bodies, based on petrological studies (Hildreth, 1981). The melt fraction would be expected to increase upward within the magma body as crystal settling occurs. Model 8 explores this effect and has a melt fraction of 15% in the upper portion and 25% in the lower portion. This produces a resistivity depth profile with a downward decrease in resistivity from 4.5 Ωm to 2.7 Ωm. Model 8 shows a difference in apparent resistivity of a factor of 1.4 compared to model 1, and a phase difference of as much as 4°.
As shown here, shear wave velocity depends primarily on the amount of melt. Therefore, if the melt fraction varies with depth, the seismic velocity will change with depth. The relation of Dvorkin (2008) predicts a velocity of 1.5 km/s for a 25% melt fraction and a velocity of 2.2 km/s for a 15% melt fraction, using the values in Figure A2. These velocities are both below the limit for the low-velocity zone of Ward et al. (2014) that was defined as Vs < 2.9 km/s.
Thus a downward increase in melt fraction will produce a downward decrease in resistivity, and make the MT-imaged magma body appear deeper. However, this increase in melt fraction, if large enough, will also increase the seismic velocity and so place the seismically defined APMB at greater depth. Thus variations in melt fraction alone may not be able to explain the greater depth of the MT-imaged APMB.
Variation in Sodium Concentration
The electrical resistivity of melt is sensitive to the sodium content (Pommier, 2014), and so changes in the sodium content between an upper and lower layer could produce a magma body with decreasing resistivity with depth, and this is investigated in Figure A6. Models 9 and 10 are identical to model 1, but with the sodium content in the upper portion of the magma body (2 wt% and 1 wt%) lower than that in the lower portion (3 wt%). The resistivity depth profile in Figure A6, section E shows a higher resistivity in the upper portion of the magma body. Model 10 shows a difference in apparent resistivity compared to model 1 of a factor of 2, and a phase difference of as much as 8°. Model 11 includes a variation in silica content (70% silica in the upper portion). Model 11 shows a difference in apparent resistivity compared to model 1 of a factor of more than 2.2, and a phase difference of as much as 9°.