Receiver function analyses of Uturuncu volcano, Bolivia and vicinity

The Central volcanic zone (CVZ) along the South American active margin is formed by subduction of the Nazca plate under the South American plate (e.g., de Silva, 1989; de Silva and Francis, 1991). The CVZ contains an active volcanic arc caused by this convergence and includes areas of Peru, Chile, southwestern Bolivia, and northwestern Argentina (e.g., de Silva, 1989; Allmendinger et al., 1997). Known as the Altiplano-Puna volcanic complex (APVC), this province within the CVZ is the result of a large, young “ignimbrite flare-up” that occurred less than 11 Ma in the Late Miocene (e.g., de Silva, 1989). The APVC comprises one of the most sizeable yet youngest active silicic volcanic regions with recently formed calderas on the planet (e.g., Baker and Francis, 1978; de Silva, 1989; de Silva et al., 2006b; Michelfelder et al., 2013).

The APVC, and thus, southwestern Bolivia, is underlain by a large sill-like structure known as the Altiplano-Puna magma body (APMB) (e.g., Zandt et al., 2003), which is evidenced by low seismic velocities (e.g., Chmielowski et al., 1999; Zandt et al., 2003) and low electrical resistivity (e.g., Schilling et al., 2006; Comeau et al., 2015), as well as high seismic attenuation (e.g., Haberland et al., 2003). The APMB region is unusually large for a magma body, ~50,000 km² between 21° and 24° S and 65.2° and 68.5° W. Considered to be the largest continental crustal magma reservoir in the world with a diameter of 200 km and a volume of 500,000 km³ (e.g., Ward et al., 2014; Perkins et al., 2016), it is thought to be generated by crustal magmatism, which, in turn, is caused by decompression melting from delamination and slab roll-back (e.g., Kay and Kay, 1993; Allmendinger and Gubbels, 1996; Allmendinger et al., 1997; Kay and Coira, 2009; Kay et al., 2010; Ward et al., 2014; Pritchard and Gregg, 2016). The thickening of the lithosphere, caused by the collision of the Nazca and South American plates, induced crustal shortening and delamination (Isacks, 1988; Kay and Kay; 1993). Less dense asthenosphere moved into the space between the foundering slab and the base of the continental lithosphere (Kay and Kay, 1993). This “pushed” the delaminating slab “backwards.” Currently, the APVC ignimbrites are speculated to have a common mid-crustal source identified as the APMB (e.g., de Silva, 1989; Ort et al., 1996; Schilling et al., 1997; Schmitz et al., 1997; Chmielowski et al., 1999; Zandt et al., 2003; de Silva et al., 2006a; de Silva et al., 2006b; Sparks et al., 2008; del Potro et al., 2013; Muir et al., 2014; Ward et al., 2014). Some studies, however, suggest that the APVC ignimibrites equilibrated at a depth that is shallower than the APMB (e.g., Lindsay et al., 2001; Schmitt et al., 2001).

Uturuncu volcano (22°15′S, 67°12′W), located in the Altiplano-Puna region of southwestern Bolivia (Fig. 1A), is a dacitic composite volcano of Pleistocene age, whose last eruption was ~270,000 years ago (e.g., Sparks et al., 2008). It is ~125 km east of the active frontal Andean volcanic arc (e.g., Michelfelder et al., 2014) and is in a region where the crust is estimated to be ~70 km thick (e.g., Sparks et al., 2008). At 6008 m above mean sea level, Uturuncu is the highest peak in southwestern Bolivia and well above the average regional elevation of 3800 m (e.g., Sparks et al., 2008; Jay et al., 2012) and above the average local elevation of ~5000 m (M. Pritchard, 2017, personal commun.). While the eruptions from Uturuncu have been mainly effusive in nature, the eruptive episodes have spanned ~620,000 years with a total of ~105 lava flows and domes identified and with repose intervals of 6000–8000 years (e.g., Sparks et al., 2008; Michelfelder et al., 2014).

Uturuncu has a high seismicity rate (~2–3 events hourly) when compared to other nearby volcanoes. This shallow seismicity generally has source depths of 3–4 km bmsl, or 7–11 km below the surface (e.g., Pritchard and Simons, 2004; Sparks et al., 2008; Jay et al., 2012; Hickey et al., 2013). Recent source depths calculated by West et al. (2013) show a more shallow distribution of events from ~1 km below Uturuncu’s summit to ~20 km below the surface (Kukarina et al., 2017; M. Pritchard, 2016 , personal commun.). The cause of the elevated seismicity is currently unknown but is thought to be either related to current uplift (Sparks et al., 2008) or to a shallow hydrothermal system heated by the APMB and is evidenced by active fumaroles near Uturuncu’s summit (e.g., Jay et al., 2012; Hickey et al., 2013).

Uturuncu has also been shown to be inflating, based on InSAR data, at nearly 1–2 cm/yr since 1992, over an area that has a width of nearly 70 km (e.g., Pritchard and Simons, 2002; Sparks, et al., 2008; Henderson and Pritchard, 2013); however, this inflation is decelerating (e.g., Henderson and Pritchard, 2017). The main deformation signal, according to geodetic studies, is centered ~4 km to the southwest of Uturuncu’s summit and is the only place within the boundaries of the APMB that has been shown to be actively deforming. This signal, which previous studies show to be sombrero-shaped with a central uplifting area and a larger surrounding “moat” (e.g., Fialko and Pearse, 2012; Walter and Motagh, 2014), has been attributed to the active and ongoing intrusion of magma into the APMB. The central uplifting region has a fracture ring with an approximate diameter of 51 km; this ring is causing a second, smaller “moat” around Uturuncu (Walter and Motagh, 2014). The larger “moat” has been subsiding at ~0.2 cm/yr over a 150 km diameter (e.g., Fialko and Pearse, 2012; Henderson and Pritchard, 2013; Hickey and Gottsmann, 2014). The volume change is estimated to be ~4 × 10¹⁰ m³ from 1992 to 2006 (e.g., Pritchard and Simons, 2004; Sparks et al., 2008; Hickey et al., 2013; Walter and Motagh, 2014).

The depth to the source of this inflation and surrounding subsidence is critical; it is what determines the spatial span of deformation (Hickey and Gottsmann, 2014). The outcome of this influx of material has yet to be predicted; it is currently unclear as to whether the intrusion will rise to the surface and erupt effusively as it has in the past, erupt more explosively with ignimbrite-forming eruptions, or whether it will simply accumulate and cool to form a pluton. The answers to these questions are critical to the hazard assessments for this region and can be found through a more detailed understanding of the magma plumbing system, including its depths, extents and/or shape, percent melt, and numbers of reservoirs and how they are connected.

Teleseismic receiver function analysis has proven to be a useful technique in the detection and characterization of crustal magma bodies such as the Socorro magma body (Sheetz and Schlue, 1992), the Aso caldera magma body (Abe et al., 2010), the Taupo volcanic zone (Bannister et al., 2004), and Mount Paektu (Kyong-Song et al., 2016), to name a few. Previous receiver function analyses of the APVC imaged a low-velocity zone at ~15–25 km depth from the surface with a 10%–20% reduction in seismic velocity from the surrounding crust (Yuan et al., 2000), with the possibility of an ultra–low-velocity zone (1 km/s) beneath Uturuncu volcano (e.g., Chmielowski et al., 1999; Leidig and Zandt, 2003; Zandt et al., 2003; Wölbern et al., 2009; Walter and Motagh, 2014). Ward et al. (2014), in a joint inversion of receiver functions and surface-wave dispersion, show these zones to have shear velocities of 2.9–3.2 km/s for the low-velocity zone, and a shear velocity of ~1.9 km/s for the ultra–low-velocity zone.

Coupled with the findings of Pritchard and Simons (2002) of a large deformation signal at this volcano, these previous studies sparked an intense interest in the crustal structure under Uturuncu. The PLUTONS—Probing Lazufre and Uturuncu Together: NSF (USA), NERC (UK), NSERC (Canada), Sergeotecmin (Bolivia), Observatoria San Calixto (Bolivia), Universidad Nacional de Salta (Argentina), Universidad Mayor San Andres (Bolivia), Universidad de Poto Si (Bolivia), Sernap (Bolivia), Chilean Seismological Service, Universidad de San Juan (Argentina)—project deployed 28 temporary broadband seismic stations around Uturuncu volcano from April 2010 to November 2012 (Fig. 2 and Supplemental Table S1¹). For this study, we have analyzed teleseismic receiver functions to better constrain the crustal structure beneath Uturuncu volcano and its surrounding area. These receiver functions have the potential to provide not only more accurate depths and thicknesses of the APMB but also to give a better idea of the shape of this feature, especially the topography of its upper and lower surfaces.

Twenty-eight three-component broadband seismometers (Guralp CMG-3T) were deployed around Uturuncu volcano as part of the PLUTONS project (Fig. 1B). The instruments were installed by direct burial at depths of 0.5–1.0 m to minimize surface effects. These stations were powered by 65 W solar panels and 100 A-hr sealed lead acid batteries. Data were recorded on Reftek RT130 data loggers at 100 sps and were collected during service runs every six months throughout the deployment. The deployment period for these stations began in April 2010 and ended in November 2012 (Table S1 [^{footnote 1}]). All data are also publicly available from the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC) Web site. Earthquake metadata with m_b = 5.5 and Δ (epicentral distance) = 28°–102° for direct P waveforms and Δ = 86°–173° for PP waveforms, were taken from the National Earthquake Information Center (NEIC) Preliminary Determination of Epicenter (PDE) database, and the waveforms extracted from the PLUTONS database (Fig. 1C and Supplemental Table S2²).

Receiver functions are an expression of the seismic impedance (the product of density and velocity) contrasts below a seismic station, providing insight to any changes in the velocity and/or density of, for example, the upper mantle and crust (e.g., Langston, 1979; Rychert et al., 2007; Frassetto and Thybo, 2013). Teleseismic events are used so that the incoming wave is near vertical as it approaches the station. The vertical component of a three-component seismogram is deconvolved from the horizontal components (N and E components rotated to radial and tangential) to isolate the response from the local crustal structure (Langston, 1979). These receiver functions contain phase arrivals that are related to velocity contrasts below the seismic stations. In this paper, we use the common conversion point (CCP) stacking method of Dueker and Sheehan (1997) to interpret the locations of these interfaces.

After the receiver functions were calculated, we manually reviewed all resulting receiver function traces for quality control. Receiver functions are removed from further processing if they are noisy or deviate substantially from other receiver functions at the same back-azimuth and/or distance range. There are many factors that affect the quality of the final receiver function, including small P arrivals and strong lateral velocity gradients. The small P arrivals lead to small conversions. Strong lateral velocity gradients can lead to off-azimuth arrivals. More importantly, lateral velocity gradient can violate an inherent assumption in receiver function analysis—that conversions come from near horizontal boundaries (Langston, 1979).

During this review, we noticed that certain stations (specifically PL07, PLHS, PLMK, and PLTP; Fig. 1B) produced consistently poor-quality receiver functions from some back azimuths. These receiver functions tended to have larger amplitude tangential components when compared to the radial components, and many had radial components with a reversed polarity for the direct P arrival. These receiver functions were removed, and the cause of this will be the motivation for other studies, but are most likely due to off-azimuth arrivals from strong lateral velocity gradients (Farrell et al., 2017).

After review, a total of 799 P and PP receiver functions from 192 separate events (shown in Fig. 1C and Table S2 [^{footnote 2}]), were inputs into the CCP-stacking code written by Dueker and Sheehan (1997) and modified by Gilbert et al. (2003) and Frassetto et al. (2010). CCP stacking is a process in which receiver function arrivals with the same conversion point are binned and summed together after corrections for move-out have been applied. This is done so that noise is stacked out (destructive interference) and the major reflections are emphasized (constructive interference) (Neuendorf et al., 2005) and to account for spatial variation in the impedance contrasts.

The CCP program does not invert for seismic velocities and thus requires a velocity model in order to place arrivals into appropriate spatial bins. In this study, the 3-D velocity model from Ward et al. (2014) was used in conjunction with the S-wave IASP91 model as modified by Ryan et al. (2016). As part of the inversion process for S-waves from surface waves, a forward model is created that results in models for S-waves, P-waves, and density, with large uncertainties for the P-waves and densities (K.M. Ward, 2017, personal commun.). Ryan et al. (2016) improved upon the P-wave velocity model by multiplying the S-wave velocities by a V_p/V_s = 1.75, taken from Wadati diagrams of local events (J. Ryan, 2017, personal commun.). Variations in V_p/V_s can lead to errors in velocities and depths; these errors, however, tend to be small (J. Ryan, 2017, personal commun.). A slightly more recent velocity model, from Shen et al. (2017), was not available at the time of analyses; however, this velocity model is very close to the Ward et al. (2014) model used in this study. In addition, the Ward et al. (2014) velocity model is three-dimensional and will thus likely provide more accurate depths than the Shen et al. (2017) model, which is a one-dimensional model. The Ward et al. (2014) velocity model was used at depths starting at the surface down to a depth of 50 km bmsl, while the modified S-wave IASP91 model is used for depths from 50 km to 150 km bmsl. The model space was parameterized into bins with a depth increment of 1.0 km from the surface to 150 km depth and horizontal bin sizes of 10 × 10 km and 15 × 15 km.

Receiver functions contain Ps conversions from major boundaries plus multiple reflections off each boundary. The weighting filter in the CCP program improves the coherence of the features by stacking the amplitudes of the first and second multiples from the strongest boundary, which happens to be the top of the APMB, of a receiver function trace into the first arrival. This parameter can range between 0.0 and 1.0. We determined, empirically, that a value of 0.25 was the value that produced the best coherence; values higher than this can create artifacts, such as spurious arrivals (examples of various weighting filters are shown in Fig. S4 of the Supplemental Figures³).

The sharing coefficient in the CCP program allows sharing of traces between bins and can have values between 1.0 and 2.0, where 1.0 is minimal sharing and 2.0 is maximum sharing. This coefficient can be thought of as a “buffer”; if a trace falls within a “buffer-zone,” then it can be shared with the nearest bin. There is a trade-off between bin size and the sharing parameter—if smaller bins are used and the sharing between them is increased, the result is a connectivity of the features for which we are looking. In other words, it smears the features together causing loss of definition (decreases the resolution). For example, with a 10 km bin radius, a sharing coefficient of 1.25 gathers all the traces in a 12.5 km radius (generally four to five traces are added along the edge of the bin). Neighboring bins share these traces along their edges, and the model is smoothed between bins to reduce any sharp distinctions as bin edges are crossed. Typically, a minimum number of traces for each bin are required to help eliminate potentially noisy bins with too few receiver functions. We set this number to be five traces for each bin for this study. We then ascertained, empirically, that a sharing coefficient value of 1.25 was ideal for both the 10 × 10 and 15 × 15 km bins, because it allows for adequate sharing between bins while also providing a low amount of smearing.

Figure 2 shows a map view of the study area with the center of each bin marked by the grid of solid blue dots. The distribution of the piercing points for each incident ray at 60 km depth bmsl (small red dots) gives a qualitative view of the receiver function coverage throughout the region. Letters with arrows point to a particular cross section and correspond to the letters in Figures 3A–3K for EW cross sections and the letters in Figures 4A–4K for the NS cross sections. These cross sections give a detailed view of the lateral and vertical changes in the earth structure under Uturuncu and the surrounding area.

Figures 3A–3K (EW cross sections) and Figures 4A–4K (the NS cross sections) show the results of the CCP-stacking program with a Gaussian width of 2.5 and bin sizes of 10 × 10 km; there are several very prominent features visible. The largest and most consistent feature in the cross sections is a strong negative trough shown in blue at ~6–11 km bmsl (Figs. 3–5). A negative pulse is produced by conversions at a low- to high-velocity boundary for the incident wave. We interpret this feature as the top of the APMB (highlighted by the top gray line). The northeast section of the APMB is deeper, with an average depth of approximately ~10–11 km bmsl.

Directly below the negative pulse (Figs. 3A–3K and 4A–4K) is a smaller but still very prominent positive pulse, shaded a bright red. This is interpreted to be the bottom of the APMB. This is also marked by a gray line (dashed where the depth of the pulse is not well constrained), to highlight that this entire structure between gray lines is the APMB. This feature also exhibits topography (Figs. 3G and 4E) and is shallowest ~13–15 km depth bmsl in the central region including Uturuncu volcano and as deep as ~22 km in the eastern part of the region (Fig. 6). The thickness of the APMB can be calculated from the depths of the two surfaces (Fig. 7). The APMB is thinnest in the central region, ~6–8 km thick, which includes Uturuncu volcano and a north-south swath to the west of the volcano. The area is bordered to the east and west by thicker regions up to 16 km thick.

Another prominent positive (red) feature is at ~60 km depth bmsl (Figs. 3, 4, and 8); we believe this feature to be the Mohorovičić discontinuity (Moho), and we have highlighted it with a green line (dashed where depth of the feature is unclear). A large degree of upwarping in this feature is apparent under Uturuncu. The depth to the Moho varies more than 20 km (40 km to >60 km depth) and is shallowest in the circular region centered just south of Uturuncu (and includes the volcano). The CCP profiles (Figs. 3 and 4) show sharp offsets at the edge of the feature in many places.

The rest of the smaller, less prominent receiver function features are likely reflections and multiples of the larger attributes. A feature seen in only a few of the cross sections (Figs. 3 and 4), near the bottom at ~140 km depth bmsl, is what we interpret to be the top of the subducting slab of the Nazca plate and is highlighted with a yellow line. A cursory examination of the PDE for our study region from January 1, 2000 to July 17, 2017 gives 249 events with the majority of hypocenters between 120 and 220 km depth, with the depth increasing eastward. The feature that we see at ~140 km depth would appear to fall near the estimated depth to the subducting slab in this region.

Results of the CCP-stacking program with values of a = 2.5 and 15 × 15 km bin sizes can be found in the Supplemental Materials (see ^{footnotes 1} and ²). Depths are the same for all major features; the only real differences between the two bin sizes are variations in complexity of the receiver function stacks. Stacks with a = 2.5 and 15 × 15 km bins are less complex than those shown in the main body of this paper due to the decreased resolution of the larger bin size compared to the 10 × 10 km bin size.

Early studies of the Altiplano-Puna region—namely, those by Wigger et al. (1994) and Chmielowski et al. (1999), in which the LVZ was named “the APMB”; Beck and Zandt (2002); Zandt et al. (2003); and Leidig and Zandt (2003) identified the existence of a low-velocity zone (APMB) and determined the size and depth of this feature. In these early studies, the top of the APMB was identified as a large negative polarity arrival on receiver functions at 2–2.5 seconds, except for Wigger et al. (1994), which found a LVZ between 10 and 20 km below the surface in a seismic-refraction study. The APMB was modeled as a thin (~1-km-thick), low-velocity layer with the top boundary at ~15 km depth bmsl. Beck and Zandt (2002) show a negative polarity arrival at ~14–20 km bmsl across the entire width of the Altiplano. Leidig and Zandt (2003) estimated the thickness of this feature to be 1–6 km. Wölbern et al. (2009), in a regional migrated receiver function study, inferred a segmented LVZ under both the Altiplano and the Puna regions in which the sections are vertically separated and have steep boundaries coincident with surface fault systems. Depths determined by Wölbern et al. (2009) include a top boundary for the APMB of ~15 km below the surface (~10–11 km bmsl) in the Altiplano with a thickness of ~20 km, and ~15 km depth in the Puna, with an unresolved thickness for this region. In addition, the aforementioned study shows the boundaries for the APMB have significant topography; an average depth of ~15 km is given, but the authors state that the low-velocity zone in this study is segmented and shows 8–10 km of vertical offset in some areas, indicative of faults dipping toward the west. More recently, Ward et al. (2014), using a joint inversion of surface-wave dispersion (Ward et al., 2013) and receiver functions, found that the top of the APMB is as shallow as 9 km below the surface with an average depth of ~15 km, and they calculated a thickness of no less than 6 km, with an average thickness of 11 km. The main difference between Ward et al.’s (2014) shallower and thicker APMB and earlier studies showing a deeper and thinner APMB is that the surface-wave dispersion required lower velocities in the upper crust.

In this study, we used the 3-D velocity model of Ward et al. (2014), in addition to the S-wave IASP91 model as modified by Ryan et al. (2016) in the processing of the CCP stacks; therefore, our depths for the top and the bottom of the APMB generally match those found by Ward et al. (2014). Depths for the top and bottom of the APMB were determined by hand; we picked the peak of large negative pulses for the top of the APMB. The bottom of the APMB was more difficult due to the diffuse nature of the positive pulses in this depth range. The depth to the top of the APMB as reported in this paper (6–12 km bmsl; Fig. 5) is shallower than most of the earlier studies (e.g., Chmielowski et al., 1999; Leidig and Zandt, 2003; Zandt et al., 2003; Wölbern et al., 2009) but matches well with those found by Ward et al. (2014). Depths to the bottom of the APMB (Fig. 6), however, are shallower than those found by Ward et al. (2014) and other previous studies. The improvement we show, as compared to the Ward et al. (2014) study, is that of more accurate local depths, versus the regional depths that have been smoothed due to the nature of a regional study and the method of surface-wave dispersion. We have calculated the thickness to range from 5 km to 15 km (9 km average thickness; Fig. 7). Farrell et al. (2017), using P-wave velocities and partial melt values to empirically calculate the average thickness of the magma body, report this thickness to be 10.2 km for the APMB in the vicinity of Uturuncu, given a partial melt percentage of 15% with 8 wt% water.

The amplitudes of CCP-stacked receiver functions are related to the thickness (or sharpness) of boundaries or transition zones as well as the impedance contrast between the layers, as shown by Dr. Charles J. Ammon on the Web site “An Overview of Receiver Function Analysis” (http://eqseis.geosc.psu.edu/~cammon/HTML/RftnDocs/seq01.html; Purevsuren, 2014). These amplitudes are influenced by lateral velocity changes—offsets in the topography of a boundary or transition zone can cause scattering of converted phases. Eagar et al. (2011) demonstrated that mapping the amplitudes of receiver functions can provide a better visualization of complicated subsurface structures. The ability to view the patterns of reduced amplitudes can lead to an improved understanding of where dipping or deformed structures are located or the shapes and extents of bodies causing the scattering or attenuation of converted seismic phases. We measured the amplitudes of the top and bottom of the APMB in the CCP stacks and calculated the ratio (amplitude of the top feature divided by the amplitude of the bottom feature).

Amplitude ratios for the APMB are plotted on Figure 9. A low ratio is indicative of similar impedances for both the top and bottom of the APMB. A high ratio is indicative of large variations in impedance contrast between the top and bottom of the APMB—in this case, a strong impedance contrast for the top of the APMB and a reduced amplitude or a very broad peak for the bottom of the APMB. In general, most of the region shows amplitude ratios of ~2 or 3. In Figure 9, the high-ratio region, near station PLRV, is due to large-amplitude, impulsive peaks for the top of the APMB and smaller amplitude but very broad peaks for the bottom of the APMB (Figs. 3F, 3G, 4A, and 4B). By contrast, the high ratios in the northwest to central region, under stations PLJR, PLSM, PLMN, and between PLQU and PLLA, are due to impulsive peaks for the top of the APMB, and slightly less impulsive (more diffuse), low-amplitude peaks for the bottom of the APMB (Figs. 3F, 3G, and 4C–4E). The high-amplitude ratios in the southernmost area of the network, between stations PLTT, PLMD, and PLSP, are due to large amplitudes for the top boundary and smaller broad peaks for the bottom boundary. In general, high ratios indicate that the bottom of the APMB is more poorly defined with a smaller impedance than the top boundary. In addition to differences in how the boundary is being defined (i.e., defining the boundary using an interface versus defining the boundary using velocity contours as in Ward et al., 2014), this diffuse boundary may be why the estimates for the bottom of the APMB do not agree as well with previous work, whereas the estimates for the top tend to correlate well with other studies.

Previous regional studies show that the Moho in this region of the Andes can have variable depths (e.g., Wigger et al., 1994; Beck and Zandt, 2002; Yuan et al., 2002; Wölbern et al., 2009; Ryan et al., 2016). Wigger et al. (1994) reported that data quality issues and attenuated phases made identification of the deeper structures, such as the Moho feature, difficult. The attenuated phases were assumed to represent zones of partial melt. When good data were available, Wigger et al. (1994) reported a possible seismic crust and/or mantle boundary at ~60 km depth below the surface in the Western Cordillera, at ~70 km depth below the surface in the Eastern Altiplano with strong lateral variations, which decrease in depth eastward to 40 km below the surface in the Subandean belt. Yuan et al. (2002) found depths as shallow as ~51 km and as deep as 80 km in the Altiplano-Puna region using teleseismic receiver functions, but they also report that for multiple smaller areas in this region, a depth to the Moho was unresolvable due to the large conversions and multiples from the APMB. Beck and Zandt (2002), using receiver functions, found that the Moho in the Central Andes varies between 80 and 32 km depth, with the depth decreasing toward the east. Specifically, in the Altiplano region, the depth of the Moho in the Beck and Zandt (2002) study is ~64 km. Wölbern et al. (2009) used migrated receiver functions that revealed a Moho with a weak signal that shallows eastward from ~80 km below the surface to ~45 km in the Altiplano, with a slight upwarping between 67°W and 69°W. Ryan et al. (2016) found similar results to Beck and Zandt (2002), providing a Moho depth of nearly 70 km bmsl in the Altiplano, but they state that under the transition from between the western plateau and the arc, the depth region of the Moho is a complex with many layers, some of which are likely multiples from the LVZ.

Moho depths as determined from the CCP stacks from this study are plotted in Figure 8. Depths were picked using the same method as for the top and bottom of the APMB, at the peak of the positive pulse; however, this was difficult because the pulses in this depth range are not only diffuse, but some traces show multiple pulses at Moho depths. Nevertheless, the Moho arrivals are fairly strong and consistent throughout most of the study area. Depths range from as deep as 65 km to as shallow as 40 km. The Moho is shallowest in a roughly circular area that extends from Uturuncu volcano (near station PLCM) to the south. This region of shallow Moho is ~50 km in diameter. The transition from shallow Moho (40–45km) to deeper Moho (55–60 km) is fairly sharp (see Figs. 3 and 4), particularly on the east and south sides where the depth increase is over 5–10 km horizontal distance (versus over 20 km lateral distance on the west and north sides) and may possibly even be fault bounded. This region of shallow Moho coincides with a region of faster velocities as shown by Ward et al. (2014).

A regional travel-time tomography study by West et al. (2013), using the same stations as this study, found a high V_p/V_s ratio (>1.9), consistent with a low S-wave velocity anomaly, which covers an area virtually the same as our shallow Moho. The low V_s anomaly is in a well-resolved area of the tomography model. These features are roughly circular (in map view) with ~50 km diameter and are most prominent at ~15–30 km bmsl between 22.2° and 22.7° S and 66.9° and 67.5° W. These anomalies occur above the Moho and encompass the APMB. S-wave velocities tend to be the most responsive to liquid phases, specifically the percentage and interconnectedness of partial melts, thus giving a low S-wave velocity for areas with melt present. The anomalies in the aforementioned tomography study are interpreted to be either a vertically elongated buoyant magma chamber or a tract of volatiles and/or partial melts that are ascending.

We were unable to resolve the Moho depth under the stations PLJR and PLSM due to the lack of distinct converted arrivals. These two stations are directly above the area of thinnest APMB and highest amplitude ratio, indicating a weak conversion at the bottom of the APMB (see Figs. 8 and 9). To clarify, the amplitude ratio quantifies the difference in sharpness of the boundaries at the top and bottom of the APMB. It is interesting to note that the Moho in this region does not produce a sharp arrival and thus is more gradational (as is the bottom of the APMB).

Cahill and Isacks (1992), based on earthquake locations, show that the subducting Nazca plate ranges in depth from 125 km to 200 km (west to east) across the southern tip of Bolivia. The recent Slab1.0 (Hayes et al., 2012), a new global 3-D model of subduction zones, shows that the subducting slab in this region of the Andes is between ~100–200 km depth. Although our model space ends at 150 km below mean sea level, it does appear that the subducting slab is visible in some of the cross sections (Figs. 3 and 4) as a positive (shaded red) pulse from ~135 km to ~150 km depth. This feature is particularly visible on the western and southern regions of the study area and is not visible in the eastern regions where we expect the subducting plate to be more than 150 km deep.

CCP stacking of teleseismic receiver functions help constrain the crustal structure of the complex area under Uturuncu volcano. Large impedance contrasts in the receiver function stacks allow for the identification of major boundaries and also the calculation of the depths of these boundaries below the study area. Mapping the depths of these features provides a determination of their extent and shape, potentially aiding in future hazard assessments for Uturuncu volcano.

The cross sections shown allow for a better understanding of the complexity of the magma system under Uturuncu volcano. A major feature that is seen in all cross sections (and in fact extends past the limits of the network in all directions) is the APMB. Depths to the top of the Altiplano-Puna magma body within the study area are as shallow as ~6 km bmsl and as deep as ~22 km, with an average of ~8–9 km depth (Fig. 5). Depths to the bottom of the magma body are as shallow as ~13 km and as deep as ~22 km bmsl, giving a thickness ranging from 6 to 15 km (Figs. 6 and 7). The APMB appears to be thinnest in a roughly north- and/or south-elongated region near the center of the study area including Uturuncu volcano. The highest relative amplitude ratios tend to fall in this region of thinner APMB indicating a weak lower boundary. Mohorovičić discontinuity depths are somewhat bimodal with a shallow central region (depths of 40–45 km) surrounded by a deeper Moho at 55–60 km below the surface. The roughly circular region with shallow Moho corresponds to a low S-wave velocity (and high V_p/V_s ratio) region in the crust from tomographic images (West et al., 2013).

Receiver functions are extremely useful in determining the depths of sharp horizontal boundaries in the crust such as the Moho and, in this study, the top and bottom of the APMB. While the depths of these boundaries are fairly well determined, there can be alternative interpretation as to the nature and process that produce each boundary. In this case, explanations for the shallow “Moho” feature beneath Uturuncu volcano such as magmatic underplating, piles of cumulates, or crustal delamination cannot be ruled out, but require additional information to distinguish between the possible models.

The authors would like to acknowledge National Earthquake Information Center (NEIC), Boulder Real Time Technologies (BRTT) and the Global Seismology and Tectonics (GSAT) group at University of Arizona for software and data provided, as well as the two anonymous reviewers for their excellent feedback and critiques. Discussions with Jochen Braunmiller, Michael West, Matthew Pritchard, Alexandra Farrell, Susan Beck, Ophelia George, Rocco Malservisi, Danielle Molisee, and Mel Rodgers aided in the writing of this paper. The National Science Foundation Continental Dynamics Program was the source of funding, with a subcontract from University of Alaska Fairbanks to University of South Florida, with the project #0909254.