New continuous GPS observations near the summit of Uturuncu volcano, Bolivia, indicate an average uplift rate of 2.4 ± 1.9 mm/yr between April 2010 and November 2015, while previous interferometric synthetic aperture radar (InSAR) observations between May 1992 and January 2011 predict an average vertical uplift rate of 7 ± 2 mm/yr. However, the GPS time series is better fit with a time-dependent function, such that the uplift rate in 2010 is less than 2 mm/yr and in 2015 may be as high as 9 mm/yr. Motivated by indications of a non-constant rate in the continuous GPS time series, we examine to what extent past InSAR measurements may have been temporally aliased. We present evidence that decreased uplift since 2004 is permitted by available InSAR and is consistent with a lower average GPS rate since 2010. We discuss how these variable rates may affect previously proposed magmatic models at Uturuncu including diapir ascent and reservoir pressurization. In particular, we explore a “dipole” reservoir model consisting of a magma source in the lower crust and sink in the middle crust for which variation in magma supply could explain sub-decadal fluctuations in deformation rate. Inversion of multiple InSAR data sets using homogeneous models constrains the deeper reservoir to 55–80 km depth (lower crust and upper mantle) and the shallower reservoir to 20–35 km. Consistent with previous work, the inferred ratio of volume change between the source:sink (i.e., deep:shallow) reservoirs is up to 10:1. We also incorporate new seismic tomography results in a 3D finite-element model to explore the combined effects of multiple sources and material heterogeneity on surface deformation. An important conclusion from our modeling efforts is that the commonly used diagnostic ratio of maximum radial to vertical surface displacements for shape of the deformation source is affected by both the presence of multiple reservoirs and subsurface heterogeneity. Ultimately, the dipole and diapir models have unresolved shortcomings given the entire suite of available geophysical data, but both provide valuable insight into the conditions for magmatic ascent in the Central Andes.

Uturuncu volcano (−67.18°W, –22.26°S, 6008 m), located in the Altiplano-Puna Plateau of the southern Bolivian Andes, has been the focus of recent research because it is located near the center of exceptionally broad (70-km-diameter), ∼10 mm/yr axisymmetric uplift (Pritchard and Simons, 2002). While the most recent eruptive product from the edifice of Uturuncu is Pleistocene in age (Sparks et al., 2008), the volcano is centrally located within a dense concentration of Miocene calderas and more than 15,000 km3 of silicic ignimbrite deposits known as the Altiplano-Puna volcanic complex (APVC) (e.g., de Silva, 1989; Salisbury et al., 2011). However, the current deformation footprint at Uturuncu does not appear to be spatially correlated with any particular caldera (Fig. 1). Uturuncu is also centrally located within the surface footprint of an extensive seismically imaged ultra-low velocity zone in the mid crust known as the Altiplano-Puna magma body (APMB) (e.g., Chmielowski et al., 1999; Ward et al., 2014). We summarize first-order geophysical features in a lithospheric cross section at the latitude of Uturuncu in Figure 2. Note that all depths in this paper will be with respect to average surface elevation, which is ∼4.6 km within 100 km of Uturuncu. In the Central Andes, petrological determination of pre-eruptive conditions implies that magma accumulates at different depths for varying periods of time (e.g., Kay et al., 2010), and the uplift at Uturuncu is believed to be related to partial melt in the APMB where andesitic magmas evolve to dacitic composition prior to ascent to shallower pre-eruptive reservoirs (Muir et al., 2014a, 2014b).

Interferometric synthetic aperture radar (InSAR) studies that synthesized data from May 1992 to January 2011 demonstrated that the zone of 10 mm/yr uplift at Uturuncu is surrounded by an approximately concentric ring of 1–4 mm/yr subsidence out to 75 km (Fialko and Pearse, 2012; Henderson and Pritchard, 2013). However, the rates reported in previous studies represent averages over the entire period of InSAR observations, and there was a suggestion of decreasing uplift toward the end of the observation period (Henderson and Pritchard, 2013). Here we present the first continuous GPS measurements near the summit of Uturuncu; these measurements indicate significant variability in vertical surface displacements since April 2010. Specifically, GPS positions indicate an average uplift of 2.4 ± 1.9 mm/yr over five years of observations, peak-to-peak annual seasonal variations of 20 mm, and evidence for acceleration to rates of 9 mm/yr in the past several years. Accurately determining the temporal variation of surface uplift is important for guiding geodynamic models that have previously attempted to explain constant uplift rates over decades due to increasing reservoir pressurization (e.g., Pritchard and Simons, 2004; Hickey et al., 2013) or buoyant diapiric ascent (Fialko and Pearse, 2012; del Potro et al., 2013).

Considering the evidence for mid-crustal melt, observed deformation at Uturuncu was previously hypothesized to be caused by magmatic intrusion into or out of the APMB (e.g., Pritchard and Simons, 2004; Hickey et al., 2013). However, these studies were undertaken prior to the discovery of co-located subsidence and therefore considered uplift only. In fact, mechanical models are generally based on pressurization of a single chamber due to added mass, with the implicit assumption that the depressurization of a deeper feeder chamber would be too deep to generate a detectable surface deformation signal. Nevertheless, one possible explanation for the unique “sombrero” pattern (axis-symmetric uplift surrounded by subsidence) is sustained deep depressurization superimposed on shallower pressurization (Fialko and Pearse, 2012; Henderson and Pritchard, 2013). We refer to this coupled reservoir model as a magmatic “dipole” in which the “source” reservoir deflates in the lower crust and the “sink” reservoir inflates in the middle crust. Fialko and Pearse (2012) considered a magmatic dipole model but ultimately favored a model of diapiric ascent occurring in the middle crust as more physically plausible. Importantly, all geodetic models to date have assumed either homogeneous crustal properties or layered properties constrained by seismic receiver functions, which prominently feature a 1–2-km-thin APMB (e.g., Leidig and Zandt, 2003; Zandt et al., 2003). An updated seismic-velocity model using both receiver functions and ambient noise tomography has significantly altered the APMB geometry from a thin sill to a thickened low-velocity (Vs < 2.9 km/s) zone between ∼9 km and 29 km depth (Ward et al., 2014).

In the following sections, we first describe geodetic data available for Uturuncu, including InSAR and two new continuous GPS stations. We show that these data sets are self-consistent if the deformation rate at Uturuncu is temporally variable, and thus the previously inferred constant uplift rate represents an average over periods of variable uplift rate. We then revisit a simple dipole reservoir model, which has the advantage of matching observed changes in deformation rates simply by coincident changes in reservoir overpressure. We also consider the effect of subsurface heterogeneity on surface deformation by incorporating the seismic tomography of Ward et al. (2014) into a finite-element model. In our discussion, we assess the merits and drawbacks of both the dipole and the diapir models, and we ultimately suggest that both models fall short of a satisfying explanation of the expanded geodetic measurements.

2.1 InSAR

Henderson and Pritchard (2013) and Fialko and Pearse (2012) present independent analyses of all available strip-map C-band InSAR data for Uturuncu from the European Remote Sensing (ERS-1/2) and Envisat satellites, spanning May 1992–January 2011; in addition, Fialko and Pearse (2012) included SCANSAR Envisat data. Observations from other SAR satellites (ALOS, Radarsat-2, TerraSAR-X, and COSMOSky-Med) operating after 2010 have few acquisitions at Uturuncu, and consequently, few interferograms are available after October 2010. Models in this paper are therefore compared against the best available data from ERS and Envisat satellites: stacked InSAR data from two ascending tracks (89 and 3) and two descending tracks (282 and 10) (Table 1). This data set provides a map of average surface ground velocities at ≈1 km2 resolution with estimated uncertainties of ±4 mm/yr (Henderson and Pritchard, 2013). The unique deformation pattern at Uturuncu is well characterized by a radially averaged profile (although with some second-order spatial variations), which shows uplift and subsidence that is not correlated with elevation (Fig. 3).

Decomposing InSAR line-of-sight (LOS) measurements into three Cartesian displacement vectors can help distinguish between reservoir geometries in the subsurface (e.g., Dieterich and Decker, 1975). The vertical or radial displacement profile alone is insufficient to identify reservoir geometry because the aspect ratio and depth of ellipsoidal sources can be varied such that profiles are indistinguishable (e.g., Fialko et al., 2001). We therefore attempted to decompose the four different LOS observations of deformation at Uturuncu into three (east-west, north-south, and vertical) components of deformation following the procedure in Wright et al. (2004). However, due to the near-polar orbits of SAR satellites, north-south deformation is poorly constrained. In fact, for the four viewing geometries at Uturuncu, ±1 mm root-mean-square error (RMSE) in LOS displacement scales to a very large uncertainty in the north-south displacement vector of ±20 mm. By assuming that surface displacements are radially symmetric, we invert for a representative profile along an east-west transect (where north-south displacements are assumed to be zero). Figure 4 shows the profile of vertical and radial displacements constrained to the region where all four InSAR tracks overlap. The four profiles from each stack are shown in Figure 5, along with the best-fit solution re-projected into the respective viewing geometries.

2.2 Continuous GPS

We installed three continuous GPS stations in April 2010 named UTUR, COLO, and QUET to help constrain the temporal deformation at Uturuncu (Table 2). Sufficient time has now elapsed to derive ground velocities from the time-series data for UTUR and COLO, but QUET was vandalized after just one year and was not repaired. Daily positions between April 2010 and November 2015 were calculated by the University of Nevada Geodetic Laboratory using the GISPY OASIS II (version 6.1.1) software package. Final, non-fiducial daily orbit products from the Jet Propulsion Laboratory (JPL) archive were used in processing. Coordinates are relative to the IGS08 reference frame, which is a GPS-only realization of ITRF2008 (Altamimi et al., 2011). Further details on processing technique can be found in Blewitt et al. (2016). Time series for UTUR and COLO are presented in Figure 6 with best-fit rates reported in Table 3. Vertical daily positions at UTUR are shown in more detail in Figure 7. Daily position estimates for both UTUR and COLO have standard deviations of 0.5 mm, 0.5 mm, and 2.0 mm on the EW, NS, and Z components, respectively. In addition to the anticipated volcanic deformation signals previously observed with InSAR, the GPS stations record interseismic deformation, seasonal hydrological loading, and other forms of noise that are further discussed in the Results section.

3.1 Surface Displacement Components from InSAR

The inversion for Cartesian components is ideally conducted with data sets spanning the same time range with similar acquisition spacing. Unfortunately, our InSAR stacks cover different time spans and have different numbers of acquisitions (Table 1), and we therefore make a critical assumption that deformation has been constant over the total extent of observations. However, in the inversion for the radial and vertical displacement components, we obtain an unrealistic result of the minimum radial displacement being shifted in comparison to the maximum vertical displacement (Fig. 4). Simple subsurface pressurization and symmetrical loading models predict that maximum vertical displacement should coincide with zero radial displacement. However, forcing this constraint causes a misfit in the ascending profiles (Fig. 5). Specifically, the predicted peak LOS displacement rate for ascending track 3 underestimates observations by ∼1 mm/yr, and the predicted rate for ascending track 89 overestimates observations by ∼1 mm/yr. These shifts could be explained by temporal variation in uplift rate since each track covers different timespans (section 3.2).

We first emphasize that the maximum vertical displacement (−67.233°W, –22.270°S) is located ∼6 ± 2 km west of the summit of Uturuncu, in agreement with offsets of ∼5 km west of the summit reported in previous studies (Sparks et al., 2008; Fialko and Pearse 2012; Henderson and Pritchard, 2013). Based on our best-fit inversion result, the peak eastern radial displacement is 3.0 mm/yr, and peak vertical displacement is 9.1 mm/yr, corresponding to a ratio of ur/uz = 0.33 (Fig. 4). However, the radial displacements show asymmetry to the west, with a peak of 2.5 mm/yr (ur/uz = 0.27), and given uncertainties in our inversion, we determine a ur/uz ratio of 0.30 ± 0.15.

3.2 Temporal Variability in InSAR Deformation Rates

As discussed in the previous section, in order to achieve aligned vertical and radial profiles, the uplift rates in ascending track 3 must be reduced by up to 1 mm/yr, and the observed rates for track 89 must be increased by up to 1 mm/yr. We believe such shifts are not unreasonable given the different temporal span of each track (Table 1) and noise in the data. In particular, track 89 spans 2004–2010, whereas all other tracks have data going back to the 1990s. Track 3 may be artificially high because only 15 interferograms were used in the stack (with only three dates before 2004), compared to hundreds of interferograms being averaged for descending data with a similar time span. Interferogram stacks by design produce average velocity values, and while a constant rate has been shown to approximate the data at Uturuncu (Henderson and Pritchard, 2013), the time series also permits a change in rate between 2003 and 2005, a period in which there were sparse observations. To demonstrate, we re-plot the time series at the location of continuous GPS station UTUR with linear fits to data before and after 2004 (Fig. 8). Line-of-sight uplift before 2004 is predicted to be on the order of 10 mm/yr (discounting track 3, which is very poorly constrained by three dates); whereas after 2004, LOS rates are predicted at less than 6 mm/yr. This observation supports a decelerating trend in uplift, which is reflected in the more recent continuous GPS measurements described next.

3.3 Continuous GPS Rates

Time-series analysis of InSAR assuming a constant uplift rate predicts ∼7 ± 2 mm/yr vertical uplift and less than 2 mm/yr horizontal motion at UTUR and vertical subsidence of –2 mm/yr at COLO (see Fig. 4). In order to check these predictions against the continuous GPS data, displacements due to interseismic motion and seasonal hydrologic signals must be taken into account. There are other non-volcanic signals in the GPS time series such as atmospheric loading, local bedrock thermal expansion, and daily positioning errors, which can vary significantly based on processing strategy (Dong et al., 2002). However, these signals should be second order compared to the hydrologic signal of loading in the Amazon Basin (e.g., Fu et al., 2013). We use a kinematic approach to remove best-fit functions from the time series of regional stations (Bevis and Brown, 2014; Blewitt et al., 2016).

3.3.1 Interseismic Signal

Stations COLO and UTUR are 350 km and 410 km, respectively, east of the Nazca–South America trench, such that constant-rate interseismic convergence is expected to be the dominant signal in horizontal displacements. We perform a simple linear regression for each GPS component to obtain first-order estimates of constant rates and 95% confidence intervals (Fig. 6 and Table 3). However, this regression does not account for co-seismic offsets or the fact that daily positions are temporally correlated. Therefore, we also use the median interannual difference adjusted for skewness (MIDAS) estimator for interseismic velocities because it is designed to account for steps and periodic annual signals in the time series (Blewitt et al., 2016). Horizontal vectors at UTUR and COLO seem consistent with the regional velocity field (–24 < latitude < –16, and −72 < longitude < −63), although our stations have slightly elevated eastward velocities (Fig. 9 and Table 3). Few stations are available beyond a distance of 350 km from the trench, and thus prediction of rates at UTUR and COLO is prone to extrapolation error. In general, MIDAS vertical velocities at regional stations vary between 0 and 5 mm/yr; so it is possible that the vertical velocity we calculated for UTUR, 2.4 ± 1.9 mm/yr, is purely due to interseismic motion (Fig. 9). Interestingly, COLO is one of the few stations with a negative MIDAS vertical velocity (–3.2 ± 3.4 mm/yr), but this best-fit rate is clearly worse than alternative estimates discussed in the next section (Fig. 6). We believe that the MIDAS estimator preforms poorly at COLO because it selects a small number of valid daily position pairs to compute a rate—only 220 pairs compared to 705 pairs at UTUR.

3.3.2 Pisagua Co-Seismic Offset and Other Interseismic Estimates

Stations UTUR and COLO are located >400 km to the SE of the Mw 8.2, 1 April 2014 Pisagua earthquake epicenter, but offsets are clearly visible in each time series (Fig. 6).

In order to quantify the offset associated with this earthquake and more accurately estimate the constant background rate for each GPS component, we fit a kinematic model with seven parameters: initial position, linear rate, co-seismic offset (as a Heaviside step function), and four Fourier series coefficients (Bevis and Brown, 2014). We fix the initial time to 15 April 2010, the time of co-seismic offset to 1 April 2014, and set the seasonal periods to 0.5 yr and 1 yr. The linear rates are compared in Table 3. In general, this model provides the best qualitative fit to the GPS time series (Fig. 6) and results in the lowest residuals (Table 3), which is not surprising given that it has more model parameters. However, the model still assumes that a constant rate applies to the entire time period analyzed. For UTUR, the constant rate is 3.2 ± 0.2 mm/yr, and the rate for COLO is –0.3 ± 0.3 mm/yr. Co-seismic offsets at regional stations due to the Pisagua earthquake are shown in Figure 10. It has been shown that large earthquakes can cause enhanced subsidence above magmatic reservoirs compared to surrounding areas (Pritchard et al., 2013). However, we do not detect significant deviations at COLO and UTUR compared to other regional stations.

3.3.3 Possible Accelerating Uplift at UTUR

Qualitative inspection of the vertical component of displacements at UTUR suggests a possibly higher uplift rate in the latter half of the time series. Indeed, residuals for a linear rate suggest some curvature indicating a constant rate may not be an appropriate model (Supplemental Material1). Polynomial trajectory models are also commonly used in GPS time-series analysis (Bevis and Brown, 2014). We therefore fit a quadratic function to the vertical time series and present residuals in Figure 7. The quadratic function results in a RMSE of 5.7 mm compared to the multi-parameter fit (section 3.3.2) RMSE of 5.4 mm with significantly fewer parameters. Importantly, the best-fit quadratic model implies a yearly increase in background rate of 2 mm/yr. After three years, the instantaneous vertical velocity is 3 mm/yr, in agreement with the 5.5 year averages determined by the MIDAS and multi-parameter methods. However, by the end of the fifth year, the vertical rate has increased to 9 mm/yr. It is worth noting that many other functional forms could be fit to the GPS time series, for example, a piecewise linear function as we have proposed for the InSAR time series. Fitting models to GPS time series spanning less than 2.5 years is susceptible to significant bias from seasonal signals (Blewitt and Lavallée, 2002); therefore, for illustration, we break the GPS time series at UTUR in the middle, 1 February 2013, such that 2.8 years of data are available before and after. We find an earlier rate of 1.5 ± 0.4 mm/yr followed by a later rate of 6.9 ± 0.5 mm/yr, with a comparable residual of 5.7 mm. Both examples illustrate a secular trend of accelerating displacement at UTUR over 5.5 years of observation.

3.3.4 Hydrologic Seasonal Signal

After removal of the quadratic trend in the time series, periodic seasonal signals in the vertical components of both UTUR and COLO are pronounced, with peak-to-peak vertical displacement amplitudes of ∼20 mm (Fig. 7). Minima are broad, but maxima tend to occur in late December. The annual variation cannot be explained by local hydrological loading, since maximum (<400 mm/yr) regional rainfall occurs mostly in the months of December, January, and February (Garreaud et al., 2003). Instead, annual signals appear to be correlated on most regional stations and have been explained by the elastic response due to hydrological loading of the Amazon Basin (Davis, 2004). Total water storage (TWS) in the Amazon Basin lags maximum rainfall and typically peaks in March-April, with a minimum in November-December (e.g., Chen et al., 2010). Minimum TWS theoretically causes maximum seasonal uplift due to the elastic unloading effect, which has been documented for the POVE GPS station with maxima in late October and minima in late April (Fig. 11). We note a lag time of ∼1 month in the POVE time series versus UTUR and COLO (Fig. 11). This lag is most likely due to the POVE station being located over 1500 km NNE of Uturuncu and maximum water storage being a non-stationary surface load as the basin drains to the Pacific Ocean. We rule out any local seasonal or nonlinear process at UTUR (perhaps due to magmatic processes), because the same maxima timing and seasonality is observed at other regional stations (Fig. 12).

We seek geodynamic models that can reproduce sub-decadal changes in surface uplift rate but with a secular uplift trend. In addition, we suspect that the new seismic tomography, in particular the change in magnitude and geometry of the APMB low-velocity anomaly (Ward et al., 2014), may affect previously proposed reservoir pressurization models. We therefore reevaluate an elastic dipole model, in which displacement rates are proportional to magma supply (volume change is equated to reservoir overpressure) both in a homogeneous half-space and incorporating subsurface heterogeneity with a finite-element mesh.

Previously, Fialko and Pearse (2012) proposed a dipole model that well explained the InSAR data; nevertheless, they rejected the model as physically untenable. We review their criticisms of this model, but we also argue that the dipole deserves further consideration.

  1. Given the evidence for crustal partial melt and the shallow brittle-ductile transition, an elastic model is not appropriate for the middle to lower crust beneath Uturuncu. Viscous creep processes take place at depth over timescales determined by temperature, strain rate, rheology, and composition. Consequently, if a material deforms viscoelastically, short-term elastic strain resulting from transient stresses should be expected for t << τ, where τ is the effective relaxation time for the material. The unknown thermal and rheological structure of the crust is therefore critical in determining the appropriate deformation for Uturuncu. A weak felsic crust predicts a shallow brittle-ductile transition in agreement with seismicity (Jay et al., 2012), but a more intermediate middle and lower crustal composition may be justified given the andesitic composition of the APMB (Comeau et al., 2015). Independent constraints on the viscosity of the lower crust in the Central Andes are limited, but one study of lake loading indicates viscosities less than 1020 Pa s below 40 km depth (Bills et al., 1994). Given an estimated shear modulus of 10–100 GPa for the lower crust, these viscosity values imply relaxation times on the order of 10–1000 years, which is similar to the timescale of ground deformation measured by geodesy and geomorphology (e.g., Perkins et al., 2016). On the other hand, much lower viscosity values of <1016 Pa s are estimated by Gottsmann et al. (2017) and <1014 Pa s in the partially molten APMB by Comeau et al. (2016). Considering the uncertainty in both viscosity and appropriate rheological law for the APMB, we think it is worth exploring elastic models to understand their implications and determine if there are other reasons these models could be discarded.

  2. There is a volume change discrepancy between lower and deeper sources. Fialko and Pearse (2012) modeled the inflating source at 25 km and the deflating source at 80 km and found that the volume of the lower source must be four times larger. We here explore a wider range of models that can match the observations and discuss mass balance between magmatic reservoirs in more detail.

  3. Connected sources of uplift and subsidence would imply a permanent magma conduit or quasi-steady dike intrusions between the middle and lower crust, which Fialko and Pearse (2012) did not think plausible. Vp/Vs tomography has been proposed as evidence for sustained magmatic pathways from the lower crust to shallow chamber elsewhere, for example at Mount St. Helens (Kiser et al., 2016). There is new geophysical evidence from earthquake tomography (West et al., 2013) for a low seismic velocity zone (with Vp/Vs > 1.9) connecting the middle- to lower-crust Uturuncu, suggesting that there could be such a necessary pathway. Although the exact depth extent of this feature is not well resolved because of the vertical smearing inherent in the tomography, based on the observations, we do not let this argument reject the dipole model a priori.

  4. The very low seismic velocities of the APMB imply a low shear modulus that may not allow the transmission of strain from the lower crust to the surface. While we agree with Fialko and Pearse (2012) that this decoupling of the lower crust is possible, it has not yet been rigorously shown. For example, the APMB of Ward et al. (2014) does not have as low a seismic velocity as the APMB model used by Fialko and Pearse (2012); thus, the amount of melt and the shear modulus may not be as low and therefore could impact surface displacements. Future work is necessary to address this question, but here we consider the implications of the scenario in which strain is transmitted from the lower crust through the APMB to the surface.

4.1 Homogeneous Half-Space Model

We obtain first-order estimates of the depth and volume of dual-point reservoirs in an elastic half-space (Mogi, 1958). We apply a Levenberg-Marquardt nonlinear least-squares algorithm to solve for source location and volume change (Fig. 13). The inversion assumes that the crust is a Poissonian elastic half-space, which is known to be an oversimplification (e.g., Masterlark, 2007). Nevertheless, the model serves as a comparison to previous studies and to more realistic heterogeneous models that we introduce in this paper.

We inverted descending interferogram stacks from tracks 282 and 10, which have a higher signal-to-noise ratio compared to ascending data because they contain hundreds more interferograms. We estimate an inflation source at 20 < Z < 35 km depth with 0.01 < ΔV < 0.1 km3/yr and a deflation source at 55 < Z < 80 km depth with -0.01 < ΔV < –1.0 km3/yr. A range in parameter values is reported based on model results obtained with variable degrees of data resampling and different initial parameter values. Furthermore, it is well known that a single Mogi source is non-unique and has a trade-off between volume and depth. Introducing a second source only compounds this issue because the sources counteract one another. Nevertheless, a map view and profile view comparing descending data to the best-fit model prediction demonstrate the ability of this model to reproduce observations (Fig. 13). The map position of the best-fit inflation source (67.22°W, 22.27°S) is shifted ∼5 km to the southwest relative to the summit of Uturuncu (67.18°W, 22.26°S) as was found in previous work (e.g., Pritchard and Simons, 2004) and in agreement with the location of maximum vertical displacement based on our inversion for surface displacement components. We did not force sources to be exactly vertically aligned in the inversion but find that the deeper deflation source is vertically aligned to within 5 km of the position of the shallower inflation source.

Inversions for multiple sources in a half-space must be interpreted with caution since the presence of a second source violates the boundary condition assumptions. Pascal et al. (2014) quantified the conditions under which multiple source inversions are valid by numerically modeling a suite of dual-reservoir configurations. According to this study, as long as vertically aligned sources are greater than 8 times the largest source radius, inversion results are not significantly biased. “Point” sources are often dismissed as unphysical; however, they represent first-order mathematical approximation of a deep (a << d) spherical reservoir with a radius determined by the relationship δV = πa3δp/µ (Mogi, 1958). If we consider an overpressure bound of 10 MPa, added volume of 0.1 km3, and shear modulus of 40 GPa, we estimate the radius for larger (lower) reservoir to be ∼5 km, and thus source interaction is non-negligible for some ranges in our homogeneous inversions. We therefore suggest that the depths and volumes in the homogenous inversions previously presented can be underestimated by up to 30% (Pascal et al., 2014). The volume and overpressure can also be modulated by assuming viscoelastic conditions (e.g., Hickey et al., 2013), and therefore we do not place a strong emphasis on a single “best-fit” model and instead emphasize the range of solutions to guide further exploration with a finite-element model.

One major implication of interacting superimposed sources is that the ratio of ur/uz is increased compared to the ratio for a single reservoir. Previous studies have shown that, in a homogeneous medium, prolate sources have ur/uz > 0.5, and oblate sources (sill-like) are characterized by ur/uz < 0.3 (Dieterich and Decker, 1975). For our best-fit model at Uturuncu, the predicted ratio is 0.46—suggesting that surface displacements could either be due to a single prolate source or a dipole. Forward modeling indicates that ur/uz increases as the source and sink reservoirs move closer together and as the ratio of source to sink volume increases. Importantly, the presence of a source reservoir does not decrease ur/uz below 0.3, and therefore inversions for a single source should not be biased toward sill-like geometries. While informative for future modeling results in other locations, this simple analysis is restricted to certain conditions that may not be broadly applicable—specifically, we have assumed a point-source approximation, vertical alignment, and elastic, homogeneous, Poissonian crust.

4.2 Heterogeneous Finite-Element Model

On elastic timescales, a weak subsurface layer is known to amplify radial and vertical strain disproportionally from pressurization at depth (e.g., Geyer and Gottsmann, 2010). In order to investigate the role of 3D heterogeneous material properties at Uturuncu, we constructed a finite-element model (FEM). Our model mesh consists of a quarter cylinder with a radius of 300 km and a depth of 200 km (Fig. 14). The quarter-cylinder design permits exploration of 3D subsurface heterogeneity using the seismic tomography of Ward et al. (2014), while maintaining reduced computational expense by taking advantage of mirror plane symmetry. Previous finite-element modeling at Uturuncu utilized an axis-symmetric formulation to explore 1D subsurface layering (Hickey et al., 2013), and a companion paper in this issue (Gottsmann et al., 2017) uses the model of Ward et al. (2014).

In our model, a sink reservoir is located at 20 km depth and the source reservoir at 70 km depth along the cylinder axis. We choose these depths because they are coincident with observed seismic anomalies (Yuan et al., 2000; Ward et al., 2014) and represent positions in the range of permissible solutions from half-space inversions (section 4.3). Boundary conditions include outward normal tractions applied to the reservoir surfaces, the upper surface is free, no displacement is allowed normal to the cut-out cylinder sides, and zero displacement is enforced on the outermost curved surface. The mesh was generated with CUBIT software and consists of 106 linear tetrahedral elements. Numerical solutions of the mechanical response are performed with Pylith software version 1.9 (Aagaard et al., 2013). Solutions are benchmarked against superimposed analytical solutions for Mogi point sources in a Poissonian half-space (Vp = 5773 m/s, Vs = 3333 m/s, ρ = 2700 kg/m3, E = 75 GPa, G = 30 GPa, ν = 0.25; see Supplemental Material [footnote 1]).

We used available seismic imaging data in order to constrain material properties in the finite-element model. Early seismic surveys indicated a thick crust of more than 70 km and low seismic velocities (Vp = 6.2, Vs = 3.45 km/s) underlying the Altiplano Plateau (James, 1971). More recently, receiver-function analysis suggested a 1-km-thick feature at 19 km depth with the following properties: Vp = 4 km/s, Vs = 1 km/s, ρ = 2400 g/cm3 (Chmielowski et al., 1999; Leidig and Zandt, 2003; Zandt et al., 2003; Fig. 15). A roughly coincident low-resistivity layer led to the hypothesis that this feature is a sill of >10% andesitic partial melt known as the Altiplano-Puna magma body (APMB) (e.g., Schilling et al., 1997). The sill-like APMB has figured strongly into previous geodynamical studies (e.g., Babeyko et al., 2002; de Silva and Gosnold, 2007; Fialko and Pearse, 2012). However, recent joint inversion of ambient noise tomography and receiver functions suggests a drastic change geometric character of the APMB. It is now believed to be a much thicker zone of low velocities (1.9 < Vs < 2.9 km/s) between 10–30 km depth, and a diameter of over 200 km (Ward et al., 2014) but with less extreme low Vs values (Fig. 15). In addition, new magnetotelluric surveys suggest low resistivity in the middle to upper crust and may represent a combination of melt fractions greater than 20% in addition to overlying brine solutions (Comeau et al., 2015; Comeau et al., 2016).

The seismic-velocity models discussed so far are synthesized in Figure 15. The new characterization of the APMB partly results from incorporating additional data from an array of 28 broadband seismometers deployed between 2010 and 2013 with 50 km of Uturuncu. Ward et al. (2014) assume a Poissonian crust (Vp/Vs = 1.73), although recent local body-wave tomography shows perturbations of Vp/Vs on the order of ±5% within a 50 km radius of Uturuncu (West et al., 2013). Such decreases in the Vp/Vs ratio are expected to slightly increase depths to certain layers in the seismic model, but detailed studies have yet to be done.

Each element in our model is assigned elastic properties based on direct conversion of Vp, Vs, and density from the tomographic velocity model in Ward et al. (2014). The primary measurement from Ward et al. (2014) is Vs, and Vp is obtained given the Poissonian assumption (Vp/Vs = 1.73). Finally, density is assumed constant at 2700 g/cm3 throughout the crust. We remark that “dynamic moduli” derived from seismic velocities are generally higher than the “static moduli” derived from laboratory rock deformation experiments. Since these discrepancies can be at least partly attributed to micro-crack density (e.g., Cheng and Johnston, 1981), better agreement is anticipated with increasing depth due to crack closure. Reductions to dynamic moduli have been proposed in some volcano deformation studies as a means to artificially simulate longer duration anelastic crustal response (e.g., Trasatti et al., 2005). Therefore we consider that the elastic material properties in our model are upper bounds and note that potential reductions to moduli amplify surface deformation, all else being equal.

For the finite-element model, we assigned material properties based on a 3D tomographic data set (Ward et al., 2014). We note that the observed velocity variations are predominately a one-dimensional velocity gradient with respect to depth (Fig. 15). The extent of the APMB low-velocity zone exceeds the deformation footprint observed at Uturuncu, and the characteristic length scale of lateral variation in seismic velocities is tens of kilometers (Ward et al., 2014).

Including this broad weak zone between the inflating reservoir and the surface causes the ratio ur/uz to decrease to 0.30, representing a significant 35% decrease compared to the homogeneous case (0.46) (Fig. 16). This value is in better agreement with the values obtained in our preferred inversion for Cartesian components in section 3.1 (0.3 ± 0.15). Maximum uplift is also increased threefold, such that the inferred volume from homogeneous inversions may be overestimated by the same amount. The effect is amplified if the layer is thicker, made of weaker material, and closer to a pressurizing reservoir. For example, moving the inflation source to a depth of 30 km decreases the ratio ur/uz by 25% compared to the homogeneous case, demonstrating how the effect of the weak layer on surface deformation diminishes for reservoirs of increasing depth. Similar findings have been reported for inversion schemes with 1D layering (e.g., Pritchard and Simons, 2004).

On elastic timescales, we confirm that a weak layer amplifies both vertical and radial strain from lower sources, but vertical displacements are affected to a greater extent (e.g., Geyer and Gottsmann, 2010; Hautmann et al., 2010). This effect is partially offset by including a deeper source reservoir. In the specific case we have tested, the APMB weak zone produces a greater reduction in ur/uz compared to the increase in ur/uz generated by inclusion of a magmatic dipole in a homogeneous space.

5.1 Mass Balance between Dipole Reservoirs

Evidence for vertically separated reservoirs has been presented for other volcanic systems (for example, Soufriére Hills, Elsworth et al., 2008; Fernandina, Bagnardi and Amelung, 2012; and Grimsvötn, Reverso et al., 2014). But the dual-reservoir model at Uturuncu is unique in that reservoirs are believed to be comparatively deep. A regional seismic Moho at 60–70 km is suspected (Yuan et al., 2000), and therefore our source reservoir falls within the range of lower crust and upper mantle. We envision a lower reservoir containing basalt-to-andesite melt associated with “MASH zone” processes above the Moho (Hildreth and Moorbath, 1988) and an upper reservoir containing melt of andesite-to-dacite in the middle crust.

We must reconcile the fact that the dipole models estimate the ratio of volume of the deep:shallow (source:sink) reservoirs is up to 10:1, which represents an upper bound from models discussed in previous sections. The same volume mismatch has been pointed out for a dual-reservoir model at Soufriére Hills volcano for reservoirs in the middle crust to shallow crust (Foroozan et al., 2010). Volume ratios of ∼1:5 have been observed for shallow dike inflation compared to adjacent chamber deflation (e.g., Wright et al., 2006). This ratio has been convincingly explained by compressibility arguments (Rivalta and Segall, 2008); however, the inverse ratio is observed for the dipole model at Uturuncu and Montserrat. Here we discuss under what conditions such a volume discrepancy could arise.

Following notation in Segall (2010), incremental mass in an intrusion is found by applying the chain rule to m = ρmV:
Equating incremental mass between deep source (δm) and shallow sink (δm+) reservoirs, it is possible to solve this for a volume ratio. We incorporate the relations for magma compressibility βm = 1/ρm * ∂ρm/∂ρ, spherical chamber compressibility is βc = 3/4µ (e.g., Tiampo et al., 2000), and reservoir volume change due to inflation δV = πa3δπ/µ (Mogi, 1958). As a result, for connected source and sink reservoirs the volume ratio becomes:

Inversion of geodetic data at Uturuncu suggests δV/δV+ > 1. In general, as depth increases, magma density increases (ρm+m < 1), but magma compressibility decreases (βm+m > 1), and reservoir compressibility also decreases (βc+c > 1). Critically, the theoretical volume ratio depends on the ratio of magma compressibility to chamber compressibility (βmc) at a given depth. Assuming andesite magma, the density should decrease as it ascends such that (ρm+m ∼0.9) (e.g., Spera, 2000). Therefore, to explain a (1 < δV/δV+ < 10), compressibility changes must offset the density effect. We estimate chamber compressibility by converting the crustal velocity model of Ward et al. (2014) to elastic moduli, such that µ ≈40 GPa and µ+ ≈12 GPa, corresponding to βc = 1.9 × 10–11 Pa–1 and βc+ = 6.3 × 10–11 Pa–1. If we assume these values for wall-rock compressibility and estimate a source magma compressibility of βm = 1 × 10–11 Pa–1 (Huppert and Woods, 2002), the corresponding sink compressibility is required to be βm+ = 3 × 10–10 Pa–1. However, we do not anticipate significant exsolution of volatiles at mid-crustal depths, and therefore we do not expect large changes in magma compressibility. In general, the volume ratio could be increased at the sink reservoir (βm+ > βc+), for example by exsolution of volatiles in a stiff reservoir. Alternatively, the volume ratio can be increased by source reservoir properties (βm < βc), for example by invoking a more ellipsoidal chamber that is much more compliant compared to the spherical assumption (e.g., Amoruso and Crescentini, 2009).

We note that only if shear modulus, density, and magma compressibility remain relatively constant, mass balance is equivalent to volume balance. Interestingly, for this particular case, vertically aligned source and sink pressure sources separated by a sufficiently small distance produce nearly equal and opposite strain fields. As a result, calculated surface deformation is small and possibly below geodetic detection thresholds (Supplemental Material [footnote 1]). As a corollary, subsidence moat signals from vertically aligned sources would only appear if the ratio of volume changes is sufficiently large. These two points are important, since they may factor into the lack of deformation observed to date in some arcs (Ebmeier et al., 2013) or the general scarcity of subsidence signals accompanying episodes of volcanic uplift. However, the assumption of vertical alignment is not always valid, and analytical predictions for surface displacements do not consider source interaction (e.g., Pascal et al., 2014).

5.2 Temporal Variation in Deformation

The discrepancy between the average uplift rate from InSAR data (7 ± 2 mm/yr, May 1992–January 2011) at UTUR and slower constant rate measured by continuous GPS (2.4 ± 1.9 mm/yr, April 2010–March 2015) is surprising at first. However, we have seen from our InSAR and GPS time-series analysis that the secular decadal uplift rate at Uturuncu averages out variation on timescales of five years. For example, daily positions from 5.5 years of daily continuous GPS solutions reflect an average vertical uplift rate on the order of 1.5 ± 0.4 mm/yr between 15 April 2010 and 1 February 2013, followed by an average vertical uplift rate of 6.9 ± 0.5 mm/yr through 24 November 2015. Importantly, we do no detect sustained periods of subsidence at UTUR.

Past studies have emphasized constant rates (e.g., Fialko and Pearse, 2012; Henderson and Pritchard, 2013); however, we note that the constant rate is representative of cumulative deformation over two decades. We have presented evidence that uplift at Uturuncu was less than 2 mm/yr between April 2010 and February 2013 but has been increasing once again in recent years (Fig. 7) such that over a decadal observation period, the average GPS uplift rate may in fact match previous InSAR observations. In summary, we suggest periods of negligible uplift over several years are real and are reflected in both the GPS and InSAR.

5.3 Modeling Temporal Variation

If temporal variations in uplift are real, they can be attributed to changes in subsurface magmatic processes. Time dependence of surface displacements for elastic models is achieved exclusively via the time dependence of reservoir pressurization. To match constant cm/yr uplift rates at Uturuncu requires that overpressure also increase at approximately ≈1–10 MPa/yr (or higher—Hickey et al., 2013). One difficulty with this model is the need for a quasi-static transport mechanism between source and sink reservoirs. Reservoir inflation via sustained pipe flow through the crust has been shown to reproduce surface uplift at other volcanic zones, such as Laguna del Maule (Le Mevel et al., 2016). Alternatively, quasi-steady intrusion of dikes toward a long-lived reservoir at the level of the APMB could explain observations (e.g., Karlstrom et al., 2009). Assuming fluctuations of surface displacement are due to transient mass addition, we require that transport must be occurring at rates on the order of tens of kilometers per year, either as quasi-steady flow or semi-regular dike intrusions focused into a reservoir zone over multiple decades.

We have limited the models presented in this paper to an elastic rheology; however, viscoelastic subsurface rheology has been used in other studies. For example, Hickey et al. (2013) invoked a standard linear solid mid-crustal viscoelastic rheology to explain quasi-steady surface uplift with either a stepped or linearly increasing reservoir overpressure. Without transient motions, it is difficult to disentangle the effects of time-dependent reservoir overpressure versus time-dependent viscous deformation. The diapiric ascent model presented in Fialko and Pearse (2012) predicted that surface uplift would gradually decelerate and expand given no magmatic recharge to the APMB. This prediction is also contingent on the entrainment of partial melt being contained to a 1-km-thick zone representing the APMB with viscosities predicted by a power-law relationship that are on the order of 1016 Pa s.

5.4 Longer-Term Outlook

A natural question is how long does magma supply persist at Uturuncu, and what is the accumulated surface displacement that can occur? A recent study by Walter and Motagh (2014) describes a “fracture girdle” that surrounds Uturuncu and is consistent with depressurization at APMB depths and, consequently, broad net subsidence at the surface over million-year timescales. On the other hand, differential GPS measurements of lake shorelines within the deformation field do not show significant tilting, suggesting the current unrest episode would have a maximum duration of 100 years (Perkins et al., 2016). Thus, there is evidence for changing rates and maybe even sustained subsidence on the century-to-millenia timescale at Uturuncu. Coupled with the changing deformation seen at other volcanoes in the Central Andes and globally (e.g., Fournier et al., 2010; Henderson and Pritchard, 2013), we should not be surprised to see temporal variations in deformation at Uturuncu over the timescale of decades.

5.5 Contrasting Ascent Dynamics

At present, multiple mechanisms have been proposed to explain the unprecedented deformation field at Uturuncu volcano. We have demonstrated how co-located uplift and subsidence at Uturuncu volcano may result from quasi-steady magmatic transport from the lower to middle crust (i.e., the “dipole model”). Mass balance arguments may explain inferred volume discrepancies from such a model, if magma compressibility or reservoir compressibility vary by one order of magnitude. The dynamics of transport between reservoirs have not been rigorously tested, but an appealing aspect of this model is that cessation of intrusions could lead to cooling and net depressurization in the middle crust, consistent with geomorphic evidence, for no net uplift on 100-year timescales, and perhaps even net subsidence on million-year timescales. However, the true rheology of the APMB and lower crust may be anelastic, even over yearly timescales. The dipole model could be extended to viscoelastic conditions in future work.

Alternatively, the model of viscoelastic diapiric ascent can reproduce the secular trend of uplift in geodetic data (Fialko and Pearse, 2012). The diapir model is appealing because it restricts observable magmatic activity to the depth of the APMB, eliminating concern over poorly constrained lower-crustal properties. However, the model predicts a secular slowing and broadening of uplift (without magmatic recharge to the APMB), and it is unclear how to explain transient accelerations in uplift rate. Further modeling could quantify the effect of distinct APMB geometries, recharge rates, and internal flow of partial melt on predicted surface displacements due to diapiric ascent. Finally, we remark that other hybrid models have been proposed that utilize the diapir-like geometry of the APMB, but we focus on viscoelastic pressurization rather than buoyant ascent (Gottsman et al., 2017).

Continuous GPS data near the summit of Uturuncu volcano show a change in the average rate of uplift between April 2010 and March 2015 compared to previous InSAR studies. InSAR observations between May 1992 and January 2011 predict an average vertical uplift rate of 7 ± 2 mm/yr at station UTUR. However, the GPS time series is better fit with a time-dependent function, such that the uplift rate in 2010 is less than 2 mm/yr and in 2015 may be as high as 9 mm/yr. We reexamined InSAR data and suggested a possible slowdown in the rate of deformation between 2003 and 2010. Variation levels of local seismicity at Uturuncu (three local events per day from April 2009 to April 2010; Jay et al., 2012) could be compared to variations in continuous GPS measurements in future work.

A finite-element model was used to evaluate effects of known mechanical heterogeneities underneath Uturuncu, such as the Altiplano-Puna magma body (APMB). An oblate weak zone such as the APMB disproportionately affects radial and vertical surface displacements, potentially biasing the inferred source geometry if homogeneity is assumed. Importantly, the superposition of multiple point sources impacts the ratio of radial to vertical displacements as well. In the case of Uturuncu, the homogeneous dipole source has a ur/uz ratio of 0.46, but accounting for the mechanically weak APMB reduces the ratio to 0.30 (recall that a single Mogi inflation source has a characteristic ratio of 0.38). We presented an inversion of four independent tracks of InSAR data to observationally constrain this ratio and found a best-fit solution of 0.33 but with considerable uncertainty (±0.15).

Finally, we have discussed several pros and cons of existing geodynamical models of magma transport at Uturuncu in light of new geodetic observations. Regular collection of continuous GPS data from the COLO and UTUR stations will provide context for the secular trend and annual variations in uplift at Uturuncu volcano observed between 2010 and 2015. Simultaneous, wide-footprint InSAR observations acquired monthly with the ESA Sentinel-1 constellation will be less prone to temporal aliasing compared to observations before 2010. The future of geodetic observations at Uturuncu is bright, promising much better temporal constraints going forward on measurements of active magmatic processes.

This work was supported by National Aeronautics and Space Administration (NASA) grants NNX08AT02G and NNX10AN57H issued through the Science Mission Directorate’s Earth Science Division. Additional funding came from National Science Foundation (NSF) grant EAR-0908281, which is part of the PLUTONS project ( Field support for GPS installation was provided by UNAVCO through the GAGE Facility with support from the NSF and NASA under NSF Cooperative Agreement No. EAR-1261833. We thank the Nevada Geodetic Laboratory for providing GPS time-series solutions. Comments from two anonymous reviewers and the editor Raymond Russo strengthened this manuscript. We thank members of the PLUTONS team for critical discussions, Kevin Ward for kindly providing the velocity model, and Yuri Fialko and Jillian Pearse for providing results from their independent InSAR analysis. Thanks as well to Sarah Doelger for support installing the GPS stations and Eric Kendrick for retrieving data from COLO and UTUR. We also benefitted from discussions with Michael Bevis and Julie Elliott regarding the GPS time series.

1Supplemental Material. Additional GPS time series and finite element modeling figures. Please visit or the full-text article on to view the Supplemental Material.
16 figures; 3 tables; 1 supplemental file

Supplementary data