The Altiplano-Puna Volcanic Complex of the Central Andes is host to an ~150-km-wide, quasi-circular ground deformation anomaly centered on Uturuncu volcano (Bolivia). The precise onset and duration of this deformation is unclear, but geomorphologic studies bracket its initiation at less than a few hundred years ago. Here we report on the deformation history over an ~50 yr period by deriving orthometric height changes from leveling and global navigation satellite system (GNSS) observations at 53 benchmarks along a regional leveling line that crosses the deformation anomaly. The comparison of interferometric synthetic aperture radar (InSAR) line-of-sight (LOS) displacements and LOS-projected orthometric ground velocities in a common reference frame reveal central uplift extending to ~35 km from Uturuncu at a maximum orthometric rate of 1.2 cm yr–1, and peripheral subsidence at a maximum rate of 0.3 cm yr–1 to ~60 km from Uturuncu. This pattern is consistent with the spatial extent and average rate of deformation observed by InSAR. Our interpretation of the data is that long-wavelength ground uplift at Uturuncu has likely occurred at a quasi-constant rate for at least half of a century. This study bridges the observational time spans between modern satellite geodetic observations (up to a few decades) and geomorphological observations (a few centuries and longer) of the recent deformation history of the continental crust in the Central Andes and adds to a select group of case studies of quantifiable long-term volcano deformation worldwide.

Satellite remote sensing techniques documented an unusually wide and spatially complex ground deformation in the Altiplano-Puna Volcanic Complex (de Silva, 1989) between A.D. 1992 and 2011 (e.g., Pritchard and Simons, 2002; Henderson and Pritchard, 2013; Fialko and Pearse, 2012) (Fig. 1), with a diameter of 150 km, a peripheral moat of subsidence, and a central area undergoing uplift at rate of up to 1 cm yr–1 at the western slopes of Uturuncu volcano (Bolivia; 22.270°S, 67.180°W) as detected in interferometric synthetic aperture radar (InSAR) line-of-sight (LOS) displacement data. The spatially and temporally complex deformation may be linked with an ~200-km-wide and ~11-km-thick geophysically anomalous upper-crustal magma reservoir—the Altiplano-Puna magma mush or magma body (Fig. 1) (Zandt et al., 2003; Fialko and Pearse, 2012; Comeau et al., 2016; Gottsmann et al., 2017). This reservoir is imaged with a volume of 300,000 km3 and is thought to contain up to 25 vol% of wet andesitic melt (Laumonier et al., 2017; Ward et al., 2014).

Situated in a large Neogene ignimbrite province where a total of >15,000 km3 of ignimbrites have been erupted over the past 15 m.y. (e.g., Salisbury et al., 2011; Kern et al., 2016), Uturuncu volcano ceased its eruptive activity ~250 k.y. ago (Muir et al., 2015; Sparks et al., 2008) and must hence be deemed extinct. However, this Pleistocene volcano shows signs of unrest with summit fumarolic activity and its western flank in the center of the uplift area. Therefore, questions remain as to the possibility of a future renewal of eruptive activity at Uturuncu volcano and the Altiplano-Puna Volcanic Complex as a whole, particularly in light of a substantial volume of melt stored within the Altiplano-Puna magma body.

While satellite remote sensing and ground-based observation only provide a temporally narrow window of opportunity to investigate this large-scale crustal deformation, geomorphic investigations have extended the window to time scales of hundreds to thousands of years. The main finding from these studies is that the observed deformation has been ongoing for no more than ~100 yr (Perkins et al., 2016), leading to the conclusion that the current deformation can be regarded as a result of transient changes in the crustal stress conditions beneath the Altiplano-Puna Volcanic Complex (Gottsmann et al., 2017). There is hence a disconnect between the annual time scale of global navigation satellite system (GNSS) observations, the decadal time scale of InSAR observations, and the time scale of geomorphic observations as to the longevity of this transient deformation and its time of onset. In contrast, a similar long-wavelength (~100-km-wide) surface deformation pattern (i.e., central uplift and peripheral subsidence) with a near-constant uplift of 0.2 cm yr–1 above the Socorro magma body in New Mexico, USA (Finnegan and Pritchard, 2009; Pearse and Fialko, 2010), has been documented geodetically for >100 yr and provides a rare glimpse into the time scales and rates of large-scale deformation from magmatic processes in the continental crust.

The purpose of this paper is to improve our understanding of the duration, wavelength, and magnitude of the ongoing deformation in the Altiplano-Puna Volcanic Complex by combining leveling data collected in 1965 with new GNSS survey data. We expect this research to inform future studies into the cause of the deformation, its spatio-temporal evolution, and likely future evolution.

Leveling Survey

We have retrieved data from the Bolivian geographical survey (Instituto Geográfico Militar de Bolivia) from one geodetic leveling line (termed BP), which were measured in 1965 west of Uturuncu in the Bolivian Altiplano region (Fig. 1). With a line length of 205 km, the original spacing of the benchmarks along BP was between 1 and 3 km. The line is classified as a first-order survey, implying that data were measured with an overall vertical accuracy of <7 cm, given by 0.5 cm × L1/2, for a class II survey, where L is the length of the survey in kilometers (Vanicek et al., 1980). Our initial field inspection in 2010 found that 53 original benchmarks along BP remained unaltered and undisturbed and fit for measuring. The rest of them were either slightly damaged, completely destroyed, or buried, i.e., deemed unfit for measuring. Figure 2 shows the spatial extent of the main BP line and its branch lines within the deformation anomaly detected in the InSAR data.

GNSS and Gravity Surveys

We obtained dual-frequency GNSS measurements in rapid-static survey mode at 53 healthy benchmarks along line BP in 2011–2012 (Table 1; Fig. 2). The first benchmark we measured is BP07, located to the north of Uturuncu volcano and a distance of 12.4 km from the start of the original BP line. From BP07, we measured a number of healthy benchmarks running south toward Uturuncu and along the road running west toward Laguna Colorada. A long segment of BP runs along a fringe of the InSAR-detected deformation anomaly (Fig. 2B), so little information on spatial deformation changes can be gleaned from that segment. We measured only two points in this segment (BP48 and BP57), and then continued with the more systematic measurements starting from BP69A with an average point spacing of 2 km. From there, the leveling line runs radially away from Uturuncu toward the southwest and reaches the border with Chile (BP100G), across and past the far side of the peripheral subsidence as identified by InSAR. Benchmark identifiers and coordinates are given in Table 1 and in the Supplemental File1.

We measured positions of benchmarks using dual-frequency Leica 1200 GNSS receivers and antennas. We used between three and six known GNSS bases located at a maximum distance of 100 km, recording between 20 min and several hours with a recording interval of 3–30 s, depending on the distance to the base stations used. Data were processed using Leica Geo Office and TopCon Tools software (;, and phase ambiguities were solved in every case, to achieve a vertical GNSS repeatability of <5 cm. We also took relative gravity measurements at most benchmarks using a Scintrex CG5 Autograv (serial number 572) as part of a larger survey (del Potro et al., 2013). For BP07 to BP14, BP16 to BP21, and BP35, where gravity was not measured on the benchmark, gravity values were extrapolated by a distance-weighted average of measurements done on nearby locations on either side of the BP benchmarks along the same road. These gravity values were reduced by taking into account the elevation difference between the benchmarks and a Bouguer slab density of 2270 kg m–3 (see del Potro et al., 2013). To set an absolute reference for the gravity measurements, we chose an arbitrary benchmark (BP15) where gravity was measured during the survey, and used the gravity value on the geoid using model EGM2008 (goBP15 = 9.7887 m s–2) to calculate the absolute gravity on the benchmark (gBP15 ) using the simplified Poincaré-Prey reduction (Hofmann-Wellenhof and Moritz, 2005) for an average crustal density of 2670 kg m–3:
where HBP15 is the orthometric height of BP15 (see below for calculation of orthometric heights). We use the so-derived value gBP15 = 9.7923 m s–2 as a reference for all relative gravity values measured in the survey. This procedure does not affect the gravity change values along the leveling line relative to BP15. See Table 2 for a list of notations used in this study.

Conversion to Orthometric Heights

GNSS data cannot be directly compared with leveled heights. The former provide Cartesian (XYZ) solutions that do not express a notion of elevation above a certain datum, whereas the latter depend on the geopotential surface along the path followed (Fig. 3). To permit direct comparison, both data sets must be converted to orthometric heights.

First, GNSS solutions are transformed into geodetic latitude, longitude, and ellipsoidal heights using the WGS84 ellipsoid model (Hofmann-Wellenhof and Moritz, 2005). Ellipsoidal heights (h) are then converted to (approximate) orthometric heights (H) by means of the geoid undulation (N; Figs. 3 and 4) (Hofmann-Wellenhof and Moritz, 2005) using:

For this we have used the national geoid model BOLGEO (Corchete et al., 2006) (Fig. 4 and the Supplemental File [footnote 1]), which provides better accuracy compared to a global model such as EGM2008. EGM2008 and BOLGEO can be reconciled by adding a vertical offset of 2.2 m to EGM2008.

Because of the nonparallelism of geopotential (level) surfaces, elevations of benchmarks from spirit leveling are not only different from heights above the geoid (such as orthometric heights) but also ambiguous; i.e., they are path dependent (Hofmann-Wellenhof and Moritz, 2005). Level heights (h′) must hence be converted to orthometric heights (H) (see Fig. 3). This is done by adding an orthometric correction (OC) to h′. Because we do not know the datum for the leveled heights, orthometric corrections are given relative to the previous point measured along the line. The correction between two benchmarks on the Earth’s surface, A and B, is given by (Hofmann-Wellenhof and Moritz, 2005):
where γ0 is the theoretical (normal) gravity calculated for 22.26°S latitude (9.784 m s–2), g is the gravity determined at the level points between A and B, and Δh′ is the leveling increment. hA and hB are the leveled heights of points A and B respectively (see Supplemental File [footnote 1]). forumla are the mean values of gravity along the plumb lines between the geoid and benchmarks A and B and are calculated as the arithmetic mean of g at the level points and their respective gravity on the geoid go (using a simplified Poincaré-Prey reduction) (Hofmann-Wellenhof and Moritz, 2005):

Because in most cases there were no gravity measurements taken between BP benchmarks, the summation of g was simplified to the average of gravity measured at two consecutive benchmarks.

The two sets of orthometric heights, derived from leveling and GNSS observations, are compatible and thus comparable (see Supplemental File [footnote 1]). Derived orthometric height changes between 1965 and 2012 along the main leveling line vary between −0.60 m and −1.15 m and are shown in Figure 4D. The constant offset is a combination of datum offsets, the choice of reference gravity, and the vertical movement of the entire region between 1965 and 2012. This implies a maximum relative height change along the leveling line of 0.55 m between the two surveys. We can derive the spatial deformation pattern from the orthometric height differences over 47 yr by plotting the differences along a radial westward trajectory from the center of the InSAR anomaly (Fig. 5).

The orthometric height change pattern shows relative uplift of up to 45 cm in the area of ground inflation detected by InSAR to a radial distance of ~40 km. From there, up to 15 cm of subsidence are observed over 20 km, which coincides with an area where InSAR data show the start of an ~20-km-wide area of LOS subsidence. The resultant orthometric deformation pattern is hence similar to the spatial deformation pattern derived from 18 yr of InSAR observations in terms of relative magnitude, wavelength, and position of areas of uplift and subsidence.

The absence of a common datum for the InSAR and orthometric deformation data introduces an unknown vertical offset between them. To address this, one possible solution is to assume that a given point has not moved at all, whereby ΔH = 0; one then uses that point as an arbitrary anchor between the two surveys. A different approach, followed here, is to use the known deformation trend from 18 yr of InSAR observations as a calibration to test whether the deformation patterns are indeed comparable and can be used to inform the deformation history since 1965. This calibration, however, has to be done with care, as both orthometric height change and InSAR LOS displacement are one-dimensional measurements of the true three-dimensional (3-D) displacement (Fig. 6). Changes in orthometric height (ΔH) represent the vertical component (uz) of the true three-component displacement vector, while InSAR LOS displacement data give the projection of the true displacement vector onto the satellite LOS.

Comparison of Orthometric Height Changes to InSAR-Observed Deformation

As shown in Figure 6, given a true displacement vector vA of benchmark BPA, composed of both upward (uz) and westward (ur) components (whereby uz = ur in the case shown), the change in orthometric height (ΔHA) does not equal the change in InSAR LOS (LOSInSAR). This is also true if ΔHA is projected onto the InSAR LOS (LOSΔHA). The mismatch between the orthometric deformation and InSAR-observed deformation becomes smaller for deformation dominated by vertical movement, but larger where horizontal deformation dominates. Moreover, two different 3-D displacement vectors, vA and vB measured at BPA and BPB, respectively, can yield the same magnitude of orthometric deformation ΔHA,B if projected onto the InSAR LOS (LOSΔHA,B) despite drastically different magnitudes of InSAR LOS displacement. This is due to the direction of the horizontal deformation component relative to the satellite’s LOS with an incidence angle of 22° for the available descending InSAR tracks (Henderson and Pritchard, 2017).

Projecting Orthometric Displacements onto InSAR LOS

The BP benchmarks located within the uplifting part of the InSAR anomaly (BP12 to BP37) are located northwest of Uturuncu. Assuming that the horizontal component of the deformation is radially away from Uturuncu, the east-west component of the deformation for a descending radar track that illuminates the area from the east would contribute to the lengthening of the radar path and yield a smaller magnitude in the InSAR LOS displacement compared with the magnitude of the projected orthometric height change (LOSΔH) (see Fig. 6).

We use the orthometric vertical displacements for each benchmark to project them into the InSAR LOS and derive a set of synthetic LOSΔH displacements that would arise from a perfectly radially symmetric uplift pattern such as one caused, e.g., by a buried pressurized sphere. Then we use the available InSAR LOS displacements data to help determine the magnitude of the vertical offset k in the orthometric height changes (ΔH) if projected onto the InSAR LOS, whereby:
in order to obtain LOS displacements in a common reference frame. Note that the magnitude of LOSΔH is modulated by parameter w whose value is dependent on the magnitude of radial displacement (ur) relative to the magnitude of vertical displacement (uz) at the leveling benchmarks. For example, for uruz = 0 (i.e., only vertical displacement and no radial displacement), w = 1 and k = LOSInSARLOSΔH. For uruz = 1 (i.e., ur = uz), w = 0.4 and k = LOSInSAR – 0.4 × LOSΔH.

In the above analysis, we limit the parameter space exploration of uruz to values up to 1, because a strong horizontal component acting radially from the center of the deformation would cause a horizontally elongated InSAR deformation pattern. This is inconsistent with the large-scale concentric InSAR pattern observed to date (e.g., Henderson and Pritchard, 2017).

Derivation of Deformation Velocities between 1965 and 2012

Using ur and uz data from an inversion of four overlapping InSAR tracks (Henderson and Pritchard, 2017), we find that uruz ratios are between −0.05 and 0.25 over a distance of ~60 km from the center of the InSAR uplift (Fig. 7A). According to Figure 7B, w hence takes values between 0.9 and 1.

The resultant vertical offset k from Equation 6 obtained by a residual minimization procedure is 2.0 ± 0.05 cm yr–1, or 94 ± 2.4 cm over 47 yr. Note, that this extrapolation assumes a constant linear evolution of the deformation between 1965 and 2012.

Adding the offset to the above LOS-projected orthometric height differences, we derive the spatial ground velocity pattern shown in Figure 8. This pattern is remarkably similar to the spatial deformation pattern from 18 yr of InSAR observations in terms of magnitudes and wavelengths of both inflation and subsidence. This implies that the uplift rate recorded along the leveling line running north-south toward Uturuncu (benchmarks BP12 to BP26) is equal to that recorded along the east-west segment (BP32 to BP45), suggesting that the vertical deformation is indeed radially symmetric from Uturuncu. A deformation source that causes a radially symmetric deformation pattern can hence reasonably explain the observations.

Extrapolating the LOSΔH ground velocities to the center of the uplift observed by InSAR and accounting for uncertainties in the derivation of velocities, we propose a maximum velocity LOSΔH of 1.1 ± 0.2 cm yr–1 . These values are within the uncertainty of maximum InSAR LOS ground velocities of 1.05 ± 0.07 cm yr–1. The orthometric uplift extends ~35 km from Uturuncu, and peripheral subsidence, at a maximum rate of ~0.3 ± 0.03 cm yr–1, extends to ~60 km from Uturuncu. The fact that the leveling line crosses the moat of subsidence to the west of the center of deformation, where InSAR observations are less reliable to detect the true magnitude of ground deflation, allows us to confirm that subsidence has indeed occurred albeit at slightly higher mean velocities than inferred from InSAR measurements.

The consistency of the spatial pattern and ground velocities between the orthometric and InSAR deformation data indicates that the long-wavelength deformation may have been constant for at least 47 yr (1965–2012). This allows us to suggest a near-constant deformation of the Altiplano-Puna Volcanic Complex over decadal time scales under the assumption that this is not a coincidence of the discrete sampling in time (1965, 1992, 2010–2011) of an oscillatory inflation-deflation process. Deriving mean velocities over such time scales from only two surveys is inherently problematic. Currently, we have no means to test our hypothesis due to the absence (to the best of our knowledge) of additional geodetic data over the same period. However, continuous GNSS observations near the center of uplift observed by InSAR indicate that the recent rate of uplift dropped to an average of 0.24 ± 0.19 cm yr–1 between 2010 and 2016 (Blewitt et al., 2016), similar to the long-term average detected above the Socorro magma body (Finnegan and Pritchard, 2009). Whether this uplift will continue to decelerate or accelerate again in the future remains to be seen. Either behavior has implications for the assessment of causative processes behind the complex spatio-temporal deformation. Their discussion is beyond the scope of this paper, and we refer interested readers to Gottsmann et al. (2017) where some of the implications of nonlinear history of the deformation anomaly are considered.

Leveling and GNSS data collected in 1965 and 2012, respectively, have allowed us to convert deformation data into orthometric height changes along a leveling line that crosses the 150-km-wide deformation anomaly centered on Uturuncu. Due to the absence of a common geodetic datum for the orthometric and InSAR displacement data, the two data sets cannot be directly compared. However, we obtain a reasonable match between the two data sets by projecting the orthometric height changes into the InSAR LOS and applying a vertical offset of 94.0 ± 2.4 cm over the 47 yr to the orthometric height changes. The ground deformation velocity pattern from 1965 to 2012 is equivalent to that observed by InSAR surveys between 1992 and 2011. We therefore propose a near-constant long-term pattern of ground deformation between 1965 and 2012 with central maximum orthometric uplift at 1.2 ± 0.2 cm yr–1 and maximum peripheral subsidence at 0.3 ± 0.03 cm yr–1 over the 150-km-wide deformation anomaly. This study adds to a select group of long-term geodetic investigations that quantitatively document protracted periods of complex spatio-temporal deformation in volcanic terrain (e.g., Finnegan and Pritchard, 2009; Parker et al., 2014). Exploitation of these data should improve our understanding of the link between protracted crustal stress changes and decadal time scales of magma transfer and reservoir growth (e.g., Druitt et al., 2012).

This work was funded by the Natural Environmental Research Council (NERC, grant NE/G01843X/1), two NERC-Geophysical Equipment Facility instrument loans (910 and 928), and the EC-FP7 VUELCO project (grant agreement 282759). M. Sunagua, SERGEOTECMIN, and SERNAP provided assistance in the field in Bolivia. The paper benefitted from thoughtful comments by two anonymous reviewers and guest associate editor Matthew Pritchard.

1Supplemental File. Benchmark locations, geoid undulations, and gravity data (text file). Please visit or the full-text article on to view the Supplemental File.
Science Editor: Raymond M. Russo
Guest Associate Editor: Matthew Pritchard
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data