Using receiver functions, Rayleigh wave phase velocity dispersion determined from ambient noise and teleseismic earthquakes, and Rayleigh wave horizontal to vertical ground motion amplitude ratios from earthquakes observed across the PLUTONS seismic array, we construct a one-dimensional (1-D) S-wave velocity (Vs) seismic model with uncertainties for Uturuncu volcano, Bolivia, located in the central Andes and overlying the eastward-subducting Nazca plate. We find a fast upper crustal lid placed upon a low-velocity zone (LVZ) in the mid-crust. By incorporating all three types of measurements with complimentary sensitivity, we also explore the average density and Vp/Vs (ratio of P-wave to S-wave velocity) structures beneath the young silicic volcanic field. We observe slightly higher Vp/Vs and a decrease in density near the LVZ, which implies a dacitic source of the partially molten magma body. We exploit the impact of the 1-D model on full moment tensor inversion for the two largest local earthquakes recorded (both magnitude ∼3), demonstrating that the 1-D model influences the waveform fits and the estimated source type for the full moment tensor. Our 1-D model can serve as a robust starting point for future efforts to determine a three-dimensional velocity model for Uturuncu volcano.

The Altiplano of the central Andes is a broad plateau with substantial relief (elevations 3.5–4.7 km) that overlies a 70-km-thick crust (Isacks, 1988). Located in the Altiplano area, Cerro Uturuncu is a stratovolcano within a regional cluster of young calderas and ignimbrite sheets termed the Altiplano-Puna volcanic complex (APVC) (de Silva, 1989). Beneath the APVC, a regional-scale mid-crustal magma reservoir termed the Altiplano-Puna magma body (APMB) was identified at depth of ∼15 km below sea level (19 km below the surface) by receiver functions (Chmielowski et al., 1999) and many other geophysical studies. These studies revealed several geophysical anomalies attributed to this magma body, which include: (1) very low seismic velocities and high seismic attenuation (Chmielowski et al., 1999; Zandt et al., 2003); (2) low electrical resistivities (Haberland et al., 2003; Comeau et al., 2015), and (3) low densities that produce negative gravity anomalies (del Potro et al., 2013). Regarding the locations, all of these anomalies revealed by these studies emerge at a depth of ∼20 km below the surface. In addition, kinematic studies using geodetic satellite observations have identified ongoing surface uplift and peripheral subsidence centered at Uturuncu volcano (Pritchard and Simons, 2002; Fialko and Pearse, 2012).

Since 2009, Uturuncu has been the focus of a multi-disciplinary project known as PLUTONS. The goal of the project is to study how magma accumulates and erupts in areas of active crustal formation, mid-crustal melt and intrusion, and volcanism. The scientific effort is rooted in data from seismology, magnetotellurics, GPS, InSAR, and gravity. The seismology component of PLUTONS comprises an array of up to 31 broadband seismic stations deployed from 2010 to 2012 (West and Christensen, 2010), the locations of which are presented in Fig. 1. Using this array, several seismic studies have begun to illuminate the structure and processes beneath the APVC and particularly beneath Uturuncu (Jay et al., 2011; Ward et al., 2013). Additionally, a more recent study using receiver functions and surface waves observed by ambient noise estimated that the APMB is ∼200 km wide, 11 km thick, at depths 14–25 km, and ∼500,000 km3 in volume (Ward et al., 2014). However, none of the studies utilized a wide spectra of seismic observables (i.e., integration of receiver function waveforms, surface wave dispersion from both ambient noise and earthquakes, and other seismic observables) to provide a comprehensive one-dimensional (1-D) seismic model with rigorous error analysis.

One-dimensional velocity models are widely used at seismic observatories for estimating source parameters including origin times, hypocenters, and source mechanisms (e.g., Ekström et al., 2012; Sipkin, 1986). Using a 1-D model, Alvizuri and Tape (2016) estimated full moment tensors for 63 events, finding a predominance of events with significant positive isotropic components. In that study, a homogeneous seismic model was chosen over the more realistic model from this study. However, the systematical investigation of the impact on the estimate of the source parameters using a more realistic 1-D model was not presented in detail in that paper, and we intend to document it here in this study.

In this paper we document a 1-D reference model beneath Uturuncu constructed using three types of seismic observables—Rayleigh wave horizontal to vertical ground motion amplitude ratio (H/V ratio) observed from teleseismic earthquakes, Rayleigh wave phase velocity dispersion from both ambient noise and earthquakes, and receiver function waveforms—through a Monte Carlo inversion algorithm. The incorporation of Rayleigh wave H/V ratio, in particular, improves the resolution of the P-wave to S-wave velocity ratio (Vp/Vs) and density structure in the upper and mid-crust (Lin et al., 2012, 2014), and the joint inversion of all three types of data within the Bayesian Monte Carlo framework provides rigorous error estimates of the inversion result by presenting the posterior distributions of model attributes (Shen et al., 2013a), which improves the usefulness of the resulting 1-D model for other applications. We describe the observation of all three types of data and the inversion algorithm in the next section. To address the effect of the more realistic 1-D model on the source estimation, we use the resulting 1-D model from the inversion to derive synthetic seismograms which are used to estimate seismic full moment tensors, and the impact of the model on estimating the moment tensors is quantitatively assessed. The discussion of the resulting 1-D model and its impact to moment tensor estimates are presented in the Discussion section. Finally, we summarize our conclusion as well as caveats of this study in the Conclusion section.

Data Processing

In this study, we incorporate three types of seismic observables: (1) P-wave receiver function (PRF) waveform; (2) Rayleigh wave phase velocity dispersion from ambient noise and earthquakes; and (3) Rayleigh wave H/V ratio measured from earthquakes. The data processing for each data set has been well documented by previous studies in other areas (Lin et al., 2008, 2009, 2012; Lin and Ritzwoller, 2011; Shen et al., 2013a, 2013b; Deng et al., 2015) and is only briefly summarized here.

For PRFs, we collect over 750 raw PRF waveforms across the PLUTONS array from a previous study (Ward et al., 2014). After normalizing the amplitudes and P-to-S moveout of these PRF waveforms to a slowness of 0.06 s/km, a comprehensive quality control is performed following the method described by Deng et al. (2015). This quality control ensures that the RFs with similar back-azimuthal angles that are coherent with each other are kept and the ones that are not coherent with others are discarded. Finally, 248 RFs pass the quality control and are stacked to obtain the average PRF. The standard deviation of the stack is taken as the error of the PRF data. Figure 2A presents the average PRF with uncertainty as a gray corridor.

For Rayleigh wave phase velocity dispersion, two separate dispersions with uncertainties are computed. First, ambient noise cross-correlations between all PLUTONS stations are calculated, and the traditional frequency-time analysis (FTAN; Levshin and Ritzwoller, 2001; Lin et al., 2008) is used to determine phase velocities between 6 and 20 s period. For each period, all phase velocity measurements satisfying the selection criteria (signal-to-noise ratio [SNR] > 8 and distance larger than three wavelengths) are averaged, and the standard deviation of the mean is used to evaluate the uncertainty. Second, teleseismic Rayleigh waves observed across the PLUTONS array are used to determine phase velocities across the array between 24 and 100 s period using eikonal tomography (Lin et al., 2009; Lin and Ritzwoller, 2011). Here, all available earthquakes with magnitude >5.0 are included. For each period, all measurements with SNR > 8 are used to calculate the averaged phase velocity, and the standard deviation of mean is used to evaluate the uncertainty. The two curves from ambient noise cross-correlation and earthquake measurements are then combined into one average dispersion, which is shown as black error bars in Figure 2B.

For Rayleigh wave H/V ratios, we use the same teleseismic events we used in our phase velocity analysis across the PLUTONS array. To mitigate the potential effect of off-great-circle propagation on the horizontal amplitude measurements, for each event and period, we first determine the wave propagation direction based on the phase front tracking method of the eikonal tomography (Lin et al., 2009; Lin and Ritzwoller, 2011). The ratio between amplitude in the direction of wave propagation and that of the vertical component is then calculated for each station. All available H/V ratio measurements across the array are then averaged, and the standard deviation of the mean is used to evaluate the uncertainty. Figure 2C summarizes the Rayleigh wave H/V ratios and uncertainties observed between 24 and 100 s period, shown as black error bars.

Model Parameterization and Model Space

Because both Rayleigh wave and PRF measurements are mainly sensitive to vertically polarized shear wave velocity (Vsv) rather than horizontally polarized shear wave velocity (Vsh), radial anisotropy is ignored in this study. Here we assume that Vsh = Vsv = Vs to simplify the 1-D inversion process. The 1-D model is parameterized by four layers: (1) a sedimentary layer described by a linear velocity gradient; (2) a fast uppermost crust lid described by a single velocity layer; (3) a 50-km-thick smooth transitional layer from the base of the upper crustal lid to the Moho with Vs described by six cubic B-splines; and (4) an uppermost mantle layer extending to 200 km depth with Vs described by five cubic B-splines. In each layer, while Vs can vary with depth, density and Vp/Vs are set as constant values. For layers 1–3, we allow density and Vp/Vs to change in the inversion, and for the mantle layer (layer 4), we fix density as 3.38 g/cm3 and Vp/Vs as 1.75. Note that we assign different Vp/Vs and densities to the upper 10 km and bottom 40 km of layer 3 to investigate potential Vp/Vs and density variations in the layer. The model space is defined by the perturbations of free parameters relative to the reference values presented in Table 1. In total, 24 parameters are allowed to change in the Monte Carlo sampling.

To reduce the model space, we introduce some prior constraints for the Monte Carlo inversion. First, we force the velocity increase in the sedimentary layer. Second, we set the first B-spline coefficient of layer 3 equal to the Vs of layer 2, and set the last B-spline coefficient of layer 3 equal to the first B-spline coefficient of layer 4. This constraint reduces our free parameters from 24 to 22. During the inversion, we use the attenuation (Q) model from a 1-D reference model AK135 (Montagner and Kennett, 1996) to correct the physical dispersion effect, and our model is normalized to 1 s.

Monte Carlo Inversion and the Resulting 1-D Model

In order to invert the seismic observables for a 1-D Vs model, a two-step Markov chain Monte Carlo (MCMC) inversion is implemented. In the first step, we follow Shen et al. (2013a) and Shen and Ritzwoller (2016) to perform a joint MCMC inversion and obtain an ensemble of 1-D models that fit the data equally well. Secondly, for each model mi in the resulting model ensemble, we compute a revised misfit functional:
in which χ(mi) is the square root of reduced χ2 misfit to all data sets for model mi, and the penalty functional Pi is defined as:
in which s is depth, D is the depth of bottom of layer 3, forumla is the reference Vp/Vs value that we set up for the crust is set to 2 for sedimentary layer (layer 1) and 1.75 for layers 2 and 3. Essentially, we are seeking the model that fits our data and also has minimum perturbation of density and Vp/Vs to our reference model.

We choose the model with minimum R(mi) as our final 1-D model, which is shown as solid black profile in Figures 3A–3B. The resulting Vp/Vs and density profiles for this model are shown as blue and red lines, respectively, in Figures 3C–3D. The fit to the data from the resulting model is shown in Figures 2A–2C. Overall, the model can fit all three measurements fairly well (square root of reduced χ2 for RF: 0.6; for surface wave dispersion: 2.22; for Rayleigh wave H/V: 1.1). The only part that is not fit well is the short-period (6–12 s) Rayleigh wave phase velocity measurements, which is likely due to the near-surface complexity. By inverting all three measurements together, we also investigate the possibility of resolving Vp/Vs and density ratio (Figs. 3C–3D). We discuss the implication of the 1-D model to structure beneath Uturuncu in the section below.

Uppermost Crustal Lid and Mid-Crustal Low-Velocity Zone

As shown in Figure 3A, three abnormal features in the Vs structure are observed: (1) a fast Vs lid (∼3.1 km/s) in the uppermost crust between 2 and 10 km depth; (2) a low-velocity zone (LVZ) beneath, between 12 and 22 km, with a minimum Vs of ∼2.1 km/s at ∼17 km depth; and (3) the absence of a Moho discontinuity beneath the LVZ. All three features are generally consistent with the previous three-dimensional seismic images on a larger scale constructed by Ward et al. (2014).

To further investigate the uncertainties of the 1-D model, we perform an additional Monte Carlo search around our final 1-D model. After the search, >3000 different models with similar misfits compared to the final model are obtained, and their full extent in Vs model space is shown as the gray corridor in Figures 3A–3B. From this ensemble of models, posterior distributions emerge. Figure 4A presents the marginal posterior distributions of Vs of the crustal lid and of the LVZ at 17 km depth. The mean and standard deviation of the posterior distributions of Vs (2.22 ± 0.08 km/s for the lid and 3.15 ± 0.11 km/s for the LVZ) suggest that the two distributions are well separated (by ∼10 σ). The extremely low Vs in the middle crust is likely related to the widespread magma body beneath the APVC (Chmielowski et al., 1999; Zandt et al., 2003; Ward et al., 2014).

Although the posterior marginal distributions for Vp/Vs and density are not completely separated as in the case of Vs (Figs. 4B–4C), substantial preference from the data can still be observed: the average Vp/Vs is ∼1.65 ± 0.09 in the crustal lid and ∼1.79 ± 0.12 in the LVZ, and the average densities are ∼2.59 ± 0.14 g/cm3 and 2.47 ± 0.14 g/cm3, respectively. We conclude that the seismic data prefer a slightly higher Vp/Vs ratio and low density in the mid-crust corresponding to the LVZ (Fig. 3C). The relatively high Vp/Vs may be due to the existence of partial melt within the LVZ. Ward et al. (2014) proposed a partial melt of as much as 25% beneath the APVC. The low density of the LVZ (mean of the posterior distribution < 2.5 g/cm3) is likely also related to the hot and partially molten magma in the middle crust consistent with a dacitic composition (Sparks et al., 2008). While low velocity and high Vp/Vs ratio in the uppermost crustal layer can reflect loose materials near the surface, the average to slightly higher density (∼2.68 g/cm3) of the layer may reflect the heavier-than-normal andesite, but the σ of the posterior distribution is too large to make further conclusion (Fig. 3C). Extra caution should be taken, however, when interpreting features related to the top layer, as our 1-D model does not fully explain the complexity of our phase velocity measurements at short period (Fig. 2B).

The absence of a clear Moho discontinuity comes naturally from the prior constraint that the bottom of layer 3 and top of layer 4 are continuous, and the fit to data supports that a sharp Moho is not required in this average 1-D model. However, we cannot rule out the possibility that a regional Moho may exist beneath part of the study area because the stacking of all receiver functions may have obscured the Moho conversion in some PRFs.

Impact on Full Moment Tensor Inversion

Alvizuri and Tape (2016) performed full moment tensor inversions for 63 events at Uturuncu volcano using both the 1-D model derived from this study and a homogeneous half-space model. The events that they examined are relatively small (mostly Mw < 1.5), and their waveforms have high frequency content (2–10 Hz). Their tests show that the homogeneous half-space model provided more reliable moment tensor solutions.

Neither body waves nor high-frequency (2–10 Hz) waveforms were used in constructing the 1-D model in this study. Therefore it is perhaps not too surprising that our 1-D model does not provide fit improvement. We would expect, however, that the 1-D model would provide improvement to moment tensor solutions and waveform fits when considering surface waves at slightly longer period. Among the 63 events studied by Alvizuri and Tape (2016), two events were large enough to generate observable surface waves. In this study we use these waveforms, filtered between periods 2-4 seconds, to invert for the moment tensor. Also as in Alvizuri and Tape (2016), we utilize first-motion polarity measurements (though not waveform differences) to help constrain the solutions.

Figures 5 and 6 compare waveform fits for one of the two events (“event A”) for moment tensors derived using a 1-D model and a homogeneous model. For this event, the 1-D model provides improved waveform fits, with the variance reduction increasing from 21.7% to 47.3%. For the second event, the improvements were slight (Supplemental Figs. S1 and S21). Our misfit summary figures, presented as Supplemental Figures, show the variation of the misfit function over the space of moment tensor source types.

The results from Alvizuri and Tape (2016) show that the integrated waveform differences for high-frequency P waveforms were larger for the 1-D model than for the homogeneous half-space model. In their study, Alvizuri and Tape (2016) also observed that the synthetic waveforms generated with the 1-D model produced artifacts that were not in synthetic waveforms generated with the half-space model (within frequencies 2–10 Hz) and were not in the observed waveforms. Still we can examine the amplitude ratios between observed and synthetic P-waves for all 63 events. Figure 7 shows a spatial map of median amplitude ratios for the vertical component of the P-wave (Figs. 7A and 7C) and for the radial component of the P-wave (Figs. 7B and 7D). Perfect agreement between observed and synthetic amplitudes would represent ln(Aobs/Asyn) = 0 (ratio of observed to synthetic amplitude), plotted as solid blue colors. Comparison of Figures 7A and 7C shows that the 1-D model has the effect of reducing the systematic amplitude ratios identified on the vertical component for the half-space model. This shows that the relative weighting of the vertical component used by Alvizuri and Tape (2016) is probably not needed when using the 1-D model.

Figure 8 is a complementary view of the amplitude ratios in Figure 7. The histograms show an improved shift toward ln(Aobs/Asyn) = 0 for the vertical component, along with a slight shift away from zero for the radial component. Because our polarity measurements are made on the vertical component, we are most interested in fitting the waveforms on this component.

The main reason for the improved amplitude ratios with the 1-D model is that synthetic amplitudes derived with the 1-D model provided a better match to the observed amplitudes, especially in the radial P-wave component. Alvizuri and Tape (2016) suggested that this might be due to a shallow layer diffracting the seismic rays as they impinge the seismic stations. In the half-space model, the synthetic radial components are too large, but the inversion tends to fit these large observed waveforms. This results in the observed P vertical component waveforms being systematically higher than the synthetic waveforms. Hence Alvizuri and Tape (2016), using the homogeneous model, down-weighted the P radial component in the inversion.

Estimating full moment tensors for the Uturuncu events is challenging, mainly due to the limited signal at longer periods. Larger events at Uturuncu volcano would generate better signals within the period range (≥5 s) used in constructing the 1-D model.

We present a 1-D seismic profile for the crust and uppermost mantle below the Uturuncu volcano using the PLUTONS seismic array. Our new velocity model contains an 8-km-thick upper crustal lid and an ∼8-km-thick LVZ, which is consistent with earlier studies on a larger scale using ambient noise and PRFs (e.g., Ward et al., 2014). By also incorporating Rayleigh wave H/V ratios observed from earthquakes, we obtain further sensitivity to density and Vp/Vs ratio in the crust. Our results suggest that the LVZ has a slightly higher Vp/Vs and lower density compared to the rest of the crust, consistent with the hypothesis that the LVZ is related to partially molten magmatic complex (Sparks et al., 2008). The results from full moment tensor inversions (Alvizuri and Tape, 2016) demonstrate the potential benefit of using the 1-D model derived from this study to better explain surface wave waveforms and P wave vertical to radial amplitude ratio for regional earthquakes. A robust 1-D model extending to the uppermost mantle constructed by incorporating multiple seismic observables suggests that a similar method can be applied to other volcanic areas where a local, dense seismic array exists (e.g., Mount Erebus in Antarctica).

Caveats of this study extend from seismic observation to model inversion. First, we only make observation of average Rayleigh wave dispersion and H/V ratio measurements and only invert for the average 1-D model beneath Uturuncu. The spatial variation of the subsurface structure (particularly at shallow depths) can introduce complexity to the data for which we were not able to fit perfectly well. This caveat can be addressed by a two-dimensional surface wave tomography and construction of a three-dimensional model that captures the complexity more accurately, although the relatively small aperture of the PLUTONS array also makes it difficult to study lateral structural variation at a greater depth. Second, in this study the potential existence of radial anisotropy is ignored, although it can serve as an important constraint on the partially molten magma complex in the crust (Xie et al., 2013; Jaxybulatov et al., 2014; Hacker et al., 2014). To address this problem, future observation of Love wave speeds from ambient noise and earthquake should be considered. Third, we use a globally averaged attenuation model to perform the physical dispersion correction due to the lack of knowledge of local crustal attenuation structure, while the local existence of partial melt can greatly reduce the crustal Q. The systematic error of using an inaccurate Q model is not quantified in this study and should be considered for further investigation when a more accurate local Q model is available.

The authors are grateful to George Zandt and an anonymous reviewer for insightful comments that improved this paper. The facilities of the Incorporated Research Institutions for Seismology (IRIS) Data Management System (DMS) and specifically the IRIS Data Management Center were used to access the waveform and metadata required in this study. The IRIS DMS is funded through the National Science Foundation and specifically the GEO Directorate through the Instrumentation and Facilities Program of the National Science Foundation (NSF) under Cooperative Agreement EAR-0552316. This research was supported by NSF grants EAR-1215959 and EAR-0909254 at the University of Alaska, NSF grant CyberSEES-1442665 (to F.-C.L.) and the King Abdullah University of Science and Technology (KAUST) under award OCRF-2014-CRG3-2300 (to F.-C.L.), and NSF grant EAR-0952154 (to W.S.) at Washington University in St. Louis.

1Supplemental Figures. Misfit summary figures show variation of misfit function over space of moment tensor source types. Please visit http://dx.doi.org/10.1130/GES01353.S1 or the full-text article on www.gsapubs.org to view the Supplemental Figures.