Abstract

Geodetic surveys of Volcán Uturuncu and the Lazufre volcanic complex in the Central Andes of South America reveal sustained surface uplift from magmatic intrusion at depth. However, the decadal timescales of geodetic surveys are short relative to the timescales of magma chamber growth. Thus, from geodesy alone, it is difficult to infer the deformation and hence magma accumulation history of these volcanoes. Here we combine data from InSAR, long-wavelength topography, GPS and high-resolution topographic surveys of lake shorelines and rivers, and lava flow morphology to constrain the spatial and temporal evolution of magmatism at Uturuncu and Lazufre. Near Uturuncu, dated lake shorelines show no evidence of tilting since ca. 16 ka, and we find no evidence of deformation in the long-wavelength topography. A lack of net surface displacement suggests that uplift related to a rising diapir must be less than a century old, or, more likely, magmatic inflation at Uturuncu is transient over millennial timescales and is therefore not recorded in the topography. At Lazufre, we also find no evidence for sustained uplift recorded in Late Pleistocene lake shorelines. However, the orientations of multiple dated lava flows suggest that the long-wavelength dome at the center of Lazufre’s uplift has persisted since at least 400 ka. Additionally, we find that the radial distribution of volcanic vents at Lazufre, coupled with the presence of an apical graben, is consistent with experimental and theoretical predictions of magmatic doming. The dome’s longevity indicates significant magma storage at depth, and therefore Lazufre is likely a highly evolved pre-caldera magmatic system. These two case studies demonstrate that combining geomorphic and geophysical data sets to extend the geodetic record back in time can help determine the style and magnitude of magma transport in volcanic systems.

1. INTRODUCTION

The relatively recent discovery of actively uplifting volcanoes in the Central Volcanic Zone of the Andes (14°–28°S) using interferometric synthetic aperture radar (InSAR) geodesy (Pritchard and Simons, 2002, 2004) has spurred much research devoted to monitoring surface deformation and understanding the magmatic inflation processes that are thought to be driving these systems (e.g., Singer et al., 2014). The uplift patterns of Volcán Uturuncu in southern Bolivia (22°16′S, 67°11′W) and the Lazufre volcanic complex on the border of Chile and Argentina (25°18′S, 68°26′W) are particularly noteworthy because the footprint of their deformation is on the spatial scale of large caldera systems in the Central Andes (Fig. 1) (Pritchard and Simons, 2002; Ruch et al., 2008; Sparks et al., 2008).

Large silicic magma chambers tend to grow and evolve over much longer timescales (104–105 yr) than can be captured by decade-scale geodesy (Lindsay et al., 2001; Vazquez and Reid, 2002), and therefore the picture of magmatic evolution gleaned solely from satellite data is often incomplete. Geomorphic processes operate over timescales similar to magma chamber growth, and thus topography has the potential to offer additional constraints on the growth history of large volcanic systems. Tectonic geomorphology as a field has developed many tools devoted to tracking the paleogeodetic history of landscapes (e.g., Burbank and Anderson, 2011), and lake shorelines and river deposits and their characteristic morphologies have proven to be effective strain markers in tracking surface deformation from magmatism specifically (e.g., Pierce et al., 2002). Additionally, the large spatial and temporal scales associated with magmatic systems imply that signatures of deformation may even appear in long-wavelength topography (e.g., Froger et al., 2007; Gregg et al., 2012). Figure 1 shows that many volcanic complexes in the Central Andes correlate with long-wavelength topographic anomalies. Notably, Pastos Grandes caldera, the Putana and Purico complexes, Cerro Panizos, Coranzuli, Cordón Puntas Negras, El Queva, Llullaillajo, Lastarria-Cordón del Azufre, the Antofalla volcanic complex, Cerro Galan, and Nevados Ojos del Salado all have long-wavelength anomalies that are suggestive of broad magmatic uplift (e.g., Lipman, 1984).

In this manuscript, we focus on the two primary actively uplifting volcanic centers in this region—Volcán Uturuncu and Lastarria-Cordón del Azufre (also known as Lazufre). Interestingly, the 70-km-diameter uplift field at Uturuncu appears to be centered on the edifice of the ∼15-km-wide volcano with no underlying long-wavelength anomaly (Fig. 2A), while the 60–70 km deformation field at Lazufre is correlated with a long-wavelength dome yet has no central stratovolcano (Fig. 2B). This simple observation suggests either that Lazufre and Uturuncu are at differing stages in their evolution or that they are operating by different magmatic processes entirely. Here we integrate data from satellite geodesy, paleogeodetic measurements from the geomorphology of lakes, rivers, and lava flows, and bulk topographic characteristics to infer the spatial and temporal scales of surface uplift at Uturuncu and Lazufre. In doing so, we hope to provide observational constraints that allow differentiation between end-member models of how magma is being transported, stored, and potentially erupted at these sites. Below we describe the background of each case study and the hypotheses that our methods will be able to test.

1.1. Volcán Uturuncu

Uturuncu is a 6008-m-high, Pleistocene stratovolcano located in southern Bolivia within the Altiplano-Puna volcanic complex (APVC), an extensive region of silicic calderas and voluminous ignimbrite sheets from 21°S to 24°S (de Silva, 1989). It was primarily active from ca. 1050 ka to ca. 250 ka and has yielded ∼50 km3 of high-K dacite and silicic andesite from effusive eruptions over its lifetime (Sparks et al., 2008; Muir et al., 2015). Recent and rapid uplift at Uturuncu was first discovered through an InSAR investigation by Pritchard and Simons (2002), who documented a peak uplift of >1 cm/yr since the onset of InSAR measurements in the region in 1992. Estimated source depths from inversion of the geodetic data suggest that the deformation could be coming from magma inflation at the top of a large regional zone of partial melt known as the Altiplano-Puna magma body (APMB), which is thought to be the magmatic source for the flare-up of ignimbrite volcanism in the APVC (Zandt et al., 2003). Recent swarms of seismicity are additionally thought to reflect pressurization of upper crustal rocks by magma inflation at depth (Jay et al., 2012).

More recent InSAR investigations using multiple look angles (Fialko and Pearse, 2012; Henderson and Pritchard, 2013) have shown evidence of a deflationary moat along the outer edges of the uplift field at Uturuncu. Two primary end-member models have been put forward to account for simultaneous inflation and peripheral deflation: lateral migration of melt and growth of a diapir from the APMB (Fialko and Pearse, 2012); or upward migration of melt out of a deeper (55–80 km) reservoir and into to the top of the APMB (Henderson and Pritchard, 2013). Both gravity (del Potro et al., 2013) and magnetotelluric (Comeau et al., 2015) surveys show a low-density, melt-rich body at the top of the APMB; however, the lack of spatial resolution of these techniques at depths below the APMB prevent them from being able to uniquely identify the mode of magma transport at Uturuncu.

Integrated over the 104–105 yr timescales that are relevant to the growth of large silicic magma chambers, however, these two models yield different predictions of long-term surface uplift. Because diapirism is a buoyancy-driven viscous process and should be relatively continuous in time (e.g., Burov et al., 2003), we expect there to be a record of long-term surface uplift at Uturuncu. Alternatively, if Uturuncu is experiencing transient pulses of inflation from magma flux between shallow and deep reservoirs, separated by periods of relative quiescence, then we may not expect to see a record of long-term net surface uplift (e.g., Pierce et al., 2002). The systems of rivers and Pleistocene lakes that flank Uturuncu (Fig. 3) provide constraints on surface deformation going back tens of thousands of years; therefore we may be able to use the rich paleogeodetic record to unravel the long-term relationship between magmatism and surface deformation here.

1.2. Lazufre

The name “Lazufre” is a portmanteau that joins together the names of Lastarria volcano (25°10′S 68°30′W) and Cordón del Azufre (25°20′S 68°31′W), two ca. 0.6–0.3 Ma volcanoes on the periphery of a rapidly uplifting ∼60-km-diameter elliptical zone that has accelerated to 2.5 cm/yr since 2005 (Naranjo and Cornejo, 1992; Froger et al., 2007; Ruch et al., 2008; Henderson and Pritchard, 2013). On the basis of InSAR inversions, the source of deformation is thought to be either an inflating sill or a flat-topped magma chamber at ∼10 km depth; however, the source depth is poorly constrained and ranges anywhere from 2 to 14 km depth (Pritchard and Simons, 2004; Froger et al., 2007; Ruch et al., 2008; Anderssohn et al., 2009; Pearse and Lundgren, 2013; Remy et al., 2014). Recent seismic investigations by Spica et al. (2015) find no evidence for a regional chamber at Lazufre above 7 km depth, which constrains this range effectively from 8 to 14 km. Observations of regional vent orientations and faults suggest the elliptical shape of Lazufre may be due to a regional maximum horizontal compressive stress oriented NW-SE (Ruch and Walter, 2010).

The uplift field at Lazufre is not centered on any particular volcano but rather on a long-wavelength, 70-km-diameter and ∼500-m-high topographic dome (Froger et al., 2007; Figs. 1 and 2), with an apical graben at its center (e.g., Pritchard and Simons, 2002). The dome is flanked by a series of smaller Pleistocene volcanoes, all of which are hypothesized to be part of a larger pre-caldera system (Froger et al., 2007; Spica et al., 2015). Implicit in this hypothesis is the idea that the magma reservoir beneath Lazufre was active during the formation of the peripheral ring volcanoes, and therefore a key test is whether these features can be genetically related to one another. A primary goal of this study is then to constrain the rate and timing of dome growth at Lazufre, as well as the structural linkages between the dome and the proximal volcanoes, using a variety of geomorphic tools in order to better understand its potential as a rapidly inflating pre-caldera system.

2. METHODS

2.1. High-Resolution DEM Generation

Mapping lake-shoreline deposits requires a digital elevation model (DEM) of sufficient resolution to resolve the subtle slope breaks that characterize these features. In areas such as the Central Andes, where LiDAR data are difficult to obtain and abundant vegetation is not an issue, satellite photogrammetry can produce digital elevation data of comparable resolution for a significant fraction of the cost. We acquired 0.5-m-resolution, stereo-paired Worldview-1 satellite imagery for a 35 × 35 km zone containing Uturuncu and its surrounding area and produced a DEM using the open-source StereoPipeline software developed by the National Aeronautics and Space Administration (NASA) (Moratto et al., 2010). The output resolution of the DEM is also 0.5 m; however, the data are more likely accurate to ∼2 m resolution. To produce an absolute DEM, we corrected relative elevations using 90 m Shuttle Radar Topography Mission (SRTM) data. To validate the DEM elevations, we then compared them to high-precision campaign GPS measurements (del Potro et al., 2013) as well as our own kinematic differential global positioning system (dGPS) surveys (Fig. 4). There is generally very good agreement between the campaign GPS data and the Worldview DEM, with a root-mean-square error (RMSE) of 1.4 m including the kinematic dGPS surveys, and 1.0 m without (Fig. 4B). Figure 4C shows a zoomed-in comparison between the kinematic GPS data and the DEM data that indicate a slight bias toward higher GPS elevations at high SRTM-corrected elevations, as well as some noise that produces horizontal excursions on the 1:1 line. These excursions are due to the fact that the photogrammetry algorithms have difficulty correlating imagery with steep slopes, which leads to the generation of noise toward the summit of Uturuncu as well as the steep-sided boundaries of Laguna Mama Khumu.

2.2. Topographic Analysis

We use three different methods of topographic analysis to measure surface deformation in response to magmatic uplift at Uturuncu and Lazufre. Because the modern uplift fields are relatively long wavelength, we examine the spatial correlation between uplift and various wavelengths of topography (e.g., Froger et al., 2007). Additionally, we utilize the well-preserved erosional and depositional landforms produced by both lakes and rivers within our study areas as topographic strain markers. Below we describe our topographic and geomorphic approaches to constrain long-term surface uplift.

2.2.1. Radial Topographic Profiles

A convenient way to visually inspect the spatial relationship between a radially symmetric uplift field and radially symmetric topography is to view the data as a function of radial distance from the center of the anomaly. In Figure 2, we plot radial line-of-sight (LOS) velocity profiles outward from the local maxima of uplift at both Uturuncu and Lazufre (e.g., Froger et al., 2007).We then compare the radial LOS velocity profiles to profiles of both the long- and short-wavelength components of topography. We do this by taking the 2D Fourier transform of a 90 m SRTM DEM using Matlab’s fft2 command and applying a low-pass filter with a cutoff wavelength of 35 km. We choose 35 km as a cutoff wavelength because it is larger than the wavelength of the vents at Lazufre as well as Uturuncu and smaller than the wavelength of the modern uplift field (∼70 km) at both study areas. Given the strikingly high degree of topographic preservation in these landscapes, a first-order hypothesis is that if the modern uplift field has been consistent over the 100 k.y. timescales that these volcanoes have been active, then there should exist a spatial correlation between the LOS velocity field and the long-wavelength topography. As Figure 2 shows, and as we discuss in the Results section, long-wavelength topography appears to be correlated at Lazufre and de-correlated at Uturuncu, suggesting distinct magmatic uplift histories.

2.2.2. Lakes as Tiltmeters

The bell-shaped uplift patterns at Uturuncu and Lazufre produce uplift gradients radially away from their peaks (e.g., Figs. 1–3). Thus, any planar surface subjected to this forcing will be tilted outward over time. Lake shorelines make excellent strain markers because they represent equipotential surfaces and are thus horizontal at the time of their formation. Using both field surveys and high-resolution topographic data from the DEM we produced, we map shorelines along the lakes within the modern uplift field at Uturuncu and Lazufre (Figs. 4–8). We then radially project shoreline elevations outward from the peak in LOS velocity.

The high resolution of our photogrammetry-derived Worldview DEM allows for remote mapping of paleoshoreline features (Fig. 7). The subtle break in slope that characterizes a shoreline inner edge (a proxy for lake level during the period of shoreline cutting) becomes apparent in a slope raster of the elevation data (e.g., Figs. 7 and 8). We use the standard slope algorithm in ArcGIS to calculate the local 2D gradient of the topography. We then map inner edges directly on the slope map for Laguna Celeste, Laguna Mama Khumu, and Laguna Chojllas (Figs. 3 and 7).

During our field campaign, we also surveyed lake-shoreline elevations using a kinematic differential GPS (dGPS) with Trimble Terrasync software. Given its favorable shoreline orientations relative to the active uplift field and the clarity of its geomorphology, we focused our survey efforts at Uturuncu along Laguna Mama Khumu (Fig. 3). Shoreline inner edges and delta fronts (another proxy for water level, e.g., Gilbert, 1890) are present at both the proximal and distal edges of the lake, allowing us to walk along these features with the dGPS (Fig. 5). To allow for direct comparison between the field and remote surveys, we extract DEM elevations from the field-mapped shoreline traces at Laguna Mama Khumu (Fig. 7B).

At Lazufre, only one lake (Laguna de la Azufrera) is proximal to the modern uplift field; however, its position at the outer edge of the uplift field means it is likely not experiencing a significant amount of deformation. Additionally, the lake is oriented unfavorably with regard to the uplift field (Fig. 6). However, it is located within the bounds of the long-wavelength topographic anomaly (Figs. 1 and 6). Much like the lakes surrounding Uturuncu, Laguna de la Azufrera has two prominent paleoshorelines (Fig. 9). Additionally, the east shore of Azufrera has a well-expressed beach ridge that runs parallel to the uplift gradient that we surveyed along its entire length in the field using kinematic differential GPS (Fig. 6B).

2.2.3. Rivers

River longitudinal profiles can also serve as long-wavelength recorders of surface deformation. At Uturuncu, three prominent tributaries of the Río Grande de Lípez (Fig. 3) flow northward past the volcano into the Uyuni basin. Valley bottoms along these drainages are typically composed of fluvial gravels, and their margins are often confined by near-vertical bedrock walls (typically volcanic deposits) (e.g., Fig. 10). Presently, however, the floodplains are typically composed of wetlands with bunch grasses and small, anastomosing channels (Fig. 10). The presence of wetlands and the small size of streams suggest that rivers are currently underfit relative to their valleys around Uturuncu. River channel change in gravel-bedded rivers is driven by the divergence in sediment transport capacity downstream, which itself is governed by channel slope. Hence, tilting of alluvial channels can drive incision or aggradation depending on the orientation of the channel relative to the tilt (Ouchi, 1985). The apparent inability of the rivers around Uturuncu to transport their underlying gravels (as indicated by wetlands at present in these channels) implies that the channels are likely not dynamically responding to any surface perturbations. In the Atacama basin to the west of Uturuncu, Nester et al. (2007) dated buried woody debris along fluvial terraces to show that material was last being transported at this location during the Tauca highstand at ca. 16–13 ka. If channels at Uturuncu have been dormant for an equivalent length of time, then the longitudinal river profiles of the Rio Quetena tributaries may be passively responding to surface uplift over this timescale without readjusting their profiles (and hence erasing the signature of deformation). Therefore, we hypothesize that their longitudinal profiles should directly reflect the uplift field and provide valuable strain indicators over this time period.

Longitudinal profiles along Uturuncu were mapped by hand in ArcGIS using our high-resolution Worldview DEM. Because the presence of wetlands within channels produces many anastomosing streams and frequently obscures channel thalwegs, the valley axis was used as a best approximation of the thalweg. Both elevation and InSAR LOS velocity data were extracted at points along the valley axis profiles. Elevation data were extracted at every meter, and InSAR data were interpolated to the smaller grid size of the DEM.

2.2.4. Mapping Lava Flow Deflections and Volcanic Vent Distributions at Lazufre

At Lazufre, the high degree of internal drainages created by overlapping lava flows coupled with the extremely arid conditions make for a conspicuous lack of throughgoing river networks. Therefore, we omit any analysis of drainage networks surrounding Lazufre and instead use the quaternary lava flows themselves to constrain the timing of uplift. Lava flows at Lazufre appear to radiate outward from the peak of the long-wavelength dome (e.g., Fig. 11), which leads to a simple hypothesis. If the lava flow directions match the aspect of the long-wavelength topography, then the topography was likely present before the flows erupted. If the direction of the flows do not correlate with the long-wavelength topography, however, that suggests the doming likely occurred after they erupted. This hypothesis is complicated by the fact that short-wavelength features such as individual volcanoes and neighboring lava flows can influence flow direction; therefore, we map as many flows as possible in order to mitigate this noise potential. We map flows within our study area from high-resolution satellite imagery in Google Earth and ArcGIS (see Supplemental File [.kmz]1). We overlaid a grid of 10 km × 10 km boxes, and chose to map only flows within each box that are longer than ∼2 km in length and that have obvious flow direction indicators. This includes flow ridges, levees, and lobe propagations. We delineate flow directions on Figure 11 from flow ridges only, and because flow ridges form perpendicular to flow direction (Fink, 1980), we trace lines orthogonally to each observed ridge train. For each mapped lava flow, we take the median direction of each mapped ridge train within the deposit. For each flow polygon, we then subtract the median flow direction from the median long-wavelength topographic aspect at each polygon location. This angular deviation then provides a measure of how close the flow direction mimics the background topography, where a low value is a close fit and a value of 180° would be the complete opposite direction.

Additionally, we map volcanic vent locations throughout Lazufre (Fig. 11). Volcanic vents provide information on magmatic pathways (and thus fractures) in the upper crust, and they may be used to infer the state of stress of the underlying magma chamber in pre-caldera systems (e.g., Martí et al., 1994; Ruch and Walter, 2010). In particular, the analog experiments of Walter and Troll (2001) show that, in inflating systems, an apical graben develops toward the center of the uplift, followed by outwardly radial fractures and circumferential thrust faults. If the dome at Lazufre is forming from pressurization by an underlying magma chamber, then its distribution of volcanic vents should reflect an underlying fracture distribution consistent with that stress field. Using the same imagery and 10 km × 10 km grid as described above, in addition to vent locations as mapped by Naranjo (2010), Naranjo et al. (2013) and de Silva and Francis (1991), we identify vents as circular depressions or points that distinctly source a lava flow, cone, or dome, and map their distribution in Figures 11 and 12. We thus map multiple apparent vent sources for individual volcanic edifices, since we are interested in using vents as a proxy for fractures in the crust. A .kmz file of our mapped vent locations and lava flows is available in our Supplemental File (see footnote 1).

2.3. Geochronology

Calculating rates of surface deformation using lake shorelines requires knowledge of the shoreline ages in addition to their displacement. Although the lake level history of the Uyuni and Poopo basins to the north of Uturuncu is well understood (e.g., Placzek et al., 2006), given the local hydrologic conditions and glacial history of Uturuncu (Blard et al., 2014), it is unclear how the large lake level changes in the Altiplano translate to these smaller, high-elevation basins.

To constrain the ages of prominent paleoshorelines along lakes that flank Uturuncu, we dated two well-preserved beach terraces on the northern shore of Laguna Loromayu using optically stimulated luminescence (OSL) (Fig. 8). Within each beach terrace deposit, we dug ∼1-m-deep pits and took two unique samples from subhorizontal beds of coarse to fine sand at different depths (Fig. 8C and Table 1). Samples were processed at the Utah State University Luminescence Laboratory to extract quartz sand for analysis following the single-aliquot regenerative dose protocol (e.g., Murray and Wintle, 2000) and calculated using the early background subtraction method of Cunningham and Wallinga (2010) and the central age model of Galbraith and Roberts (2012). Samples for environmental dose rate determination were analyzed for radioisotope concentration using inductively coupled plasma–mass spectrometry (ICP-MS) techniques (Table 1) and converted to dose rate following the conversion factors of Guérin et al. (2011). Sediments collected in air-tight containers were dried to determine in situ water content. Contribution of cosmic radiation to the dose rate was calculated using sample depth, elevation, and latitude/longitude following Prescott and Hutton (1994). Total dose rates were calculated based on water content, radioisotope concentration, and cosmic contribution (Aitken, 1998).

3. RESULTS

3.1. Uturuncu

3.1.1. Shoreline OSL Ages at Laguna Loromayu

The OSL ages for the two prominent beach terraces at Laguna Loromayu range from 13 to 17 ka, consistent with the Tauca highstand seen in the Poopo and Uyuni basins of the Altiplano (e.g., Placzek et al., 2006) (Fig. 13 and Table 1). This timing also corresponds to a period of rapid glacial retreat signaled by a rising equilibrium line altitude (ELA) at Uturuncu (Blard et al., 2014) (Fig. 13). Interestingly, the higher terrace at 4730 m has a younger combined age (14.3 ± 2.8 ka, n = 48) than the lower 4715 m terrace (16.6 ± 3.3 ka, n = 40). Although these ages are within error of each other, if indeed the lower terrace is geomorphically older than the upper terrace, the rise and fall of lake levels must have occurred rapidly in order to preserve the lower shoreline. This is consistent with the rapid shoreline ascent and descent seen during the Tauca highstand in the Uyuni basin (Placzek et al., 2006) (Fig. 13).

3.1.2. Deformation Timescales at Uturuncu from Lake Surveys

At Uturuncu, none of the surveyed shorelines show evidence for substantial, sustained tilting. Figure 7 shows mapped shorelines and radial projections from Laguna Celeste, Laguna Mama Khumu, and Laguna Chojllas. We omit shoreline surveys at Laguna Loromayu because it is located along the seam of smaller DEM tiles that were mosaicked together and thus has a slight elevation discrepancy that makes shoreline correlation more difficult. The pink lines within the lower panels of Figure 7 represent the amount of tilting expected at each lake assuming the modern uplift rate has been consistent since the Tauca highstand (ca. 16 ka).

Table 2 shows the average measured slopes in the direction of tilting and for shorelines at Laguna Celeste, Laguna Mama Khumu, and Laguna Chojllas. On average, there is a very small amount of observed tilting away from the summit. Comparing the tilt of measured shoreline profiles to the gradient of the surface uplift along the lakes allows us to estimate the longevity of the modern uplift at Uturuncu. Combining all our shoreline measurements, we calculate an average total tilt (in units of slope) for our measured profiles of ∼10–5 (m/m) away from the summit of Uturuncu (Table 2). The average tilt rate calculated from InSAR across the lakes is ∼10–7 m/m/yr; so dividing the measured shoreline tilt by the average geodetically observed tilt rate yields a rough estimate of ∼100 years for the persistence of the present uplift episode. This timescale also assumes a constant rate of uplift back in time, which may not be accurate. Regardless, since only two of the three lakes actually tilt away from the summit of Uturuncu (Mama Khumu appears to tilt negligibly in the other direction), the time constraints this estimate provides should be taken with a concordant amount of uncertainty.

We can also place more conservative upper bounds on our timing constraints by estimating the inherent uncertainty in measuring shoreline elevations. Here we use our GPS surveys of lake shorelines to estimate their vertical geomorphic uncertainty, and we find that on average a strandline elevation has a standard deviation of ±0.5 m. It is therefore possible to estimate the amount of time it would take for a clear deformation signal to be measurable above this uncertainty. For example, given the modern tilt rate of –3.2 × 10–7 yr–1 over the 3.9 km distance of Laguna Mama Khumu, it would take ∼800 years to achieve a 1 m differential elevation of shorelines on either side of the lake. Therefore, from this order of magnitude analysis of the topography and survey data, we estimate the age of the modern uplift episode at Uturuncu to likely be less than 100 years and almost certainly less than 1000 years.

3.1.3. River Profiles along Uturuncu

The profiles of alluvial rivers that circumnavigate the periphery of Uturuncu do not correlate with any gradients in surface uplift (Fig. 14). The best candidate for measuring changes in surface slope related to deformation is the mainstem of the Río Grande de Lipez (bold lines in Fig. 14). Here the uplift gradient linearly decreases from km 13 to km 36 along the mainstem, resulting in a tilt rate of ∼1.1 × 10–7 yr–1. Integrating the present tilting rates forward from the Tauca highstand would result in a doubling of surface slope along this section; however, the slope of the Río Grande de Lipez remains unchanged through the uplift field.

Along each of the tributaries, there are also two prominent knickpoints, which can often result from either spatial gradients in surface uplift or impulsive changes in uplift rate, but these do not appear to be correlated with an LOS velocity anomaly (Fig. 14). One knickpoint occurs at the junction of the northern and southern tributaries, which commonly results from a change in sediment load between the two rivers (e.g., Finnegan and Pritchard, 2009). The large knickpoint along the northern flank of Uturuncu is located at the northwestern base of Cerro San Antonio where the river appears to be diverting around the volcano (Fig. 3). Cerro San Antonio’s heavily dissected morphology indicates that it is significantly older than Uturuncu (ca. 3 Ma, de Silva and Francis, 1991); thus, this knickpoint is very likely unrelated to the magmatic processes beneath Uturuncu.

3.2. Lazufre

3.2.1. Constraints on Deformation at Lazufre

At Lazufre, Laguna de la Azufrera is at the distal margins of the observed uplift but within the bounds of the long-wavelength topographic anomaly (Figs. 1 and 6). However, like the lakes that surround Uturuncu, its shorelines are remarkably flat. Figure 8 shows survey data collected along the shorelines of Laguna de la Azufrera and a radial projection of their elevations away from the peak LOS velocity anomaly. In this projection, the upper shoreline has a negative slope (tilting away from the uplift) of ∼10–4; however, the lower shoreline that we surveyed for over a km has a positive slope (tilting down toward the uplift) of 10–4. Because the mean slope of the differing shorelines is nearly flat and averages approximately zero, there does not appear to be any signal of long-term deformation in the shorelines preserved at Laguna de la Azufrera. The horizontal lake shorelines are superimposed on the long-wavelength topographic bulge, which has a topographic slope of 3% across the lake (Fig. 6). This implies that the majority of surface uplift across the dome was accomplished before the highstands at Laguna de la Azufrera, which are likely Tauca in age.

3.2.2. Uplift Constraints from Lava Flow Deflections at Lazufre

Mapped lava flows in the vicinity of Lazufre, particularly along the western and southern sectors of the complex, tend to radiate outward from the peak long-wavelength topographic high and modern uplift center. Figure 11 shows a map of lava flow directions, vent locations, and line-of-sight velocity measurements at Lazufre. Although this type of analysis is inherently noisy—because short-wavelength topography like individual volcanoes can greatly influence the local direction of lava flows—most of the flow orientations are down gradient of the long-wavelength topography (Fig.11). This suggests that at the time of eruption for the quaternary volcanoes proximal to Lazufre, there was already a topographic high roughly where we find one today. There are two notable exceptions where lava flow directions clearly do not match the direction of the long wavelength topography. The flows associated with the apical graben, which is too small to be reflected in the long-wavelength topography, are inwardly dipping toward the center of the dome (Fig. 11). Additionally, the flows along the eastern flank of the dome, which are morphologically older than those to the west (de Silva and Francis, 1991), do not track with the long-wavelength topography.

Unfortunately there is not excellent age control for most of these lava flows; but constraints exist at Lastarria (LS), Volcan La Moyra (LM) at the northern segment of Cordon del Azufre (CA), and Cerro Bayo (CB), from K-Ar ages (Naranjo and Cornejo, 1992; Naranjo et al., 2013). The oldest deflected flows at Lastarria range from 0.6 ± 0.3 Ma to ca. 40 ka, 0.3 ± 0.3 Ma at Volcan La Moyra, and 1.6 ± 0.4 Ma at Cerro Bayo (Naranjo and Cornejo, 1992; Naranjo et al., 2013). The best constraint comes from a complex of elongated flows at Volcán Negriales (VN), the oldest of which is 400 ka (Naranjo et al., 2013). Because the longest flows at Lastarria and La Moyra appear to reflect the long-wavelength topography, it is reasonable to suggest that there was a topographic high at Lazufre since at least before ca. 0.4 Ma. The morphologically older (e.g., de Silva and Francis, 1991) volcanic edifices along the eastern edge of Lazufre do not appear to reflect the topographic forcing of the present dome; therefore, they likely pre-date the surface uplift. The youngest dated flow of these older structures is 0.9 ± 0.3 Ma from the southernmost flow of Volcan Río Grande (RG) (Richards and Villeneuve, 2002). Given these constraints, we can place approximate bounds on the onset of tumescence from 0.9 to 0.4 Ma. However, better geochronology is needed in order to test these constraints with any degree of certainty.

A >0.4 Ma magma chamber beneath Lazufre is consistent with both structural and volcanological evidence for a long-lived stress orientation similar to that observed today (Ruch and Walter, 2010; Spica et al., 2015) as well as magnetotelluric imaging of a partially molten body emanating from the upper mantle into the mid-crust (Budach et al., 2013). Thus, the magma body beneath Lazufre has possibly been causing punctuated surface inflation for hundreds of thousands of years.

3.2.3. Vent Distributions at Lazufre

In addition to mapping lava flow orientations, we also mapped visible vent locations throughout Lazufre using the high-resolution digital imagery within Google Earth (Fig. 11 and .kmz in the Supplemental File [see footnote 1]). This mapping complements previous work looking at the spatial distribution of vents within the area (e.g., Froger et al., 2007; Ruch and Walter, 2010) and provides additional detail from the high-resolution imagery. A key prediction from numerical and analog experiments of magmatic uplift of an inflating dome is a stress field that results in radial extensional fractures emanating from the summit of the dome and increasingly divergent circumferential reverse faults toward the apex that meet at approximately right angles (e.g., Komuro et al., 1984; Gudmundsson et al., 1997; Walter and Troll, 2001). We show an example from one such analog experiment by Walter and Troll in Figure 11B. Specifically to Lazufre, this pattern of fractures is predicted from the work of Spica et al. (2015), who find that a deep inflationary source could generate such a stress field and consequently lead to the formation of the smaller-scale magmatic anomaly at Lastarria volcano.

Visually, we see in Figure 11 that the distribution of vent locations highlights a pattern of radial cracks emanating from the peak of the dome and/or deformation source, an outer ring of vents at the margins of the uplift field, and a small inner ring that corresponds to the location of the apical graben (highlighted by the inwardly oriented lava flows) (Fig. 11). We attempt to quantify this spatial variability of vent distributions by plotting them as a function of radial and angular distance from the peak of the uplift field (Fig. 12). Visualized in this radial coordinate system, a radial trace of vents will show up as a vertical line, and a circumferential ring of vents will show up as a horizontal line (Fig. 12). Thus, the predicted distribution of vents should form a waffle-like pattern in the radial coordinate system, and indeed this is roughly what is shown in Figure 12. Along each axis is a histogram showing the relative abundance of vents for each dimension.

The histograms along the angular coordinate show peak concentrations for chains of vents radiating out in the SSW direction, the ESE direction, and a much more subdued signal approximately due N (Fig. 12). The strong concentration of ESE-oriented vents may be due to the presence of preexisting NW-SE–oriented structural lineaments such as the Archibarca lineament that runs from Llullaillaco to Cerro Galan (Richards and Villeneuve, 2001; Ruch and Walter, 2010), although the respective trends of these lineations appear oblique to one another by ∼15°–30°. Furthermore, it is thought that the lineation itself only migrates through the northern segment of the study area (e.g., Richards and Villeneuve, 2002, their fig. 2), where the vent concentration is remarkably low (Fig. 11). Similarly, the SSW-oriented vents may exploit the structural grain from the SSW/NNE–trending Pedernales-Arizaro thrust fault (Naranjo et al., 2014), which runs through the eastern edge of the uplift field. However, this mid-Miocene fault is buried by the younger volcanics at Lazure (Naranjo et al., 2014); so there is significant uncertainty as to its exact location in the study area. Additionally, there are a few vent lineations in the northern quadrants that appear to nucleate at the NW and NE corners of the apical graben and trend oblique to the radial direction (Fig. 11), and these chains form the diagonal clusters seen in the lower section of Figure 12. These off-axis sub-radial fractures are also seen in the analog experiments of Komuro et al. (1984) and Walter and Troll (2001) (Fig. 11B).

Along the radial dimension, there is a strong cluster of peaks at the approximate radius of the apical graben as initially observed by Pritchard and Simons (2002). There is also a secondary cluster toward the margins of the active uplift field, but this distribution appears more diffuse because the zone of deformation at Lazufre is elliptical and thus not radially symmetric. At distances out toward the margins of the long-wavelength dome (30–40 km), there is a general decrease in the circumferential distribution of vents, except where the outer edge of the dome intersects the northern, SSW, and ESE lineaments. This makes sense in the context of doming at Lazufre, because the peripheral ring fractures at the edge of the dome should be in a compressional regime and thus less likely to provide accommodation space for magma (e.g., Gudmundsson et al., 1997; Walter and Troll, 2001).

4. DISCUSSION

The long-term constraints on surface deformation provided by the local geomorphology and volcanology of both Uturuncu and Lazufre may help to discriminate between geophysical models of the plumbing systems beneath each region. Below we address how our data fit with various interpretations of subsurface magma movement.

4.1. Diapirism and Deflation at Uturuncu

The non-uniqueness of geodetic inversions at Uturuncu makes interpreting the source properties of the surface deformation difficult (e.g., Yun et al., 2006; Henderson and Pritchard, 2013). However, with our long-term topographic constraints, we can provide additional context for various model interpretations. Fialko and Pearse (2012) use geodetic observations of a peripheral deflationary moat around Uturuncu to justify a model of a rising diapir from the zone of partial melt beneath the volcano. Though their model matches well with the 18-year time series of surface deformation, they do not consider the effect of a rising diapir over centuries to millennia.

In a numerical modeling study, Burov et al. (2003) simulate topographic effects from brittle (plastic), viscous, and elastic deformation by a buoyantly rising magma body over an ∼1 m.y. model run. Model scenarios of early diapiric rise from 10 to 15 km depths, as is thought to be the case at Uturuncu (e.g., Fialko and Pearse, 2012), show a slowly evolving 200–600-m-high topographic dome as the magma body ascends. This topographic effect diminishes over time in their purely viscous scenarios as the diapir impinges on a viscosity contrast in the upper crust and begins to spread laterally, but given that geophysical (del Potro et al., 2013; Comeau et al., 2015) and geodetic (Fialko and Pearse, 2012; Henderson and Pritchard, 2013; Hickey et al., 2013) inversions suggest a source depth of ∼15 km, lateral spreading of the diapir at shallow crustal levels is likely not occurring. Furthermore, given the presence of seismicity thought to result from pressurization by the underlying magma chamber (Jay et al., 2012), we can infer that Uturuncu is susceptible to brittle deformation. These modeling results then suggest that topographic doming should accompany a rising diapir at Uturuncu, which we do not see in background topography of (Fig. 2). However, it is worth pointing out that the late-stage topographic profiles seen in the purely viscoelastic cases of Burov et al. (2003; see their Fig. 9) show little actual topographic deformation and, in some cases, even a topographic depression. If this were the case, it might suggest that the observed brittle deformation beneath Uturuncu is extremely small relative to the viscoelastic deformation it is experiencing. More constraints on the upper crustal rheology beneath Uturuncu should help shed light on this issue.

An alternative model involving deflation of a deep source and inflation of a shallow reservoir at Uturuncu has been proposed by others (e.g., Henderson and Pritchard, 2013). This model can also resolve the observed deflationary moat, and although its viability relies on assuming a much lower viscosity in the lower crust than the mid-crust to account for mass balance discrepancies between source and sink (S. Henderson, personal commun., 2015), it suggests inflation is periodic and transient rather than continuous over long time intervals. Viscous relaxation time in the crust may be on order of 102 years at Uturuncu (S. Henderson, personal commun., 2015); therefore if inflation episodes are transient, as geodetic observations suggest, then long-term deformation could relax viscously and hence go unrecorded in the topography. Alternatively, Walter and Motagh (2014) use mapped lineations proximal to Uturuncu to infer a model of long-term subsidence from normal faulting triggered by deflation at depth; however, interpreting the origin of lineations at Uturuncu is difficult given the presence of linear geomorphic features such as rivers and wind streaks, the plan-view orientations for which aren’t necessarily dictated by surface fracturing (e.g., Bailey et al., 2007). However, both ideas are consistent with the lack of observed long-term surface uplift at Uturuncu.

Taking the observational and theoretical evidence together, it is most likely that Uturuncu experiences pulses of transient inflation and either viscous relaxation or active deflation over millennial timescales with no net surface uplift, possibly from a lower-crustal reservoir feeding magma into the APMB (Fig. 14). This long-term “heavy breathing” is seen in other active silicic magmatic systems such as Yellowstone (e.g., Pierce et al., 2002); however, at Yellowstone there is a complex relationship between magmatic intrusion, migration of hydrothermal fluids, and surface deformation (Wicks et al., 2006). There does exist evidence of sustained surface uplift over thousand-year timescales at both Campi Flegrei (e.g., Cinque et al., 1985, where a small section has experienced a net inflation of 40 m since 4–8 ka, and Laguna del Maule, Chile, the shorelines of which have been differentially uplifted ∼67 m since the Holocene (Singer et al., 2015). Alternatively, if we are observing a rising diapir from the upper surface of the APMB, then we must be witnessing the nascent mobilization of the magma body because there hasn’t been sufficient growth to be observed in either the long-wavelength topography or shoreline deflections.

4.2. Pre-Eruption Doming at Lazufre

The long-lived topographic dome, lack of observed short-term peripheral subsidence, and shallower source depth make Lazufre’s deformation history unique from Uturuncu and more likely associated with a pre-caldera magmatic system (e.g., Ruch et al., 2008). As such, there are a number of questions that arise regarding the nature of long-term deformation here.

4.2.1. Why Is the Wavelength of the Long-Lived Topographic Dome Different from the Observed Short-Term Deformation?

Figure 2 shows that the radius of the topographic dome is wider by ∼10–15 km than the deformation footprint measured by InSAR (Figs. 2 and 14). The surface deformation pattern at Lazufre can be fit equally well by intrusion from either a thin sill or a flat-topped magma chamber (Ruch et al., 2008; Pearse and Lundgren, 2013; Remy et al., 2014). This non-unique result is predicted from modeling (Yun et al., 2006); therefore, evaluating the magma chamber source geometry requires additional constraints. The fact that the topographic dome is >400 ka at Lazufre suggests the presence of a much larger magma chamber than a single sill beneath the inflation source. Hence, the topography and apparent longevity of the dome at Lazufre may reflect chamber overpressure from accumulation of magma over many such transient sill intrusion events (e.g., Remy et al., 2014).

An additional mechanism for the difference in observed wavelength is the possible time-dependent deformation of a magma chamber with a viscoelastic shell. Given the topographic evidence for a long-lived central magma chamber beneath Lazufre, it’s likely that the rocks adjacent to the magma chamber have been thermally weakened beyond the brittle-ductile transition (e.g., Remy et al., 2014). In this scenario, over short timescales the pattern of surface deformation would reflect the elastic component of deformation from the inner magma chamber. Over time as the stress is transmitted outward toward the edge of the shell, the magnitude of radial displacement increases but at an exponentially decaying rate (Dragoni and Magnanensi, 1989; Jellinek and DePaolo, 2003; Karlstrom et al., 2010). Dragoni and Magnanensi (1989) estimate the characteristic response time τ as a function of the shear modulus μ, Poisson’s ratio ν, the viscosity η, and the ratio of the shell radius R2 to the chamber radius R1. Assuming a Poisson’s ratio of ν = 0.25, the relationship simplifies to τ = 1.8(η/μ)(R2/R1)3. Few constraints exist for estimating wall-rock viscosity, but Newman et al. (2001) calculate a range of values from 1015 to 1018 Pa s in the wall rocks of Long Valley caldera (e.g., Jellinek and DePaolo, 2003). Over this range of viscosities, a shear modulus of 30 GPa for the upper crust at 10 km depth (from the shear-velocity model of Ward et al., 2014), and a ratio of shell thickness to chamber radius of 0.625, the relaxation timescale may range from as little as a week to ∼5 years. If the magmatism is very young, a higher crustal viscosity of 1020 Pa s would yield a substantially longer response time (∼500 years). Given the topographic and volcanological evidence for long-lived magmatism (Froger et al., 2007; Ruch and Walter, 2010; this study), however, it is likely that the viscosities are toward the lower end of our estimates. Therefore, the lack of outward growth at Lazufre in the geodetic observations suggests that the presence of a viscoelastic shell may not be a primary mechanism to account for the differences in wavelength between the active zone of deformation and the long-wavelength topographic dome.

The simplest explanation for this discrepancy may come from the work of Komuro et al. (1984), who show through finite element modeling of an uplifting magmatic dome in an elasto-plastic medium that the diameter of the dome increases as a function of the chamber depth. Thus, the larger diameter dome may be a result of plastic deformation from a relatively deep sourced magma chamber beneath Lazufre, and this plastic (permanent) deformation is consistent with the inferred presence of brittle fractures from the volcanic vent lineaments.

Additionally, from their analysis, Komuro et al. (1984) find that chamber depth can be estimated for a given ratio of chamber width to dome width. Assuming that chamber width is approximately at the boundaries of the current uplift field, a chamber width:dome width ratio of ∼2/3 at Lazufre (e.g., Fig. 2) would yield an approximate chamber depth of ∼18 km (Komuro et al., 1984, their Fig. 9). Interestingly, this depth is similar to the tomographically imaged depth of the Altiplano-Puna magma body (e.g., Ward et al., 2014), but it is lower than the maximum predicted depth of 14 km by Remy et al. (2014).

4.2.2. How Does Lazufre’s Spatiotemporal Deformation Fit in with Our Understanding of Caldera Magmatism?

Recent thermomechanical models of caldera-scale magma chamber pressurization show uplift signals of ∼1 km as systems progress toward eruption (Gregg et al., 2012). Existence of the dome during the formation of peripheral volcanoes along Lazufre since at least 400 ka provides further evidence that they could be related to zones of high stress at the margins of a large silicic chamber (e.g., Froger et al., 2007), and thus it seems increasingly likely that the peripheral volcanoes at Lazufre are genetically related to a single pre-caldera magmatic system (e.g., Lipman, 1984; Froger et al., 2007; Karlstrom et al., 2015). This is further borne out by our analysis of vent locations that show superimposed radial and circumferential fracture patterns spatially correlated with the observed uplift and topographic anomalies, consistent with long-term inflation (e.g., Walter and Troll, 2001). As Froger et al. (2007) state, examining the bulk geochemistry of the peripheral volcanism may confirm this hypothesis that is strongly backed by the morphologic evidence.

If Lazufre’s dome began inflating at ca. 500 ka, then it would have a long-term average uplift rate of ∼1 mm/yr, which likely represents periods of quiescence separated by episodes of rapid inflation like we see today. If sill intrusions tend to cause surface uplift at the geodetically observed rate of ∼2.5 cm/yr, that implies that Lazufre has been quiescent over ∼96% of its lifetime. Presently only ∼1% of the ∼1100 volcanoes that comprise the Central Andes are actively deforming (see list at http://www.geo.cornell.edu/eas/PeoplePlaces/Faculty/matt/volcano_table.html). Depending on the longevity of active uplift pulses, which an extended time series of InSAR measurements will help uncover in the future, this could provide an explanation for why we only see active uplift at a small percentage of volcanic systems in the Central Andes.

If Lazufre is indeed a pre-caldera system, then an interesting question arises as to whether we can use measurements of the system geometry to learn anything about its potential eruption style. A critical metric for determining the mode of failure in models of caldera eruption is the caldera roof aspect ratio R, which is the ratio of magma chamber width to depth below the surface (Karlstrom et al., 2010; Gregg et al., 2012). As Remy et al. (2014) show, source depths from inversion of the geodetically observed deformation are poorly constrained and can range from 2 to 14 km depth. However, recent seismic imaging by Spica et al. (2015) finds no large chamber at Lazufre above 7 km, which places a hard constraint on the source depth inversions. With a modeled source geometry of ∼35 km width and depths ranging from 8 to 13 km for the inflating sill or chamber roof (Ruch et al., 2008; Remy et al., 2014), the magmatic system at Lazufre has a roof aspect ratio of R ∼0.22–0.36, which puts it near the threshold of critical failure via roof collapse for modest overpressures, if the chamber is sufficiently large (≥103 km3) (Gregg et al., 2012). Using a truncated cone model as a deformation source for their InSAR inversion, Remy et al. (2014) estimate a chamber volume of ∼2800 km3, well above the 1000 km3 threshold for the roof-collapse regime of Greg et al. (2012). However, an estimate of the entire magma chamber volume beneath Lazufre remains elusive. If the modeled deformation source corresponds to the top of a thicker chamber rather than a thin sill, a chamber height of 1–2 km would yield a volume of 700–1400 km3 for a 30-km-diameter magma body, indicating that overpressures not significantly greater than the 10 MPa presently observed at Lazufre (Ruch et al., 2008) could potentially lead to eventual failure of the overlying roof. Geophysical imaging of the magmatic plumbing system at Lazufre (e.g., Spica et al., 2015) may help determine what the true geometry of the magma chamber is and, hence, what mechanical conditions are necessary for a future eruption at the caldera scale.

5. CONCLUSIONS

Our observations using geomorphic markers to infer the paleogeodetic history of Uturuncu and Lazufre reveal similar characteristics yet distinct deformation histories for each region. Uturuncu shows no evidence for net inflation since ca. 16 ka from our lake shoreline and river profile measurements and has no signal of long-wavelength uplift that is predicted for a buoyantly rising diapir. This is likely because either: (1) the inflation is due to transient upward migration of melt from a deep source to a shallow sink (the APMB), and either dissipates or deflates over time; or (2) we are witnessing the first century of its development. At Lazufre, we find that the observed ∼500-m-amplitude, long-wavelength dome has likely existed as a local topographic high since at least 400 ka. Though the rate of change of surface elevation may be too slow to measure using thousand-year timescale measurements such as paleoshorelines, the longevity of deformation suggests the presence of a large magma chamber below the active sill-like intrusion at ∼10 km depth. Knowledge of surface displacements over the long timescales associated with magma intrusion can therefore be a useful tool in understanding the processes that shape large silicic magmatic systems.

We would like to thank Shan de Silva and Kerri Johnson for field guidance, assistance, and many conversations about volcanoes and geomorphology. We would also like to thank our drivers and Linda Hauge Kossatikoff from Lipiko Tours for their logistical support and Jose Antonio Naranjo for providing geologic map data. We thank Rodrigo del Potro for providing campaign GPS data. This manuscript greatly benefited from suggestions and criticisms by Guest Editor Matthew Pritchard and reviewers Jean-Luc Froger and Leif Karlstrom. The work for this project was supported by National Science Foundation grant EAR-0908850 to NJF.

1Supplemental File. Mapped flows within the study area in .kmz from high-resolution satellite imagery in Google Earth and ArcGIS. Please visit http://dx.doi.org/10.1130/GES01278.S1 or the full-text article on www.gsapubs.org to view the Supplemental File.