A set of 14 teleseismic earthquakes was studied to determine how wave propagation was affected by a presumed magma body beneath Uturuncu volcano, Bolivia. Teleseisms are suitable for study because they are relatively long period, contain purely P waves, and have near-vertical incidence angles. The number of events is small but the events have good signal-to-noise ratios and very similar waveforms for each event so that reliable measurements could be made of arrival times and amplitudes. Attenuation of amplitudes occurs in a NW-SE trend beneath the volcano, 14 by 34 km (long axis NW-SE). Calculated values of the quality factor Qp are an average of 12.4, with extreme values as low as 1.8. These calculations are based on the assumption that the highest amplitude observed is the “true” amplitude, and all others have been attenuated. The average thickness of the anomaly is 10.2 km, and the center is ∼20 km SE of the summit, within the area of surface uplift measured geodetically. Time delays of up to 0.8 s were also observed. The pattern of attenuation and relative time delays together showed four trends: fast and not attenuated (normal crust), slow and attenuated (partial melt), fast and attenuated (likely high fracture density), and slow but not attenuated (possible deep low Vp structure).

Back azimuth differences of up to 60° were observed. In nearly all cases, azimuths were rotated into directions parallel to local rock fabric, suggesting that shallow crustal properties affected near-surface wave propagation. Overall results suggest partial melt as high as 10%–20% in a region of varying thickness, low Bouguer gravity and resistivity, high Vp/Vs, persistent seismicity, and overlapping a locus of recent uplift.

Uturuncu volcano, Bolivia, has attracted considerable scientific attention over the past few years. It was one of only four volcanic centers to show deformation based on an interferometric synthetic aperture radar (InSAR) survey of over 900 South American volcanoes by Pritchard and Simons (2002). This survey was superseded by a later study showing nine deforming centers in the Central Andes (Henderson and Pritchard, 2013). Modeling showed that the volcano was inflating 1–2 cm/yr, though this rate has decreased with time (Henderson and Pritchard, 2014), from a source at mid-crustal depths of 17–20 km below sea level. A reconnaissance survey in 2003 showed a high rate of seismicity even though the volcano has not erupted for 270 ka (Sparks et al., 2008). A joint geophysical experiment between scientists from the United States and the United Kingdom, funded by National Science Foundation and Natural Environmental Research Council (NERC), respectively, was conducted from 2009 to 2014. Called PLUTONS (shortened from PLUTONNNSSSSSSSS; Probing Lazufre and Uturuncu TOgether: NSF, NERC, NSERC, Sergeotecmin, Sernageomin, Observatorio San Calixto, Universidad Nacional de Salta, Universidad Mayor San Andres, Universidad de PotoSi, SERNAP, Chilean Seismological Service, Universidad de San Juan), the project included geologic (Sparks et el., 2008; Perkins et al., 2016), gravity (del Potro et al., 2013), magnetotelluric (Comeau et al., 2015), geodetic (Pritchard and Simons, 2002; Henderson and Pritchard, 2013), and seismic (Jay et al., 2012; Ward et al., 2014; Kukarina et al., 2014) components. Here we present new results based on seismic investigations of teleseisms, distant earthquakes whose long-period P waves travel through the crust at near-vertical incidence.

Determining the interaction between seismic waves and the seismic low-velocity zone (and possible magma body) beneath Uturuncu volcano is important because such interactions affect the results of seismic imaging methods and can be used as additional sources of information to constrain the location and properties of partial melt. Previous studies show that velocity decreases and attenuation increases as the solidus is approached and partial melting begins to occur in laboratory samples (e.g., Sato et al., 1989). We anticipate that seismic waves will be slowed down as well as attenuated, if partial melt is encountered in the crust. Attenuation at volcanic centers can be caused by several factors such as presence of magma, hydrothermal systems, heterogeneous material properties, or a combination of these (e.g., Schurr et al., 2003). Therefore, if the inflation center is an area of magma injection, we expect seismic waves that have passed through the inflation source to become attenuated, and therefore have diminished amplitudes compared to those that have passed through normal crust. Likewise, a crustal magma body should decrease the values of the crustal P-wave quality factor, Qp, both directly in the body and possibly within the heated host rock surrounding the magma. While seismic velocity is decreased by the presence of partial melt, especially interconnected melt, it is not particularly sensitive to temperature and therefore won’t be significantly affected by rocks at near-solidus temperatures; this is where attenuation, which is sensitive to temperature, can add information (Lees, 2007). It has been documented that the presence of magma within the crust can bend seismic waves, thus resulting in distorted seismic raypaths given the locations of the hypocenters with respect to the stations (e.g., Steck and Prothero, 1993; Yeguas et al., 2011; Galluzzo and La Rocca, 2013). This can modify the results of seismic studies that rely upon detailed knowledge of the raypaths. In this paper, we examine velocities, attenuation, and directional features of teleseismic waveforms to place additional constraints on magma body properties, location, and shape.

Uturuncu is located in SW Bolivia at latitude 22.27°S and longitude 67.18°W. The volcano is also at the SW part of two major provinces, one based on geology (the Altiplano-Puna volcanic complex [APVC]) and the other based on the presence of a seismic feature at depth (the Altiplano-Puna magma body [APMB]). These are subfeatures of the much larger Altiplano-Puna plateau.

The Bolivian Altiplano is part of the 1800-km-long, 350–400-km-wide Altiplano-Puna plateau, enigmatically formed in the absence of continental collision (e.g., Allmendinger and Gubbels, 1996). Uplift began ca. 25 Ma, as the convergence rate between the Nazca and South American plates increased and subduction shallowed (Allmendinger et al., 1997). Presently, the subduction angle is 30° beneath the Altiplano, which has a lithospheric thickness, attributable to crustal thickening, of up to 150 km (Allmendinger et al., 1997). The average Vp/Vs ratio of the southern Altiplano is 1.80, suggesting that the seismic velocities have been affected by high temperatures and regions of partial melting (Yuan et al., 2000). These mainly lower the Vs (e.g., Nakajima et al., 2001).

The APVC and APMB are two distinct features of the Central Andes that are affected by stresses in addition to subduction. The APVC is a major silicic volcanic province, located between 21° and 24°S (333 km) and 66° and 69°W (307 km), and resulting from an ignimbrite “flare-up” beginning in the late Miocene. It is deemed to be active because of Late Pleistocene and younger volcanic activity and active low-temperature geothermal fields (de Silva, 1989; de Silva et al., 1994). It is the youngest but largest of several such ignimbrite fields in the Central Andes volcanic zone, covering an area of ∼60,000 km2 located at the transition from the ∼4 km high Altiplano south through the ∼5-km-high Puna (de Silva, 1989; Whitman et al., 1996).

The APMB was identified as a ∼20-km-deep, 1-km-thick zone of low seismic velocity, for which Vs < 1.0 km/s (Zandt et al., 2003). It mostly underlies the <7 Ma ignimbrite complexes of the APVC, among which Uturuncu is located, rather than the Quaternary arc volcanoes to the west. The APMB is coincident with an ultra-low velocity zone (ALVZ), a region characterized as causing about a 10%–20% reduction in seismic-wave velocity over a thickness of 10–20 km (Yuan et al., 2000). Schilling et al. (1997) discovered a region of high conductivity (1 S/m) from depths of ∼20 km to at least 60 km below the Western Cordillera, which they modeled as the result of a body of 14–27 vol% interconnected partial melt. Using temperature modeling, heat-flow densities, and eruptive history, Schilling and Partzsch (2001) determined that this melt most likely has a granite-like composition and a crustal origin. It is likely overlain by a thin film of saline fluid (Schilling et al., 2006). Gravity data analyzed by Schmitz et al. (1997) displayed a region below the magmatic arc of density 0.16–0.25 g/cm3 lower than that of the surrounding crust. They interpreted it as a magma body of 10–20 vol% melt that reduced P-wave velocity by 0.8–1.2 km/s. These geophysical studies all agree on ∼10–27 vol% partial melt.

The crust above the magma body is strongly anisotropic, with 20%–30% anisotropy in a 3-km-thick surface layer and 15%–20% anisotropy in the remaining crust between this surface layer and the magma body; the anisotropy of both crustal regions has a strike of between 300° and 330° with a tilt of 45° from vertical (slow S-wave orientation), producing significant azimuthal variations in the radial receiver function components (Leidig and Zandt, 2003; Zandt et al., 2003). Anisotropy likely results from alignment of a system of fluid-filled cracks oriented ENE/WSW and dipping to the north from 45° to 80° from horizontal (Leidig and Zandt, 2003), or a crustal strain partitioning into subhorizontal zones of high and low strain (Schilling et al., 2006). The presence of a large crustal magma body in the Andes, an area of strong compressive tectonics, suggests that the APVC may be experiencing a transition from vertical thickening to horizontal extension along the arc (Riller et al., 2001). Vertically, this magma body also forms at the boundary between upper-crustal imbrication and lower-crustal thickening (Yuan et al., 2000).

There are three existing seismic models for the geometry of the APMB. The first is that of Zandt et al. (2003), based on receiver functions, in which the magma body is at 19–20 km depth below the surface, 1 km thick, and is a compositionally uniform lens (see Fig. 1). Here, the APMB underlies an ∼60,000 km2 area of 3° in longitude and 2° in latitude (Zandt et al., 2003). The second model is the joint ambient noise tomography and receiver function result of Ward et al. (2014), in which the magma body is ∼150 km across and centered beneath the center of surface InSAR-detected deformation (Fig. 1). The thickness of the central, strongest part of the anomaly is ∼11 km, though there is an anomaly from 4 to 25 km below sea level. In this model, the APMB is shallowest (4 km depth to surface) below Uturuncu, and the volume of the magma body is given as ∼500,000 km3. In this model, the APMB contains ∼4%–25% partial melt (Ward et al., 2014), a value that is consistent with resistivity estimates of at least 20% partial melt (Comeau et al., 2015). The third is a model of local thickening in two rectangular blocks of ∼40 km deep by 20 km wide (shallower block) and ∼30 km deep by 10 km wide (deeper block) beneath the volcano, as recovered by Vp/Vs tomography from Kukarina et al. (2014; shown schematically in Fig. 1).

Local thickening of the APMB beneath Uturuncu volcano, as observed by Comeau et al. (2015), may be the result of a diapir (Fialko and Pearse, 2012). However, magnetotelluric data suggest that the magma that sources this APMB upwelling exhibits a greater resistivity than that of the APMB proper. This may reflect a different composition, differing melt connectivity, or a lower melt fraction in this source magma (Comeau et al., 2015). In addition, above 10 km below sea level, magnetotelluric studies show that crustal magma distribution is not symmetric and occurs in discrete bodies (Comeau et al., 2015; Comeau et al., 2016). In the magnetotelluric study, the APMB is located at 16–20 km below the surface, with a resistivity of <3 Ω m centered ∼3 km W of the summit and extending off eastward. There is a vertically elongated low-resistivity body directly beneath the volcano that can be attributed to the magma chamber discussed in Muir et al. (2014). The resistivity of this body suggests saline aqueous fluids in addition to 15% melt, as even pure melt is insufficiently conductive (Comeau et al., 2015; Comeau et al., 2016). Highlighting narrow, vertically elongated structures in the upper crust above the APMB, del Potro et al. (2013) modeled the APMB as a body of non-uniform thickness with diapirs extending toward the surface from the top of the magma body (Fig. 1). The body thickens beneath these diapirs and thins away from them. The model they favor shows a 25% melt fraction, though they agree that the APMB more likely contains varying melt fractions.


Uturuncu volcano (Fig. 2) is a 6008-m-tall stratovolcano located in the Central Volcanic Zone of the Andes. Uturuncu is located upon ∼70-km-thick crust ∼150 km above the top of the subducting slab, 50 km east of the axis of the Andean volcanic arc (Allmendinger et al., 1997). Local seismicity has been recorded since seismic networks were installed in 2009 (details below). The volcano was last active from 890 to 270 ka (Sparks et al., 2008), erupting magmas of dacitic and andesitic compositions. Current local, low-level thermal activity is characterized by two active fumarole fields near the volcano’s summit (Sparks et al., 2008). Despite its dormancy, recent InSAR studies show a maximum inflation rate of 1.5–2 cm/yr over an area of 70 km width (Pritchard and Simons, 2002) and which is surrounded by subsidence to create a sombrero-shaped deformation pattern (Fialko and Pearse, 2012; Henderson and Pritchard, 2013; Hickey et al., 2013). This inflation was initially modeled as a point source with a depth of 15–17 km below sea level located 3 km to the southwest of Uturuncu’s summit (Pritchard and Simons, 2002). However, current interpretations are that the deformation signature, with the associated ring of subsidence, is the result of buoyant magma at a depth of 19 km below the surface, bulging upward at the center of deformation, displacing hot and ductile crustal rocks which then flow aside and downward (Fialko and Pearse, 2012). In this model, magma is pulled laterally into the diaper at depth. Current studies show a source volume of decreased density with respect to the crust, a high Vp/Vs ratio of >1.9, and a depth to the top of the layer of ∼6 km below the surface (e.g., del Potro et al., 2013; Kukarina et al., 2014; McFarlin et al., 2014). The greater depth (16–20 km below surface) of the APMB’s magnetotelluric signature is attributable to possible layering of different compositions (i.e., magmas of differing resistivities) in the magma body, which would push the magnetotelluric depths to greater depths than the actual top of the magma body (Comeau et al., 2016). Ward et al. (2014) used the results of a combined receiver function and surface wave dispersion study (Ward et al., 2013) to determine the location and size of the APMB. The Ward et al. (2014) study is regional scale, smooths the anomalies over a broad lateral area while preserving vertical resolution, and focuses on results derived from velocities. Here we seek more detailed information on the smaller region around Uturuncu, as well as adding to the understanding of how the anomaly affects seismic amplitudes via attenuation.

Uturuncu volcano exhibits a remarkable lack of variation in composition throughout its eruptive history (Sparks et al., 2008; Muir et al., 2014). Although the volcano is compositionally dacitic, it shows no evidence of having erupted explosively (Sparks et al., 2008). This follows reasoning by de Silva (1989) that intrusion rate into the crust has a stronger effect on eruptive style than magma composition or volatile content. The composition of erupted products at Uturuncu is likely the result of both the formation of dacite by fractional crystallization of andesitic parent lava within the APMB and magma mixing in a shallow storage zone 5–7 km below Uturuncu’s summit (Muir et al., 2014). Their study informs on both the depth at which magma is stored preeruptively in the Andean crust beneath Uturuncu volcano and the composition of the melt, which impacts factors such as the density and bulk modulus of magma being encountered by seismic waves. Geothermometry indicates that temperature varies temporally and spatially within the magma body feeding Uturuncu and that the magma may experience large temperature fluctuations prior to eruption (Muir et al., 2014), suggesting that magma storage conditions are non-uniform.

Seismicity at Uturuncu consists of shallow, near sea level (∼5 km beneath the volcano’s summit), almost exclusively volcano-tectonic events; however, some events with low frequencies appear to be deeper (>20 km) events whose higher-frequency components have been attenuated (Jay et al., 2012; Kukarina et al., 2014). Seismic swarms were observed a few times a month, consisting of 5–60 events occurring over a time span of minutes to hours. One swarm appears to have been triggered by the 2010 Mw = 8.8 Maule earthquake. This suggests that the hydrothermal system at the volcano is metastable, being affected by stresses caused by magma accumulation at ∼20 km depth within the mid-crust (Jay et al., 2012). The calculated b-value is low, within the range for tectonic earthquakes at 0.64 ± 0.04 (Jay et al., 2012). Analysis of regional slab earthquakes showed variable amounts of S-wave attenuation (M. West, 2016, written commun.) suggesting that waves from different earthquakes experience varying amounts of attenuation in their respective paths (Jay et al., 2012). Jay et al. (2012) also used heat-flow density to determine that the depth to the brittle-ductile transition zone beneath Uturuncu is at approximately sea level (6 km below the summit), though it could be as deep as 8.5 km below sea level. Because Jay et al. (2012) favor the shallower depth, earthquake activity apparently persists into the inferred brittle-ductile transition zone.

The PLUTONS network (Fig. 2 and Table 1) consisted of 28 seismometers deployed between April 2010 and October 2012. Some stations of the PLUTONS array were co-located with previous ANDIVOLC stations (Jay et al., 2012). The first half of the network was installed in April 2010 and removed in March 2012; the second half was installed in April 2011 and removed in October 2012. Dates of operation and numbers of stations are given in Figure 3. Each station consisted of a Guralp CMG3T three-component broadband seismometer, a separate GPS clock, and a RefTek RT 130 Datalogger. Data were collected at sampling rates of 100 samples/s and 1 sample/s (same data but with different sampling frequencies) and were recovered during service runs every six months.

In the current work, we studied 14 teleseismic events from several subduction zones at different azimuths with respect to the network in Bolivia (Fig. 4 and Table 2). Deep events were chosen to ensure that the seismic waves have only traveled once through the crust and asthenosphere in the vicinity of Uturuncu and, thus, were not affected by the crust or asthenosphere near the earthquake sources. Slab events of magnitude 5 and greater are visible in the Uturuncu seismic data with good signal-to-noise ratios. The teleseismic waveforms generally have periods of several seconds; thus, the wavelengths are of the order of 10 km and not sensitive to fine structure. The incidence angles are steep (Fig. 5); so the P waves mainly sampled vertical structures. The teleseisms show systematic changes between waveform peak-to-peak amplitudes and time residuals as a function of different earthquake back azimuths. We used events from the European subduction zones (ESZ; one event, the Calabria subduction zone), >80 km depth for the Kermadec-Tonga subduction zones (KTSZ; four events), >400 km depth for the Japan subduction zone (JSZ; four events), and >100 km for the South Sandwich Islands subduction zone (SSSZ; five events). Results from the ESZ were difficult to constrain given relatively low signal-to-noise ratio and therefore only included in the discussion of amplitudes.

We used a combination of MATLAB and Antelope (Boulder Real Time Technologies [BRTT]; www.brtt.com) to calculate the P-wave peak-to-peak and zero-to-peak amplitudes of the vertical component at each station for each earthquake. For an event from Japan (JSZ2a and JSZ2b; Table 2), we analyzed both PKIKP and PKP phases to determine if a differing angle of incidence and travel path gave different results. To ensure that the same P-wave phase of each teleseismic earthquake was analyzed for every station, we matched the phase by aligning waveforms using cross correlations of a 10–20 s window of data filtered from 0.375 Hz to 1.5 Hz across the network and common waveform characteristics, such as wave shape and frequency. Traces in which the phase could not be identified were discarded manually. Correlation coefficients ranged mostly from 0.72 to 0.949 (Table 3 and Table S11) with 1 being perfect correlation. One station, PLCO, had systematically low values. We then compared our calculated amplitudes across the network to determine if, for a given source, there is a corresponding area of diminished amplitude where the waves have passed through an attenuating body. As a control, and because few studies like this are known to us, we also performed this method using a teleseismic event from Bolivia on Transportable Array (TA) stations in Florida (Fig. S1 [see footnote 1]) to show how larger-scale variations could be resolved with this technique. Details are given in the discussion below and in the Supplemental Materials (see footnote 1). We identified seismic phases through TauP modeling (University of South Carolina [USC]; www.seis.sc.edu/taup/) of arrival times. The phases are shown in Table 2. Examples of the seismograms for one event are shown in Figure 6; seismograms from all other events are shown in the Supplemental Materials (see footnote 1).

To calculate the P-wave quality factor (Qp) values within the upper 30 km of each ray’s travel path, we used the formula:
where R is distance traveled in the attenuating medium (20 km for this study), f is frequency of the wave, V is velocity (average of local velocity model of 5.15 km/s), A is the amplitude at each station, and A0 is the maximum amplitude measured for the earthquake. We assume the wavefronts are planar at teleseismic distances; hence, we ignore geometric spreading. The calculations are based on the assumption that the highest amplitude observed is the “true” amplitude, and all others have been attenuated. Thus, this formula results in an apparent infinite Qp value for the station with largest amplitude (A = A0); so that value was not included in subsequent analyses. To determine if crustal Qp values are frequency dependent, we analyzed the amplitudes of each waveform at different filters in octave steps from 0.375 Hz to 1.5 Hz (0.375–0.75 Hz and 0.75–1.5 Hz), and then over the full range of 0.375–1.5 Hz. Then, we used the relationship between Qp, ln(A/A0), and f (the center frequency of the bandpass filter) to determine how a change in f affects Qp. This evaluation was done using the waveform suite for MATLAB (Reyes and West, 2011).

The P waves travel through the region of interest at depth (presumed to be a magma body) but also travel through the local shallow crust. H. McFarlin (2014, personal commun.) followed the methodology of Frankel (1982) and used local shallow earthquake data to determine the average Qp value for the local crust and obtained a range of values from 60 to 770, with most values (16 out of 26 total measurements) falling within the range of 130–240. There were no systematic trends observed in these data. Hence the shallow crust causes a modest reduction in all amplitudes, but we infer that the main deviations are caused by the presence of a magma body at depth.

We also determined relative arrival times for each event, assuming straight rays through a 2-D velocity model. To determine delay times, here defined as the deviation in relative arrival times from a plane wave arriving at the geographically closest station earliest, we first created a plane using the coordinates of the station closest to the earthquake epicenter as well as the known back azimuth and angle of incidence. Then, we calculated the shortest distance between the plane and each station using the point-plane distance calculated from the Hessian Normal Form of the plane. Using a 2-D P-wave velocity model modified from Ward et al. (2014), we determined the time expected for this plane wave to reach each station as our expected times. We subtracted these expected times from the actual arrival times. This procedure returned delay values that are positive for a late arrival, negative for an early arrival, and always zero at the station closest to the epicenter. The absolute value of the lowest negative value was then added to each delay time to determine the delay with respect to the station at which the earthquake arrived the earliest compared to the expected time. We interpret this station to be the one least affected by the magma body, which has been imaged to underlie the entire array.

From these delay values, we determined percent partial melt using the velocities of P waves in melt and in solid rock. We used a density of 2300 kg/m3, the value determined as the density of the region of partial melt by del Potro et al. (2013), and the partial melt values using this density are reported in this paper. We also used the crustal density value of 2800 kg/m3 (Perkins et al., 2016), though this value only decreased the partial melt percentage necessary by 1%–3% for a thickness of 20 km (increasing slightly with increasing percent melt). We employed typical crustal values of rigidity (26 GPa) and bulk modulus (45 GPa), which give velocities that agree with Uturuncu’s velocity model. For the melt, we used a bulk modulus of 45 GPa but a rigidity of 0.5 GPa, and for the model including water, we calculated water as 8% of the melt fraction, with a bulk modulus of 2.2 GPa and a rigidity of 0.1 GPa. These values are for pure melt and pure water in a silicate melt, allowing us to mix the moduli for melt (with water) and calculate velocities for a range from a purely solid crust to 100% melt. To determine the velocity for each percent partial melt, we computed the Hashin-Shtrikman-Wapole bounds (Hashin and Shtrikman, 1963; Mavko et al., 2009) to calculate the minimum and maximum rigidity and bulk moduli bounds as a function of percent melt. These values, combined with the density of the melt body, were used to calculate P-wave velocities for each integer melt percentage. Using this calculation, we obtained velocity values from 5.6 km/s for 0% melt and 4.5 km/s for 100% melt for a dry melt. Incorporating 8 wt% water into the magma, we obtained velocities between 5.6 km/s and 2.9 km/s. We used the lower bounds because crustal melts are likely to be interconnected (Schilling and Partzsch, 2001; Comeau et al., 2016). This was done for a mixture of melt and solid rock, then for this combination but with of water constituting 8 wt% of the melt fraction (Sparks et al., 2008). Using rigidities of 0 GPa for the melt and water fractions results in a lower rigidity bound of 0 GPa and, therefore, when using this to calculate velocities, causes the velocity to sharply decrease with the inclusion of any the fluid. This is our reason for using nonzero rigidities for the water and melt. Then, the thickness of the magma body was given (we tested values of 1, 10, and 20 km), and we calculated the time taken to travel this distance with a velocity equivalent to that in solid rock plus the time delay to determine our control values. Then, we used the velocities through melt to calculate the travel time expected for the above distance at each percent partial melt. For each station, the travel time value that matched to within a tolerance threshold (half of the difference between travel times for adjacent percent partial melt values) of the control value for that station gave the average percent partial melt value of the magma body encountered along that raypath. The results are given in Table 4. Additionally, we were able to solve for magma body thickness given percent partial melt values by using the equivalent partial melt velocities. As a test, we used depth values from sea level to the upper surface of the magma body calculated by receiver functions (McFarlin et al., 2014) and used the bottom depth of the Ward et al. (2014) model (25 km below sea level) to determine a magma body thickness and calculate percent partial melt values for this magma body geometry (Table 4).

Initial angles of incidence were calculated using the locations of the earthquake hypocenters with relation to Uturuncu volcano in the center of the seismic array. However, when calculating delay times for events from the KTSZ and SSSZ, we discovered a dependence on lag time with distance, suggesting that the calculated angle of incidence differs from the actual incoming raypath. Using the correlation matrix methodology of Jurkevics (1988), this issue was resolved, and therefore these calculated angles of incidence were favored in the analysis.

Finally, we determined the back azimuth—the angle clockwise from horizontal from the station to the event—for each of the above teleseismic events. We calculated the peak-to-peak amplitudes of each component for all stations. Using these values, we calculated each back azimuth as arctan (N/E) for events from the ESZ and KTSZ and arctan (E/N) for events from the JSZ and SSSZ, where N is the amplitude on the north-south component, and E is the amplitude on the east-west component. Then, we compared these measured azimuths to the expected back azimuth given the earthquake hypocentral location relative to the volcano.

We discovered a persistent region of reduced amplitudes and delayed arrival times centered near the station PLAR, SE of the volcano (Fig. 7) when analyzing teleseisms from the JSZ and SSSZ, from the NW and SE of the volcano, respectively. For these events, the attenuating zone is a 14 by 34 km (preferred values based on station distribution) elliptical region southeast of the summit centered between stations PLAR and PLSS (Fig. 7). This is 31° clockwise and 6° counterclockwise off of the median back azimuth for JSZ back-azimuth sets AJSZ and BJSZ, respectively, and 5° counterclockwise and 27° clockwise off of those for the SSSZ sets ASSSZ and BSSSZ (Fig. 7). In addition, there is a region of low amplitude continuing to the northwest of the volcano in the events coming from the SSSZ. This is ∼20 km by 23 km, centered on a line between station PLJR and Uturuncu’s peak, and is 3° clockwise off of the median back azimuth for SSSZ set ASSSZ and 33° clockwise off for set BSSSZ (Fig. 7). The signals analyzed had frequencies of 0.4–1.33 Hz, with corresponding wavelengths on the order of 4–12 km; so no topographic effects were considered. The incoming azimuthal angle differences are large (Fig. 8) and suggest significant anisotropy or preferred fabric either from the subduction zone, the mantle wedge, or the crust. The differences between the maximum and minimum peak-to-peak amplitudes varied between an amplitude loss of 46% (KTSZ2) and 90% (SSSZ4).

Qp values (Table 5 and Table S2 [see footnote 1]) were calculated for a 20 km and 25 km thick attenuating layer. The following Qp results are from the 20-km-thick assumption. They follow the same spatial pattern as the amplitudes, with an overall mean (excluding infinite values) of 12.4. The three stations with consistently small Qp values (between 2.3 and 6.5, 2.9–6.7, and 1.8–8.9, respectively) are PLSS, PLAR, and PLBR, all to the south and/or southeast of the volcano. Stations PLLL and PL07, also to the S/SE of the volcano, respectively, give mean Qp values of 7.1 and 5.5, both well below the overall mean. The waveforms for stations with the smallest values of Qp have undergone the most attenuation and therefore have the smallest signal relative to the maximum signal amplitude.

Analyzing two phases for one event—PKIKP and PKP phases for event JSZ2—did not show any significant differences in the distribution of amplitudes. The pattern was still the same, with a region of low amplitudes centered to the southeast of the volcano. For both phases, this region—including stations PLAR, PLSS, and PLLL—contained the smallest amplitudes in the network.

We infer that the number of arrival time picks is insufficient to gain much insight from a full inversion; however, relative arrivals help determine the geometry and properties of the potential magma body. The accuracy of relative arrival time differences is high at a few hundredths of a second for these long-period waves. Relative residuals vary from 0 to 0.8 s. We plot normalized amplitude versus time residual for event KTSZ2 (incoming from the SW) in Figure 9A. There is no trend; however, this showed four quadrants, divided by arbitrary reference lines: fast and not attenuated (upper left), slow and attenuated (lower right), fast and attenuated (upper right), and slow but not attenuated (lower left). Stations from each quadrant show spatial clustering in map view. Figure 9B shows that the cluster of stations S/SE of the vent, PLBR, PLSS, PLLL, PLO7, and PLAR, all show both delayed arrival times and amplitude attenuation. This particular combination is consistent with the presence of partial melt (e.g., Schmeling, 1985; Sato et al., 1989; Lees, 2007), and this location S/SE of the vent shows evidence for slow and attenuated waves for teleseisms from the JSZ, KTSZ, and SSSZ. Stations to the NW show early arrivals and moderate attenuation (likely high fracture density; see below), whereas those SW show early arrivals and no attenuation, suggesting intact country rock. Finally, the group of stations to the NE shows time delays but no attenuation. Several effects are likely occurring (see discussion below). Similar pairs of plots for other azimuths are shown in the Supplemental Materials (see footnote 1).

The combination of high attenuation (low Qp) and velocity reduction strongly suggests that partial melt is present and agrees with laboratory studies (e.g., Sato et al., 1989) and theoretical work (e.g., Chu et al., 2010). Velocity can be modeled independently as a first-order check; this assumes either a solid (compressibility and rigidity) or melt (compressibility only) and varying percentages of each. Results are shown in Table 4 where partial melt values are given for different thicknesses to account for the observed velocity reductions. Using the results of Chu et al. (2010) for a fluid-saturated porous material (granite, rhyolite melt, water, and CO2) yields 15.2% and 14.7% velocity reductions for stations PLLL and PLAR, respectively, corresponding to ∼8% and 7% (19 and 13, if excluding the influence of water) partial melt (Table 4). Laboratory studies of lherzolite and peridotite with no fluids (Sato et al., 1989) show low values of Qp of 10 or less depending on pressure at high temperatures >1155°C. Combining velocity and Qp requires temperatures of 1250°C at 0.5 GPa to give a partial melt of ∼15%.

The back-azimuth calculations (Fig. 8) showed remarkable variations between the expected azimuths (those calculated from the hypocentral locations of the events relative to the location of the volcano) and those that we determined from the waveform data. The maximum variation was 60° (Fig. 8), and the minimum was 1°. The means for the JSZ, KTSZ, ESZ, and SSSZ are –19.6°, +19.9°, –1.2°, and +17.6° (Table 6), respectively (+ is clockwise and – is counterclockwise). Although part of this variation can be the result of 1°–5° uncertainties in alignment between the N-S component and true north, the observed variations are on average an order of magnitude larger. The events from the NE (ESZ) and SW (KTSZ) both rotate into closer alignment with the NE trend of crustal fabric (perpendicular to anisotropy), suggesting that whatever is causing the S-wave anisotropy is also affecting the wave transmission for the long-period teleseismic P phases. Anisotropy studies use S waves, but they reveal rock fabric that affects P waves as well. The events coming from the NW (JSZ) and SE (SSSZ) exhibit preferential rotations of the calculated azimuths to the NE with respect to the expected azimuths. We interpret this to be caused mainly by properties of the shallow crust above the magma body, specifically a system of fluid-filled cracks oriented NNE-SSW (e.g., Leidig and Zandt, 2003).

We infer the presence of a magma body based on prior results by others including seismic Vp/Vs ratios, Bouguer gravity anomalies, magnetotelluric and deformation data, receiver function analysis, and the already established APMB. Our data can best be used to place additional constraints on the geometry and physical properties of the inferred magma body, including percent partial melt.

Teleseisms versus Regional Events

There are several reasons that we chose to study teleseismic events. The first is that the events are far away; hence the arrivals can be represented as plane waves. This greatly simplifies the geometry. Regional slab events, by contrast, produce curved wavefronts that are more difficult to model. Second, the teleseismic locations are known and are independent of the uncertain local velocity model, which affects slab event locations. Third, the teleseisms sample a very small part of the radiation pattern, whereas radiation pattern effects would be strong for slab events because of the closer distances and greater range of take-off angles. Fourth, the long wavelengths (4–12 km) of teleseisms means they are not sensitive to small-scale (<1 km) structures. Seismograms from slab events with shorter wavelengths show ringing from local structures. Finally, the teleseisms are pure P waves and have nearly vertical incidence, again providing simplicity. The study of slab events would be complementary and should be done. They would potentially include both P and S waves and may have favorable geometry and/or higher resolution assuming the modeling factors above are properly taken into account.

Attenuation and Qp Values

There are several mechanisms involved in attenuation, among them damped resonance or dissipation due to relative movements along grain boundaries, viscous relaxation and fluid flow in pores, and scattering at inhomogeneities (Jackson and Anderson, 1970; Richards and Menke, 1983; Sato et al., 1989; Watanabe and Sassa, 1996). Both attenuation and velocity of seismic waves change near the solidus as well as when partial melt is present, with attenuation increasing (Qp decreasing) and velocity decreasing. This condition corresponds to the lower right part of Figure 9A (red symbols). Conversely, intact crust has low attenuation (high Qp) and higher velocities as in the upper left of Figure 9A (blue symbols). Scattering by inhomogeneities may explain high attenuation but with higher velocities retained (e.g., Richards and Menke, 1980), and further, faults and fractured zones have been identified using attenuation tomography (Watanabe and Sassa, 1996). This corresponds to the lower left of Figure 9A (green symbols). Finally, the upper right portion of Figure 9A (black symbols) shows low velocities but low attenuation. This may be caused by lithology or other causes, but the low attenuation suggests that little partial melt is present.

Schurr et al. (2003) performed a regional-scale study of the subduction zone in the Central Andes. They found low values of Qp (∼80) in the vicinity of the active volcanoes to the west, but resolution near Uturuncu was poor. Low Qp (also ∼80) was also observed beneath Cerro Tuzgle to the SSE. Similarly, Haberland and Rietbrock (2001) determined Qp values of ∼100 near the border between Bolivia and Chile (to the west and southwest of our study), which they interpreted to be the magma source feeding the APVC. Our calculated values of Qp are quite low, an average of 12.4 with a few single-digit values. These are consistent with loss of amplitude of 46%–90% over just a few wavelengths. Similar low values have been observed in the laboratory (Sato et al., 1989) associated with partial melting of peridotite.

The distribution of amplitudes at Uturuncu can be viewed in several different ways to best identify anomalous regions. These amplitudes were shown above in Figure 7 with symbols corresponding to the direct observations. To highlight the amplitude differences regardless of source azimuth, a summary plot showing averages of the best three events from each direction (except the ESZ) is shown in Figure 10. The same highly attenuating zone SE of Uturuncu is observed.

Although our calculated Qp values are quite low, values such as these are not unprecedented. Schlue et al. (1996) have determined Q values for both P and S of ∼30 at frequencies of ∼1 Hz for the Socorro magma body. Similarly, Wilcock et al. (1995) determined Q values of, at minimum, 10–20 at the East Pacific Rise, which they attributed to thickening of a layer that is interpreted to be high-porosity lava flows and pillows (e.g., Hooft et al., 1996). Similarly, a shallow intrusion of gas-rich magma into the flank of Mount Etna exhibited Qp values of 10–30 (Martínez-Arévalo et al., 2005). Wilcock et al. (1995) ascribe Q values of 20–30 to “no more than a few percent melt” in this basaltic environment. There are also non-magmatic examples of fluid in pore spaces bringing apparent Qp to single-digit values (e.g., Korneev et al., 2004). Even without partial melt, raising the temperature of basalts and gabbros to near the solidus can give Q values of 20–40 at 10 Hz, and Q decreases from there with the introduction of partial melt (e.g., Wilcock et al., 1992, 1995). Additionally, increasing the thickness of the attenuating layer would increase our calculated values of Qp.

Validation of the Method

We are aware of only a few studies of teleseismic amplitudes similar to this one (e.g., Ward and Young, 1980). Therefore, we tested the efficacy of our teleseismic amplitude method by applying it to a completely independent region. Using the Temporary Array (TA) stations in Florida (Fig. S1 [see footnote 1]), we can show that this method of comparing raw amplitudes of teleseismic body wave arrivals is sensitive to large-scale variations as well as smaller-scale ones, such as those observed at our network at Uturuncu volcano. Because the Florida stations are more geographically distributed than those in Bolivia, we first checked to determine that the Florida rays didn’t cross any nodal planes in the focal sphere or result from focusing of different waves on a travel-time curve. All the waves are Pdiff from 106° to 112°. Additionally, we determined that the distances between source and stations wouldn’t result in interference from multiple contemporaneous arrivals, which could increase waveform amplitude. In the case of Florida, the peninsula is composed of relatively uniform thick Quaternary and Tertiary limestone units with small amounts of sandstone, whereas the northern region—the region showing larger amplitudes—is mainly composed of Miocene siliciclastics and carbonate sediments (South Florida Information Access; sofia.usgs.gov). Seismic energy is retained traveling through older, well-consolidated rocks to the north and northwest, along the panhandle, and lost traveling through the layered rocks of the peninsula. The greater attenuation in the peninsula results from both scattering as the seismic energy travels through the layered rock and encountering fluids, such as the numerous aquifers in the peninsula of the state. This result is strong enough to overcome the effect of distance on seismic amplitude, as the stations located on the peninsula are the closest to the epicenter but show the lowest amplitudes. Additionally, the majority of the sinkholes in Florida occur from the middle of the state (i.e., Orlando) to the north, with few occurring on the panhandle (Florida Department of Environmental Protection; dep.state.fl.us). Therefore, karst features, such as sinkholes and caves, may cause scattering as well. As such, while our Uturuncu network shows large variations in an area of diameter <100 km, the majority of the peninsula of Florida shows a lack of variation in amplitudes and low overall amplitudes over a much larger area (Fig. S1 [see footnote 1]). However, these observations reveal that the method allows determination of large-scale features that are correlated with crustal geology.

Frequency Dependence

Crustal Qp for teleseismic earthquakes may show frequency dependence in the 0.5–3.0 Hz band (Bache et al., 1986); all processing in this paper was done using a filter of 0.375–1.5 Hz, therefore intersecting this band. This suggests that within this band, Q increases with increasing frequency (Morozov, 2008). We see this frequency dependence when we compare the calculated percentage of energy difference between the minimum and maximum amplitudes for each teleseism to the period of the waveform, though it manifests as just a slight difference (within ±7 values of Qp and decreasing to <2 with decreasing Qp). To limit this effect, estimates of Qp (Tables 5 and S2 [see footnote 1]) from our attenuation data are only presented processed under a filter of 0.375–1.5 Hz, and only the dominant frequency of each waveform within this band is used in calculations of Qp. This is to ensure that variations in Qp within each event are mostly attributed to the physical parameters of the crust beneath Uturuncu.

Estimates of Partial Melt

Many previous investigations of attenuation and partial melt have mainly focused on the upper mantle and mid-ocean ridge systems, as well as laboratory studies. These benefit from the fact that upper mantle and mid-ocean ridge basalt (MORB) compositions are fairly uniform worldwide. In contrast, volcanoes have significant structural and compositional variations, as well as significant concentrations of fluids. Nevertheless, it is common practice to extrapolate from laboratory studies to real-world situations (e.g., Sato et al., 1989). Velocity studies, particularly of Vp/Vs ratios, often rely on partial melt to explain anomalies, as well as presence of fluids such as water or CO2. We did not take into account varying melt geometries (with the exception of assuming interconnected melt) or dynamic melting (Li and Weidner, 2013).

We believe that the partial melt and thickness estimates of the magma body derived from incorporating 8 wt% water into the melt provide the best model for this case. While the water content may vary spatially, there is ample evidence to suggest that crustal magmas in the APVC contain water (e.g., Sparks et al., 2008; Muir et al., 2014), with Laumonier et al. (2017) claiming a minimum of 8 wt%. Additionally, there is evidence of a hydrothermal system at Uturuncu (Sparks et al., 2008; Jay et al., 2012), and the observed magnetotelluric anomalies require the presence of another conducting fluid (i.e., saline aqueous fluid), as an unrealistic melt fraction would be necessary to result in a 3–7 Ωm anomaly (Comeau et al., 2016).

Comparing the estimates for partial melt and magma body thicknesses (Table S3 [see footnote 1]) with the depths to the top of the APMB given by receiver functions (McFarlin et al., 2014; Table S4 [see footnote 1]), we see some correlation. While the station overlying the shallowest part of the APMB, PLAN, commonly plots in the upper left quadrant (i.e., Fig. 9), stations PL07, PLAR, PLBR, PLLL, and PLSS all overlie a shallower-than-average part of the modeled APMB. In addition, station PLAR is over one standard deviation shallower than the average depth to the surface of the magma body, and station PL07 is just slightly within one standard deviation (by 0.011) from the mean. The fact that these two data sets don’t completely agree on the thickness of the melt layer could mean that in addition to a slightly thicker APMB beneath the stations to the SE of the volcano, there could be a larger percent melt and/or greater water content.

Directional Effects

The upper surface of the APMB has been modeled to be irregular (McFarlin et al., 2014); hence, we cannot rule out the possible effects of focusing and/or defocusing of rays. Because the center of the deformation signature, located ∼3 km to the SW of Uturuncu’s summit, may show convex-upward curvature as the center of a diapiric bulge of the APMB, defocusing of seismic energy as the wavefront passes from the APMB to the surrounding rock is a possibility for this location (e.g., Sheriff, 1975). This would reduce the amplitude for stations above the bulge because of diverging rays. Station PLSS, SSE of the vent, shows strong amplitude reduction, and PLBR, SW of the vent, shows modest reduction. The data from these two stations may be affected by bending. However, none of the teleseisms studied shows similar effects at all four close stations PLBR, PLCM, PLMK, and PLSS; therefore, we downweight these possible effects.

Our attenuation data show evidence of anisotropy in the magma body and vicinity, a feature also revealed in the receiver functions of Zandt et al. (2003). There are differences in attenuation patterns between waveforms coming from the NW/SE and those coming from the NE/SW (Fig. 7), both in shape and magnitude. The teleseisms sourced to the NW and SE (JSZ and SSSZ, respectively) show smaller regions of greatly diminished amplitudes—a region to the southwest of the volcano and a corridor centered on the volcano and extending in the NW-SE direction, respectively. This differs from the broad zone of slightly reduced amplitudes sweeping counterclockwise around the volcano from the west through the southeast generated from teleseisms sourced to the NE and SW (ESZ and KTSZ, respectively; Fig. 7). This effect suggests that the waves arriving from different azimuths encountered different subsurface properties in the crust beneath Uturuncu. This could be the result of a preferential orientation of fluid-filled cracks, as is commonly the case in the crust (e.g., Savage, 1999; Leidig and Zandt, 2003). This orientation agrees with the NW-SE alignment of earthquake locations (Jay et al., 2012) as well as the orientation of the Lipez-Coranzuli lineament (Lema and Choque, 1996) and the 330° upper-crustal anisotropy (Leidig and Zandt, 2003; Zandt et al., 2003). Thus, our teleseismic data are in agreement with these independently observed trends.

Teleseismic Waveform Bending

We observe that the predicted and calculated back azimuths for the teleseismic events differ (Fig. 8). This suggests that these waves interact with one or more strong refractors along the raypath. Some variation between theoretical and calculated values of back azimuth may be the result of misorientation of seismometers (Koch and Kradolfer, 1999), uncertainties in earthquake locations, errors in observation at the station, errors in data processing, or geology along the raypath—i.e., lateral velocity variations, as from dipping interfaces (Krüger and Weber, 1992; Lin and Roecker, 1996; Koch and Kradolfer, 1997). Seismic waves are refracted out of high-velocity slabs (Vidale, 1987), causing a systematic change in raypath for seismic energy that traveled through the slab. We believe that this is what gives the broad counterclockwise-clockwise trends in our azimuth study. Because the seismic energy is traveling from great distance, we don’t expect the waves to be influenced by heterogeneities in the slab, which would distort the waveforms and cause amplitude variations (Sleep, 1973; Vidale, 1987; Sekiguchi, 1992). Therefore, the large variability between the two azimuths expressed at each station is likely the result of interaction with a refractor that is shallower than the subducting slab. Two candidates are the base of crust and the bottom of the APMB; each would give similar effects to the local rock fabric mentioned above.

Gravity and Vp/Vs Tomography Data

We compare our amplitude and arrival time observations with other geophysical data to place further constraints on possible magma bodies. Gravity data, analyzed by del Potro et al. (2013), show results that correlate with the teleseismic amplitude data. The area with the largest negative gravity anomaly (∼375 mGal, the largest area of negative density contrast) is located to the southeast of the volcano (Fig. 11) and also corresponds to areas of low seismic P-wave velocity from the tomographic inversion of Kukarina et al. (2014). Likewise, the general shape of the negative gravity anomaly is reflected in the teleseismic event amplitudes for wave arrivals coming from the northeast and southwest. Station PLAR, located within this most significant negative Bouguer anomaly, consistently shows the second smallest Qp values of any station in the network and is also consistently delayed in time, suggesting that it is both experiencing high attenuation and significant reduction in velocity. This same Bouguer anomaly extends, though at a reduced amplitude, beneath stations PL07, PLBR, PLLL, and PLSS. These same stations have low seismic Qp values as well as delayed times. This suggests that the attenuation, time, and gravity anomalies are similarly sourced. This relationship has also been seen in other volcanic terrains, such as Krafla volcano in Iceland and at Yellowstone (e.g., Le Mével, 2009; Miller and Smith, 1999).

Varying the density contrast (Δρ) among end member values considering realistic melt fraction values changes the size and shape of the observed gravity anomaly (del Potro et al., 2013). The anomaly has a thick, wide appearance at Δρ = –70 kg/m3, which is along the liquidus for a granitoid composition. However, for Δρ = –400 kg/m3, which is associated with 95 vol% melt, the anomaly is tighter and less distributed (del Potro et al., 2013). Therefore, the effect on seismic amplitudes would be less pronounced in the first case but more widespread, whereas the second case would result in a spatially tight, high-amplitude anomaly. Assuming that melt fraction contributes to the anomalies, it is likely that a combination of both low and moderate melt fraction contributes to the amplitude anomaly distribution that we observe. We have calculated values of percent partial melt (Table 4) that vary between 0 and 10 (0 and 45, if discounting the effects of water). If we look at the partial melt percentages calculated for a melt with no water, values above 35% are not likely given the lack of an S-wave shadow zone (Chu et al., 2010). To account for this, the anomaly would need to have a greater thickness along the raypath for the stations showing 35% or more, or perhaps these paths are more influenced by the presence of water in the melt. Likely, it is a combination of these two factors—varying thickness and percent of water in melt. Given 15% partial melt, the maximum thickness needed would be 16 km (24 km for dehydrated melt), which is within the realm of possibility given the model of Kukarina et al. (2014; Fig. 1) and similar to the thickness of the Ward et al. (2014) model (from 4 to 25 km, or a maximum of 21 km thick) and the thickness values determined using receiver function data (McFarlin et al., 2014). However, the mean value of thickness given 15% partial melt is 10.2 km for hydrated melt and 15.3 km for dehydrated. Our data favor a moderate percent of water-rich partial melt and, therefore, a strongly attenuating body of varying thicknesses.

Comparison with the Socorro Magma Body

The results that we see for seismic waves interacting with the APMB can be compared to seismic waves interacting with another large crustal magma body—the Socorro magma body (SMB) near the Rio Grande rift in New Mexico. These two mid-crustal magma bodies have several similarities despite their very different tectonic settings, with the APMB in a zone of orogeny and crustal thickening and the SMB in a region of extensional continental rift. The land above the SMB has been uplifting at a rate of ∼2 mm/yr since 1912 (Pearse and Fialko, 2010), exhibiting shallow seismicity and swarms at 6–7 km depth (Stankova et al., 2008). Unlike several of the most recent models of the surface of the APMB, the SMB has a relatively flat upper surface (Rinehart et al., 1979; Balch et al., 1997). The SMB has been modeled to a depth of 19 km (Sanford et al., 1977), similar to the depth of the APMB at 4–25 km depth below sea level (Zandt et al., 2003; Ward et al., 2014; Fig. 1). Additionally, the SMB values for Qp of ∼30 (Schlue et al., 1996) fall within our range of Qp values but are higher than our smallest Qp values of as low as 1.8. This suggests a lower percent partial melt in the SMB than in some, but not all, parts of the APMB. The SMB, therefore, is similar to the APMB—it has been modeled to a comparable depth, has comparable physical properties, and has been experiencing uplift. These similarities may be the necessary conditions to permit the existence of long-lived mid-crustal magma bodies.


An important large-scale limitation of this study has to do with whether or not the APMB is a continuous feature (Fig. 1). This would mean that all waves from teleseismic events must travel through it, and are all slowed and attenuated. Thus all our results are based on relative measurements.

In this paper, partial melt is determined using relative travel times through the magma body. This makes the assumption that the fastest raypath was unaffected by encountering partial melt. Therefore, the values that we calculate are minima. However, using relative travel-time delays highlights the differences seen within the seismic network. Also, given the fact that the gravity model of the magma structure under Uturuncu (del Potro et al., 2013) doesn’t necessitate that magma underlies the entire area beneath the PLUTONS seismic network, this assumption may be valid. We make a similar assumption in our amplitude analysis—that the waveform with the greatest amplitude has been unattenuated and is therefore the baseline amplitude value (A0) in our Qp calculations. This could mean that we are underestimating the attenuation occurring. However, given the logarithmic nature of the calculation, even increasing A0 by a factor of 10 resulted in very little change in the lowest Qp values (a decrease of <3 in Qp values ≤10). The greatest underestimates would occur for those values of Qp that are 30 and above.

There were several other limitations to our study, three of which deal with the geometry of our network and of teleseismic earthquake distribution. Firstly, our study was limited by the aperture of the 62–97-km-diameter seismic network, whereas the feature of interest may extend beyond the network. However, the center of the inflation source is ∼3 km from the center of the network, and the Bouguer gravity anomaly shows significant variation within the magma body beneath our network. Thus, the network straddles the strong part of the anomalous zone. Secondly, the study was limited by the distribution and number of distant large subduction zone events. The threshold between clear, unfiltered events and those needing to be filtered to clearly see the teleseism differs slightly by distance from the earthquake to the network but generally is the equivalent of Mw = 5. We analyzed deep teleseisms because we wanted to limit the waveforms to those traveling through the asthenosphere and crust only once in the vicinity of the target volcano, therefore ensuring that variations seen in the waveform resulted from the crust beneath Uturuncu. Also, the epicentral distances ensure a relatively steep angle of incidence, with a range between 4° and 24° with respect to the vertical. Our study would have benefitted from large teleseismic events from similar depths and well distributed around the volcano. Because this is not the case, our analysis may miss or underrepresent structures that are not aligned with the azimuths of the seismic data.

Detailed analyses of waveforms from 14 teleseismic events show that the peak-to-peak amplitudes, zero-to-peak amplitudes, and relative travel-time residuals change across the network in a way that agrees with, and better constrains, models of the thickening of the APMB in the vicinity of Uturuncu volcano. There is also azimuthally dependent variability oriented in a NW-SE direction, suggesting significant crustal changes. We are not the first to observe varying values of Q encountered during different seismic raypaths under Uturuncu (Jay et al., 2012); however, our results are more quantitative and cover a larger area. Also, the attenuation results correlate well with a Bouguer gravity anomaly, suggesting that the source of the gravity anomaly and that of the attenuation anomalies are related.

Our results are complementary to studies of gravity, receiver functions, magnetotellurics, and Vp/Vs tomography and other studies and provide additional constraints on the magma body beneath Uturuncu. Specifically, our data suggest a magma body with <10%–20% partial melt. This magma body is 14 by 34 km in size, with long axis NW-SE. It has an average thickness of 10.2 km given a percent partial melt of 15% and 8 wt% water and is centered ∼20 km SE of Uturuncu’s summit.

Seismic attenuation, time delays, and raypath bending of teleseisms beneath Uturuncu volcano, Bolivia

Alexandra K. Farrell, Stephen R. McNutt, and Glenn Thompson

(first published on 30 March 2017, doi:10.1130/GES01354.1)

When this article was originally published, one author was inadvertently left out of the author list. Glenn Thompson has been added to the author list on the first page of the paper.

We acknowledge the following sources of data and software: National Earthquake Information Center (NEIC), Boulder Real Time Technologies (BRTT), TauP, South Florida Information Access (sofia.usgs.gov), Florida Department of Environmental Protection (dep.state.fl.us), and the Waveform suite. The paper benefitted from discussions with M. West, R. del Potro, G. Thompson, J. Braunmiller, and H. McFarlin. Funding was provided by the National Science Foundation Continental Dynamics program via a subcontract to University of South Florida from the University of Alaska Fairbanks, project #0909254.

1Supplemental Materials. A map showing the results of testing our method of determining amplitudes in a different location, additional seismograms, and additional amplitude vs. time graphs. Please visit http://doi.org/10.1130/GES01354.S1 or the full-text article on www.gsapubs.org to view the Supplemental Materials.

Supplementary data