We have carried out a tomographic inversion for seismic velocity in the vicinity of Uturuncu volcano (Bolivia) based on a 33-station temporary seismic network deployment. We combine travel times from earthquakes in the shallow crust with those from earthquakes on the subducting Nazca plate to broadly constrain velocities throughout the crust using the LOTOS tomography algorithm. The reliability and resolution of the tomography is verified using a series of tests on real and synthetic data. The resulting three-dimensional distributions of Vp, Vs, and Vp/Vs reveal a large tooth-shaped anomaly rooted in the deep crust and stopping abruptly 6 km below the surface. This feature exhibits very high values of Vp/Vs (up to 2.0) extending to ∼80 km depth. To explain the relationship of this anomaly with the surface uplift observed in interferometric synthetic aperture radar (InSAR) data, we propose two scenarios. In the first, the feature is a pathway for liquid volatiles that convert to gas, due to decompression, at ∼6 km depth, causing a volume increase. This expansion drives seismicity in the overlying crust. In the second model, this anomaly is a buoyant pulse of magma within the batholith, ascending due to gravitational instability. We propose a simplified numerical simulation to demonstrate how this second model generally supports many of the observations. We conclude that both of these scenarios might be valid and complement each other for the Uturuncu case. Based on joint analysis of the tomography results and available geochemical and petrological information, we have constructed a model of the Uturuncu magma system that illustrates the main stages of phase transitions and melting.

Massive caldera-forming eruptions represent rare but catastrophic events. Strong explosive eruptions in recent human history—including Mount Tambora in Indonesia (A.D. 1815) (Stothers, 1984) and Huaynaputina in Peru (A.D. 1601) (de Silva and Zielinski, 1998)—were locally devastating and impacted lives in other countries. Though the actual risk from eruptions of this type is modest, the scale of these eruptions makes them significant. They imprint geology around the world, capture the imagination of the public, and represent a formative process in continental earth science.

Perhaps more importantly, the scale of these eruptions indicates the existence of prodigious magmatic systems in the crust beneath them. Decades of field studies have demonstrated that, indeed, surface volcanism is but one manifestation of a larger subsurface process that is ultimately responsible for constructing much of the Earth’s continental crust. Crustal magmatism consists of both intrusive and extrusive processes. While there is little debate that both intrusive and extrusive magmatism contribute to the formation of continental crust, the relationship between the two is less clear (e.g., Bachmann et al., 2007). Specifically, it is unclear whether large ignimbrite eruptions are necessarily paired with the formation of new plutons (e.g., Glazner et al., 2004). Alternatively, perhaps large intrusions are less prone to eruption because they thermally weaken the crust, making it more accommodating to in situ magmas (e.g., Jellinek and DePaolo, 2003). All of these reasons provide relevance for the study of volcanic systems capable of hosting large caldera-forming eruptions.

Uturuncu volcano (Bolivia) is one such location. It is part of the Altiplano-Puna volcanic complex (white dashed line in Fig. 1), above the central Andes subduction zone and the downgoing Nazca plate. It sits near the border juncture between Chile, Argentina, and Bolivia. In the late Miocene, a “flare-up” of ignimbrite volcanism is recorded throughout the Altiplano-Puna volcanic complex. Since 10 Ma, >30 caldera forming eruptions have occurred here. At least ten of these eruptions left a mark on the global scale (Salisbury et al., 2010). These ignimbrite strata cover most of the area around Uturuncu providing a clear, if ominous, reminder of eruptions in the geologic past.

Uturuncu lavas consist of andesites and dacites. Rock compositions show that dacites are likely formed during andesite fractional crystallization with segregation of noritic cumulates, and zonal phenocrysts of orthopyroxene and plagioclase indicate the subsequent mixing of dacitic and andesitic magmas (Sparks et al., 2008).

Argon-argon (39Ar/40Ar) dating show that the volcano was active 890–270 k.y. ago, and the eruption activity included lava doming and flows (Sparks et al., 2008; Salisbury et al., 2010). The influence of glaciers is seen on the summit lavas, confirming the absence of Holocene effusive activity. The lavas contain various xenoliths: gabbroids, cumulates, and basement rocks (hornfels, sandstones, limestones). Two fumarole fields are active at the Uturuncu summit. They produce significant amounts of sulfur, with gas temperatures <80 °C (Sparks et al., 2008). With the absence of modern cooling lava domes, the fumaroles are a primary indicator for an active magmatic reservoir.

Uturuncu has received considerable attention in recent years following the recognition that it was inflating at 1–2 cm/yr during the late 1990s in a manner consistent with magma accumulation at depths of 10–20 km (Pritchard and Simons, 2002). The radially symmetric inflation, surrounded by a ring of slight deflation, was later interpreted by Fialko and Pearse (2012) as evidence for a diapir sourced by the mid-crustal magma layer. The authors posited that this diapir accumulates magma beneath Uturuncu by drawing it out of the surrounding region, thus creating the long-wavelength subsidence as well as the narrower uplift feature.

This combination of subsidence and uplift changes strain patterns on existing fractures and in hydrothermal systems near sea level, generating volcano-tectonic earthquakes that often manifest as swarms. This was supported by Jay et al. (2012) who investigated the distribution of local seismicity and the volcano response to major regional and global earthquakes.

Based on receiver function analysis, Chmielowski et al. (1999) suggested an extensive magma body beneath the Altiplano-Puna volcanic complex. Using the same method, but newer station networks, Zandt et al. (2003) provided more arguments for the presence of the magma body and determined that it is associated with a very low-velocity layer. They suggested that this anomalous layer is a sill-structured magma body beneath the entire Altiplano-Puna volcanic complex, and proposed a special term for it—the Altiplano-Puna magma body. Based on gravity inversion, integrated with the available geophysical, geological, and petrological observations, del Potro et al. (2013) proposed that partially molten felsic bodies ascend as diapirs through the hot ductile middle to upper crust. They claimed specifically that these features are likely rooted in, and extend up from, the Altiplano-Puna magma body near Uturuncu. A joint inversion of surface wave data and receiver functions (Ward et al., 2014) suggests a more distributed zone on the order of 10 km thick and an estimated potential volume of ∼500,000 km3. Ambient noise tomography (Jay et al., 2012) suggested the potential for magma above sea level just north of the Uturuncu summit.

Each of these techniques provides one view of the magmatic system, shaped by the selective resolution of that particular approach. For example, ambient noise tomography images structure in the upper crust, but cannot retrieve features at deeper levels where the major magma sources are expected. Receiver function methods provide information on the interface geometry, but suffer from depth tradeoffs with velocity parameters.

The analysis presented in this study complements these previous efforts. Travel-time tomography, based on local earthquakes, provides three-dimensional (3-D) images of the P- and S-wave velocity (Vp and Vs) in a volume of crust, around and beneath Uturuncu. By including data from earthquakes occurring on the underlying subducting slab, this tomography is able to provide resolution in the mid- and lower crust. When evaluated in the context of existing studies, the tomography presented here provides a glimpse of the deep crustal roots of the Uturuncu system.

The seismic data were recorded as a part of the PLUTONS project from April 2010 to October 2012. The network consisted of 33 stations surrounding the volcano with an aperture of ∼100 km (Fig. 2). The network was operated in two parts, each for roughly 18 mo, such that the full network was in place for one year. All stations comprised three-component Guralp CMG-3T (120s) sensors digitized on Reftek RT-130 dataloggers at 100 samples per second.

An initial catalog of earthquakes for this time period was compiled using the automated approach described by Thompson and West (2010). Inclusion in the catalog required an earthquake to have detections consistent with P-wave phase arrivals on the vertical component seismograms at five or more stations. To emphasize local earthquakes and minimize the inclusion of teleseisms, detections were carried out on traces bandpass filtered on 3–20 Hz. Each set of arrivals was compared, via grid search, to the travel times predicted from a grid of hypocenters extended 0.5° in each direction from Uturuncu and to a depth of 200 km. If the detections matched predicted arrivals within a tolerance of 1 s, then the event was assigned a hypocenter corresponding to the best-fit grid point. Events were declared based solely on the presence of consistent P-wave arrivals. However, S waves detected on horizontal component seismograms using the same filters and detection criteria were associated with earthquakes if their arrival times fit the matching criteria. Automated magnitudes were assigned to each event based on a local magnitude (ML) algorithm. Readers are referred to Thompson and West (2010) for a complete treatment of the automated earthquake detection process.

This automated catalog of earthquakes was the starting point for creating a manually reviewed earthquake catalog. Because the tomography relies on high-quality arrival times, all phase arrivals were reviewed, and adjusted or deleted manually. P and S phases were added, where appropriate, for additional stations. New hypocenters were calculated using the GENLOC algorithm of Pavlis et al. (2004). Shallow earthquakes >80 km from Uturuncu are excluded from the catalog because they fall outside the aperture of the network and are poorly constrained. The crustal earthquakes cluster primarily in the top 15 km below sea level (Hutchinson, 2015). We augmented this catalog with earthquakes occurring on the subducting slab. The Nazca slab is a prolific source of seismicity across the study area. To integrate these earthquakes, we used the hypocenters assigned by the Observatorio San Calixto (OSC; Bolivia) in their national earthquake catalog for Bolivia. We then added P- and S-phase picks for the PLUTONS stations. While the Uturuncu stations do not constrain the slab events well (especially in depth), the steep and deep ray paths from these earthquakes hold the potential to image the crust much deeper than the local earthquakes. The locations assigned are based on the regional seismic network and are generally well constrained in location and depth. We included a selection of slab events within 150 km of Uturuncu. These earthquakes were manually selected from the OSC catalog to produce a geographically distributed set of slab earthquakes. Because the Nazca slab dips ∼35° to the east in this area, these earthquakes increase from 100 to 300 km depth eastward across the study area (Fig. 2).

For tomography, we further subset these earthquakes to include only those with at least eight total P- and S-wave picks. After performing preliminary earthquake locations, we rejected travel times with residuals >1 s and >1.5 s for the P and S waves, respectively. These threshold values represent possible residuals estimated using the expected values of the anomalies and maximum length of the rays. The selection resulted in a total of 7186 P-wave and 4579 S-wave ray paths from 677 earthquakes.

Earthquake location and the tomographic inversion for the P- and S-wave velocities were performed using the LOTOS code (Koulakov, 2009), which is freely available at http://www.ivan-art.com/science/LOTOS. The main input data for the algorithm are the arrival times of P and S waves and station coordinates. Additional requirements include a one-dimensional (1-D) seismic velocity reference model and several user-defined parameters for different stages of calculations.

The first step is the absolute location of earthquakes within an optimum 1-D velocity model. We use the grid search method described by Koulakov and Sobolev (2006), with a fixed Vp/Vs ratio and P-wave velocities defined at a variety of depth nodes. Between these levels, velocity values are linearly interpolated. Using a grid search algorithm to locate the earthquakes is much slower than inverting the travel times, but is also considerably more stable. This helps retain as much of the original data as possible. We started searching for an optimal 1-D model from P-velocity distributions derived by Jay et al. (2012) and Ward et al. (2014). Then we perturbed velocity and Vp/Vs values to minimize the travel-time residuals and retain as many phase picks as possible in the data set with residuals smaller than the predefined thresholds. The best model obtained after several iterations is presented in Table 1. This model serves as the reference velocity structure for the 3-D tomography. Note that during this optimization procedure, we considered models with and without low-velocity layers mimicking the Altiplano-Puna magma body. We did not however detect any significant difference in residuals. This suggested that the body-wave tomography would be relatively insensitive to broad 1-D variations in mid-crustal structure.

The main part of the nonlinear tomographic process consisted of iterative repetition of the source location and matrix inversion steps. The source location was performed in an updated 3-D velocity model using the bending algorithm of ray tracing (Um and Thurber, 1987). Unlike the preliminary 1-D location procedure, at this stage we solved for hypocenter locations using an iterative gradient descent method. This inversion is faster and more precise, but less stable, than the grid search method. The velocity structure was parameterized on a grid of nodes with spacing corresponding to the density of rays. In map view, the nodes were distributed regularly with a spacing of 3 km in areas with sufficient ray density (>10% of the average ray density value). In the vertical direction, the node spacing was inversely dependent on the ray density, but nowhere <3 km. To minimize artifacts related to the grid configuration, the inversions were performed independently on four separate grids, rotated in azimuth by 0°, 22°, 45°, and 67°. These four results were then averaged into a single model.

The first-derivative matrix includes the elements responsible for the P- and S-wave velocity distributions, source coordinates, and origin times. The inversion was performed using the LSQR method (Paige and Saunders, 1982; Nolet, 1987). The flattening of the velocity model was controlled by an additional matrix block damping the difference between velocity parameters in neighboring nodes. Increasing the weight of this block makes the anomalies smoother. To balance the flattening and the number of iterations, we fixed the number of iterations at five and tuned the stability of the solution by changing the flattening coefficients. The most important parameters used for inversion are presented in Table 2. These values were determined through testing on synthetic models.


We used a variety of synthetic tests to assess the reliability and resolution of the data set and parameterization. We tested the impact of random noise using a so-called odd-even test. We divided the entire experimental data set into two subsets (e.g., with odd- and even- numbered events) and performed independent inversions for both subsets. To the extent that random noise is an influential factor in the data, it should cause discrepancies between the two versions of the model. The results of this test for depth slices of 15 km and 30 km for the P- and S-wave anomalies are shown in Figure 3. For the P-wave velocity model, the resolved patterns are nearly identical. For the S-wave velocity model, we can identify some modest differences. For example, a small positive anomaly to the southwest of Uturuncu at 15 km depth is visible in one data set and is not restored in the other one. This indicates that the anomaly might not be stable and is quite possibly an artifact of noise in the travel-time data set.

Another set of tests measures the ability to recover synthetic anomalies. These tests are useful not only for checking spatial resolution but also for tuning the parameters that control the inversion (Table 2). Synthetic travel times were computed for the original source-receiver pairs using the 3-D ray-tracing algorithm based on the bending method (Um and Thurber, 1987). The synthetic travel times were then perturbed with random noise (in this case 0.1 s and 0.15 s for the P- and S-wave data). After computing the synthetic data, we “forgot” all information related to the velocity model and source locations. The reconstruction was conducted using exactly the same workflow as for the real data, including the determination of an initial source location in the reference 1-D velocity model.

To assess spatial variations in horizontal resolution, we used checkerboard tests with alternating positive and negative anomalies. In the version presented in Figure 4, 15 × 15 km anomalies are separated with unperturbed areas 5 km in width. The amplitudes are ±7% and have opposite signs for the P- and S-wave anomalies, creating contrasting patterns in the Vp/Vs ratio. These anomalies extend in depth as columnar features. The features are relatively well resolved at 15 km depth, recovering the basic checkerboard geometry across much of the region with amplitudes diminished by roughly 30%. At a depth of 30 km, we observe a fair degree of lateral smearing in the northwest-southeast direction. The Vp/Vs images, somewhat surprisingly, appear to be more robust than the P- and S-wave velocity images.

It is common in tomographic studies to encounter poorer resolution with depth as a result of tradeoffs between source depth and model velocity. For these data, we expect additional smearing in the vertical direction because of the generally steep rays provided by the slab earthquakes. Though Figure 2 demonstrates considerable ray density beneath Uturuncu to depths of ∼100 km, these rays are generally steep and subparallel. As a result, we anticipate model features being potentially elongated in the vertical direction. To assess the impact of vertical smearing, we consider two tests with synthetic models defined in the vertical section. In the first model, we consider two anomalies with the size of 15 × 15 km located in the depth interval between 5 km and 20 km depth (Fig. 5). The recovered models show that the upper limit of these anomalies at 5 km depth is fairly well recovered, whereas the lower limit at 20 km depth is strongly smeared. In the second case, we consider four anomalies located in two depth rows. The transition between the rows occurs between 10 and 15 km depth. We see that the anomalies are better recovered for the P-wave model; in the case of the S-wave model, the high-velocity anomaly is smeared with a dip to the west, and we observe the same for Vp/Vs ratio.

These tests indicate a fair ability to resolve features with length scales of 15 km or more, with the caveat that smearing may be a factor in some directions.

Results and Discussion

The final model reduced P-wave travel-time residuals from 0.143 to 0.089 s (38%) and the S-wave residuals from 0.282 to 0.155 s (45%). The higher reduction in S-wave travel-time residuals stands out considering the higher noise levels for S-wave phase picks. This can be explained by the higher sensitivity of the S waves to crustal properties. The travel-time residuals are larger, and the resulting anomalies are larger.

The P- and S-wave velocity anomalies, as well as Vp/Vs ratio, are presented in two horizontal sections at 15 and 30 km depth corresponding to the upper and mid-crust (Fig. 6) and in one vertical section cutting west-east through the Uturuncu edifice (Fig. 7). These images are relative to the starting 1-D model. In most areas, we observe that the Vp and Vs features generally correlate with each other, though with different amplitudes. For example, at 15 km depth, to the north of Uturuncu, we observe a positive (fast) velocity anomaly, which is stronger in the Vs model and less prominent for the Vp model. At 15 km depth, the Uturuncu complex seems to be associated with a circular low-velocity anomaly surrounding a local higher-velocity anomaly, which is especially clear for the Vp model.

In the vertical section, starting from the depth of ∼2 km below sea level (bsl), a local high-Vp pattern coexists with a strong low-Vs feature. This creates a very strong high in the Vp/Vs ratio reaching 2.0 at 15–20 km bsl. This feature is roughly shaped like a tooth, having a flat upper boundary and roots down to ∼80 km depth. As we see in synthetic tests in Figure 5, the limited vertical resolution might be responsible for downward smearing of crustal anomalies; thus, the observed tooth-shaped anomaly might be partially caused by this effect. At the same time, when comparing the results in Figures 5 and 7, we see that the downward propagation of synthetic anomalies is less important than that of the real ones, which indicates that this anomaly does really exist in the mantle. A flat upper boundary appears to be consistent with the flat-topped magma body suggested by Walter and Motagh (2014). The width of this feature is ∼20 km. The depth of the most-pronounced velocity anomalies corresponds to the location of the magma body responsible for the ground deformation as estimated by Fialko and Pearse (2012).

This combination of high Vp and very low Vs anomaly has been observed beneath several active volcanoes, including Mount Spurr (Alaska, USA; Koulakov et al., 2013b), Mount Redoubt (Alaska; Kasatkina et al., 2014), and Klyuchevskoy volcano (Kamchatka Peninsula, Russia; Koulakov et al., 2013a). This relationship might be explained by the intrusion of mafic magmas with a high degree of partial melting and volatiles (e.g., Takei, 2002). P-wave velocity is quite sensitive to composition, tending to manifest mafic magmas as high velocity (Mechie et al., 1994; Trampert et al., 2001; Behn and Kelemen, 2003). In contrast, the S-wave velocity is quite sensitive to the percentage and connectedness of any liquid phases (Karato, 1995; Hammond and Humphreys, 2000; Schutt and Lesher, 2006). As a result, magmas tend to manifest as negative Vs anomalies, regardless of composition.

On top of the large tooth-like feature, at ∼2–6 km bsl, there is a cluster of earthquakes. We observe this cluster on the resulting velocity anomaly model; earthquakes were registered during the seismic network was operated, in years 2010–2012. Above this cluster, the tomographic results show a decrease in both Vp and Vs, though the amplitude of the Vp anomaly is stronger. This results in a low Vp/Vs ratio. Note that a very similar observation exists at Mount Spurr corresponding to a period of unrest in A.D. 2004–2005 (Koulakov et al., 2013b). Beneath Mount Spurr, a deep columnar feature with high Vp/Vs ratio was overlain by a shallow zone of low Vp/Vs ratio. The transition between these two features occurred at ∼2 km bsl, coincident with a vigorous cluster of earthquakes. At Mount Spurr, this Vp/Vs change was explained by the decompression transition of liquid volatiles (high Vp/Vs) to gases (low Vp/Vs). As shown by experimental studies (e.g., Takei, 2002) and some field experiments (e.g., Chatterjee et al., 1985, De Siena et al., 2010), the Vp/Vs ratio is strongly sensitive to the pore content. Gas phases create low Vp/Vs ratio (e.g., Julian et al., 1996). The shallow seismicity at a few kilometers below sea level beneath both Mount Spurr and Uturuncu may reflect the considerable volume change as volatiles flash into a gas phase. In both the Mount Spurr and Uturuncu cases, the presence of the shallow low Vp/Vs anomaly and seismicity coincides with fumarole activity.

A second hypothesis explaining the relationship between the Vp/Vs and seismicity patterns might be associated with the growth of the magma body beneath Uturuncu. In this case, the deep high-Vp/Vs feature may represent very hot crust that contains a significant percentage of volatile-rich melts. Vertical ascent of such magmas would cause the overlying rock to deform mechanically, leading to increased seismicity. This hypothesis is supported by the observed surface inflation in the Uturuncu area measured with interferometric synthetic aperture radar (InSAR) data (Pritchard and Simons, 2002; Fialko and Pearse, 2012) and by the observed earthquake focal mechanisms (Hutchinson, 2015; Alvizuri and Tape, 2016).

To assess the reliability of these scenarios, we present a simplified mechanical model simulating the stress and deformation effect of a gravitationally unstable magma body. The calculations were performed using the COMSOL Multiphysics 4.4 software (https://www.comsol.com). We consider an axisymmetric linear elastic model with a radius of 100 km and a depth of 104 km. We define zero normal stresses on the top surface and zero normal displacements on the side and bottom surfaces. Density is set according to the geologic information on the Altiplano-Puna region (Prezzi et al., 2009). The elastic Lamé parameters are defined according to the P- and S-wave velocities (Table 1). The shape of the magmatic chamber is modeled as a smoothed cylinder with a radius of 10 km, height of 30 km, and density of 2650 kg/m3.

The displacement suggested by the model is qualitatively similar to observed and modeled surface uplift patterns (Fialko and Pearse, 2012). Polyansky et al. (2010) considered a thermo-mechanical viscoplastic dynamic model for a similarly sized magmatic diapir to obtain a cumulative surface uplift of hundreds of meters. The fracture areas (suggested by von Mises stress; Figs. 8, 9) are located not only in the uplift center but also in the peripheral ring-shaped zone. A similar ring-shaped distribution of the third principal stress was obtained by Walter and Motagh (2014) for a deflating source. Our simple elastic-only model neglects changes in crustal properties with depth. The high heat (Jay et al., 2012) in this region would quickly increase viscous and creep effects with depth, leading to a decrease in von Mises stresses. The resulting decrease in seismicity that would be expected with depth is consistent with the observed seismicity at Uturuncu.

The low Vs, high Vp/Vs, and lack of seismicity strongly suggest the presence of partial melt. Based on this, we suggest that the deformation source does, in fact, represent some type of localized magmatic reservoir. At the same time, the opposite signs of the Vp and Vs anomalies strongly suggest compositional differences between magma body and the surrounding rocks. Comeau et al. (2015) made the same conclusion based on magnetoteluric models. The main question then is the relationship between this local accumulation of magma and the broad regional feature of the Altiplano-Puna magma body. Some research (Lipman, 2007; de Silva and Gosnold, 2007; Annen, 2009) suggests a model of batholith formation as an amalgamation of magmatic chambers in volcanic provinces. According to this model, the Uturuncu magma body is part of the batholith. Any magmatic system beneath Uturuncu is dynamic and interacting continuously with the surrounding crust and other magmatic components. Geochemical and geochronological constrains by Muir et al. (2015) show that late-stage (<485 ka) erupted lavas were derived from a source that experienced extensive magma mixing or contamination. These authors envisaged mixing of a dacite magma with two different andesites. One of the andesitic components is characterized by high MgO and Cr contents, suggesting an influx of more-mafic magma. The isotopic composition of the Uturuncu lava domes and flows also implies that they originate from mingling between very different, but coeval, magmas. These magmas could be created by differentiating primitive island-arc basalts and through strong contamination by upper-crustal material (Michelfelder et al., 2014). This is evidence for long-lived replenishment of the Uturuncu magma reservoir. The interaction of high-temperature lower-crust magmas with dacitic and andesitic mush or with already solidified parts of the reservoir may lead to a decrease of the mush crystallinity or even partial melting of solidified parts of the reservoir (Fig. 9). The hot mix of intruding magma and siliceous partial melts at the base of the magmatic column become buoyant and penetrates upward into the magmatic mush (Huppert and Sparks, 1988; Couch et al., 2001; Burgisser and Bergantz, 2011; Muir et al., 2015). These remobilized buoyant magmas generate pressure on the reservoir roof. This combination of buoyant force and accompanying degassing offer an explanation for recent unrest at Uturuncu volcano (Sparks et al., 2008). The transient pressure associated with these remobilized buoyant magmas would create temporary uplift. If this is the case, then it is possible that subsidence might later be expected as these magmas solidify and lose buoyancy. In this model, the shallow layer of seismicity is caused by (1) transient stresses above the magma body caused by the arrival of buoyant magmas, and (2) pore-pressure increases caused by the degassing of these magmas. It is interesting to note that these processes may occur on different time scales. Seismicity related to degassing might occur long after magmas arrive and as they begin to cool. If so, then there is no reason to expect that uplift and seismicity should always be concurrent. As an example, though deformation observed in the 1990s had ceased by the time of this experiment, the observed seismicity was much more vigorous than in surrounding areas. It is possible that this seismicity is an after-effect, tied to degassing, of the inflation over the previous decade or more.

The mid-crustal feature observed here, which is suggestive of a magma feeding zone, appears very localized to Uturuncu. Uturuncu lies between large Miocene caldera complexes (Fig. 9), implying that the voluminous caldera eruptions in its vicinity have been common in the past (Salisbury et al., 2010). This kind of eruption is associated with fluid-rich magmas. At the same time, the eruptions of Uturuncu have never been voluminous. The total volume of the edifice is estimated at 50–89 km3 (Muir et al., 2015; Michelfelder et al., 2014). The presence of biotite and amphibole indicates that volatiles were present in magma, but the melts were probably not fluid rich. Although the parameters discussed and the scale of processes indicated are very large (the area of surface deformation, the volume of the magma and volatiles needed to support this deformation, the size and shape of the reservoir from the tomographic images), we cannot state unequivocally that these amount to evidence for a so-called super-eruption similar to those that took place in the late Miocene (Gregg et al., 2015). However, these parameters do suggest that a more modest eruption is at least plausible at some point in the future, and a better estimate of possible eruption volume would entail an analysis of all structural complexity, including time variations and approximation of the future system development (e.g., Gottsmann et al. [2017]). The numerical model proposed here demonstrates how such an eruption might lead to caldera formation (Fig. 8). The lateral distribution of von Mises stress shows a ring-like zone where increased fracturing should be anticipated. This feature is roughly 40–80 km in diameter, roughly consistent with calderas observed throughout the region. The shallow crust in this zone is where fractures should be most anticipated during an eruption. These ring fractures may facilitate caldera collapse and are a likely zone of melt penetration prior to eruptions.

We have cataloged earthquake activity and constructed a seismic velocity model with 3-D distributions of P- and S-wave velocities and Vp/Vs ratio beneath the Uturuncu area in the crust and uppermost mantle. The earthquake distributions reveal two independent features. The deep earthquakes associated with the subducting Nazca plate are located at depths from 100 km to 300 km beneath the study area. Travel times from these earthquakes provide constraints on seismic velocity structures in the mid- and lower crust. A second group of earthquakes is located in the upper crust and adds resolution in the uppermost part of the model. Based on a series of tests of real and synthetic travel-time data, we can postulate that the derived seismic models are reliable and horizontal resolution is good, whereas in the vertical direction, anomalies are smeared.

The most prominent feature is a large tooth-shaped anomaly with Vp/Vs >2.0 and with roots extending down to a depth of ∼80 km. The upper limit of this anomaly is flat and is located at ∼2 km bsl (∼6 km below the surrounding area, and ∼8 km below the Uturuncu summit). Above this feature, we observe low Vp/Vs and a vigorous cluster of laterally distributed earthquakes. The location of this anomaly coincides with the surface uplift identified in InSAR observations. We propose two possible scenarios explaining the link between the tooth-like body and surface deformations. In the first scenario, the columnar anomaly of high Vp/Vs represents a pathway of ascending partial melts. At a depth of ∼6 km below the surface, volatiles are released from the melt as a result of decompression. This causes considerable volume increase, surface uplift, and elevated seismicity in the uppermost crust. The second scenario posits a magma body ascending due to gravitational instability. We use a simple model to simulate the stress field caused by a buoyant body of realistic shape to derive a qualitative fit with the observed ground deformations. It is likely that components of both of these models exist at Uturuncu.

This study is supported by the U.S. National Science Foundation (grant 0909254) and the Russian Science Foundation (grant 14-17-00430). The authors are grateful to O.P. Polyansky for the consultation in numerical modeling of the observed body. We thank N.L. Dobretsov for suggestions and improvements. We are grateful to the Supercomputing Center of Novosibirsk State University for computational resources. We thank two anonymous reviewers for helpful and constructive feedback.