Abstract

Interferometric synthetic aperture radar (InSAR) and GPS measurements beyond 2010 are presented for the first time for the Lazufre volcanic center in the Central Andes. Vertical uplift at Lazufre was known to affect an area >50 km in diameter at rates exceeding 3 cm/yr between 1997 and 2010. Analysis of new InSAR data through August 2016 indicates that the spatial pattern of uplift is relatively unchanged but the amplitude of uplift has significantly decreased to <1.5 cm/yr since at least December 2011. We present a time-series inversion for InSAR data between 1996 and 2016 that is well fit by a double exponential model, with an inflection point occurring in 2006. For two continuous GPS stations installed within the deformation footprint in November 2010, we have determined vertical velocities through 2014 or 2015 (depending on the station) that agree with contemporaneous InSAR-derived velocities. Velocities from campaign GPS benchmarks established in November 2011 and reoccupied in March 2014 are also presented. We use a previously proposed model of an inflating sill at 10 km depth to explain geodetically observed displacements. Opening rates are halved (6.8 ± 1.25 × 106 m3/yr) compared to inferred values using data prior to 2010. Subsurface heterogeneity is accounted for by assigning elastic parameters based on local seismic tomography in a finite-element model. Surface displacements (or inferred volume change estimates) for heterogeneous models compared to homogeneous models are amplified by up to 7% within a 10 km radius of the center of uplift.

INTRODUCTION

The Lazufre volcanic uplift signal is centered between the summits of Lastarria volcano (25.168°S, 68.507°W, 5706 m elevation) and Cordon del Azufre volcano (25.336°S, 68.521°W, 5481 m elevation) in the Central Andes Volcanic Zone (CVZ) of northern Chile and Argentina (Fig. 1). These volcanoes were identified as having had Holocene activity based on remotely sensed morphological indicators such as post-glacial eruptive features (Francis and de Silva, 1989). Subsequently, extensive regional geological mapping has been conducted, which indicated ∼120 km3 of erupted material from vents in the vicinity of uplift since the late Pliocene (Naranjo, 2010). Although the majority of dated eruptive products have Pleistocene ages, the most recent lava flow at Lastarria has been dated at ∼2500 yr old (Naranjo, 2010), and at Cordon del Azufre flows have been dated to 0.3 ± 0.3 Ma (Wilder, 2015). The strike axis of current uplift is NNE-SSW and aligned with structural lineaments, and the spatial footprint contains a high concentration of volcanic vents (Ruch and Walter, 2010). Furthermore, topographical analysis combined with dated flows suggests persistent tumescence at the location of uplift since at least 400 ka (Perkins et al., 2016).

A striking characteristic of the Lazufre system is extremely vigorous degassing from the Lastarria edifice. In fact, localized uplift of the edifice of Lastarria volcano has prompted the hypothesis that Lastarria may act as a “pressure valve” for a deeper magmatic plumbing system (Froger et al., 2007; Ruch et al., 2009). Furthermore, recent in situ studies of gas composition at Lastarria have suggested a possible transition from a hydrothermal character between A.D. 2006 and 2009 (Aguilera et al., 2012) to a shallow magmatic source in 2012 (Tamburello et al., 2014). Remote sensing measurements of SO2 emissions at Lastarria from 2004 onwards also suggest a possible time dependence with a peak in 2006 (Carn et al., 2013). The recent temporal changes in gas emissions at Lastarria provide evidence for changes to the subsurface magmatic system, which motivates examining recent geodetic measurements for transient signals.

Uplift began at Lazufre between late 1997 and 2000, accelerating to a maximum line-of-sight (LOS) velocity of 3.5 cm/yr and affecting a broad region >50 km in diameter (e.g., Henderson and Pritchard, 2013). A second, smaller region (∼2 km across) of uplift of about ∼0.5 cm/yr starting in 2002 has been observed on top of Lastarria volcano (e.g., Froger et al., 2007; Ruch et al., 2009) but is not discussed in this paper because it is not well resolved by the data sets used in this study. Owing to the remote location of Lazufre, deformation measurements have mostly been made with interferometric synthetic aperture radar (InSAR); however, due to the end of Envisat (European Space Agency [ESA]) satellite observations in 2010, there was a gap in observations. For the first time, we present post-2010 InSAR observations along with continuous GPS observations at Lazufre to address the question of whether previous spatiotemporal surface uplift trends have continued in the past 6 yr.

To date, many source models of deformation at Lazufre have been constrained by a single LOS viewing geometry, utilizing data from European Remote Sensing [ERS] Satellites 1 and 2 from ESA and Envisat tracks 282 spanning July 1995 to May 2010 (Pritchard and Simons, 2002, 2004; Froger et al., 2007; Ruch et al., 2008). Owing to the nonunique nature of geodetic models, there has been some disagreement as to whether the source was a laterally propagating sill (e.g., Ruch et al., 2008; Anderssohn et al., 2009) or had a fixed geometry with variable opening rate. These models were based on a rectangular elastic dislocation model (Okada, 1985). All available ascending and descending data through 2010 were used recently by Pearse and Lundgren (2013) to address the question of lateral growth of the sill, and the best-fitting source was concluded to have a fixed geometry with a depth of 8 km, strike of 10° east of north, dip of 10°, length of 30 km, width of 20 km, and maximum opening rate of 5 cm/yr (implying volume change on the order of 0.01 km3/yr). Note that the depth in this case describes the middle of the sill below the local average elevation of 4.8 km. Remy et al. (2014) modeled, with different methodologies, the same InSAR data set, adding campaign GPS observations spanning 2006–2008, and pointed out that either a truncated cone-shaped reservoir or a sill with a geometry similar to that found by Pearse and Lundgren (2013) can explain the observations. Models generally agree that the source must be ellipsoidal or rectangular at ∼10 km depth with a northerly strike (0°–25°) and slight dip to the east (0°–10°E). Both recent studies agreed that there was no evidence for a propagating sill and that a fixed geometry was consistent with the data (e.g., Pearse and Lundgren, 2013; Remy et al., 2014).

All geodetic models to date for Lazufre have assumed that the crustal structure is homogeneous and elastic. In recent years however, seismic tomography and conductivity models have been published, which help constrain subsurface heterogeneity under the Lazufre region. In particular, magnetotelluric data revealed a strong conductor under the deformation anomaly dipping eastward to the base of the crust. Modeling suggests that the feature could be due to 5–8 vol% partial melt which is feeding current sill inflation in the upper crust (Budach et al., 2013). A follow-up study with additional local instrumentation defined several distinct low-conductivity anomalies—one directly under Lastarria volcano (1–10 Ω·m) to 1 km depth, another south of the edifice (5–10 Ω·m) at 7–8 km depth, and finally an anomaly (0.1–1 Ω·m) at 5–15 km depth under the deformation centroid (Diaz et al., 2015). These results match quite well with local seismic tomography in which roughly coincident low S-wave velocities were observed: 1.2–1.8 km/s, 1.5–2 km/s, and 2.3 km/s, respectively (Spica et al., 2015). The prevailing interpretation is that magmatic fluids migrate upward to ∼6 km from a reservoir at 10 km depth, thereby providing heat to the extensive and shallow hydrothermal system (Spica et al., 2015).

In the following sections we present up-to-date InSAR data analysis, and measurements from recent continuous and campaign GPS surveys, extending geodetic measurements through August 2016. We validate measurements by comparing against past and contemporaneous data sets and constrain the time dependence of uplift at Lazufre using InSAR time-series analysis. InSAR measurements after 2010 are then jointly inverted for sill opening rates in order to assess implications for source dynamics and volume changes. Finally, we examine the effect of subsurface elastic heterogeneity on surface uplift and source models by incorporating newly available seismic tomography into a finite-element model.

DATA and METHODS

InSAR Measurements

To densify InSAR observations prior to 2010, we processed Advanced Land Observing Satellite (ALOS; operated by the Japan Space Agency [JAXA]) and Envisat scanning synthetic aperture radar (ScanSAR) data sets with acquisitions between October 2005 and February 2011. In addition, we processed recent InSAR data sets spanning May 2010 through August 2016 (Table 1). Strip-map data from the ALOS-1 and TerraSAR-X (operated by the German Space Agency [DLR]) (hereafter TSX) satellites were processed with ROI_PAC software (Rosen et al., 2004). Envisat ScanSAR data were processed with a modified ROI_PAC workflow (Liang et al., 2013). Data from COSMO (Constellation of Small Satellites for Mediterranean basin Observation, operated by the Italian Space Agency [ASI])-SkyMed and RADARSAT-2 (Canadian Space Agency [CSA]) (hereafter CSK and RS2) were processed with ISCE 2.0 software (Rosen et al., 2012). Data from Sentinel-1 were processed with ISCE release 201604, which included a workflow specific to the unique Terrain Observation with Progressive Sans (TOPS) observation mode of the Sentinel-1 sensor. We removed the topographic signal with the Shuttle Radar Topography Mission (SRTM; NASA) digital elevation model (Farr et al., 2007). For ALOS and Envisat we used topography with a 90 m posting, and for CSK, TSX, and RS2 processing we utilized the recently released 30 m native-resolution topographic data. All interferograms were masked (pixels with coherence < 0.3), spatially filtered (Goldstein and Werner, 1998), and finally unwrapped with the SNAPHU program (Chen and Zebker, 2001), before being converted to surface displacements. Upon geocoding, images have a 90 m posting; however, spatial filtering in the processing chain creates an effective ground resolution of ∼200 m.

C Band Data: Envisat ScanSAR, RADARSAT-2, and Sentinel-1

Envisat ScanSAR data have been largely overlooked, as only Anderssohn et al. (2009) were able to calculate a single coherent interferogram. We processed the ScanSAR subswath containing Lazufre for scenes with a burst overlap of at least 10%. This included 29 interferograms from 19 dates in Wide Swath track 361 (October 2005–November 2009) and eight interferograms from nine dates in WS89 (December 2005–October 2009). We stacked these interferograms to generate maps of average ground velocity (see Henderson and Pritchard [2013] and Supplemental Material1). The spatial footprint and signal amplitude in stacks agree with those from Envisat strip-map data (Fig. 2). We did not reprocess the Envisat strip-map interferograms, but used the results from Remy et al. (2014) and Henderson and Pritchard (2013).

We calculated an RS2 stack with 15 descending interferograms from six acquisitions covering 2.25 yr between April 2014 and August 2016 (Fig. 2G). Due to the persistence of topography-correlated phase delays in half of the interferograms, we used an empirical correction in which the linear trend of unwrapped phase versus elevation in non-deforming far-field areas was subtracted before stacking (e.g., Froger et al., 2007). Deforming pixels were masked by defining a region within a 50 km radius of the center of maximum uplift. This normalization resulted in a stacked result with a mean velocity centered on zero and standard deviation of 1.5 mm/yr (Fig. S6 [footnote 1]).

Sentinel-1 is a constellation of two C-band SAR satellites launched by ESA in April 2014 and April 2016. Although we don’t discuss the Sentinel-1 data set in this paper, we include a single interferogram spanning December 2014–August 2016 (Fig. S7 [footnote 1]) to demonstrate the promise of measurements from this system going forward.

L Band Data: ALOS and UAVSAR

There are observations of Lazufre from U.S. National Space and Aeronautics Administration’s (NASA) airborne Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR) sensor spanning 2013–2015 (24 March 2013, 29 April 2014, 01 April 2015). Interferograms are processed by the Jet Propulsion Laboratory (JPL; Pasadena, California) and made available via the website uavsar.jpl.nasa.gov. For a single L-band interferogram, 16.8 cm of cumulative deformation would be represented by a single cycle of phase, but at most we expect 3 cm of relative deformation between 2013 and 2015 at Lazufre. Consequently, the uplift signal does not stand out from background noise in currently available interferograms, but we mention the UAVSAR data set here because future acquisitions could prove useful.

There are 16 ALOS L-band (23.6 cm wavelength) acquisitions between 09 February 2007 and 20 February 2011, providing additional observations toward the end of the Envisat mission. Unfortunately, many of these acquisitions show strong ionospheric signals, which have been noted in other regional studies (e.g., Fournier et al., 2010). We stacked 19 successful interferograms and found that the spatial and temporal patterns of deformation agree with those of previous data sets (Fig. 2).

X Band Data: TerraSAR-X and COSMO-SkyMed

TSX provides X-band (3.1 cm wavelength) observations of Lazufre from three different orbits (o134, o20, o111). Unfortunately, the smaller footprint of TSX acquisitions results in only partial coverage of the complete uplift footprint. Observations with TSX began in April 2008 (o134), but most acquisitions are after 2011 with the most recent in August 2016. We created stacks from the time span 2014–2016 to directly compare with C-band measurements (Fig. 2).

CSK is also an X-band platform, and we generated data stacks using two ascending and descending satellite passes. Due to a change in incidence angle and polarization in 2012, we are unable to create stacks from a single set of interferograms. Instead, we divided acquisitions into two sets before and after May 2012, creating four stacks in total (Supplemental Fig. S4 [footnote 1]).

InSAR Time Series

Stacks represent the average velocity over the elapsed time of observations, but more detailed temporal evolution can be determined by a least-squares time-series inversion using the set of overlapping interferograms (e.g., Berardino et al., 2002; Schmidt, 2003). In brief, we assumed constant velocity between consecutive dates and solved for cumulative displacement relative to the first date for each pixel (see Supplemental Material [see footnote 1] for further detail). We calculated an InSAR time series for RS2 and TSX o20 data sets only, because they have very similar viewing geometry compared to the long-term 1995–2010 ERS–ENVISAT track 282 time series (Table 1). Each data set was inverted independently, and there is a gap with no connected observations between May 2010 and December 2011. To create a continuous time series between 1996 and 2016, we fit linear rates before and after the gap in observations (specifically before and after February 2011). These rates are interpolated across the gap in observations in order to obtain the final 1995–2016 times series shown in Figure 3.

GPS Observations

Three continuous GPS stations were installed in the deforming region in November 2010 and November 2011. These stations could only be installed in the country of Chile, and consequently measurements are limited to the western side of the deformation footprint (Fig. 1). Station LCEN was placed as close as possible to the center of deformation given the existing network of roads, but is ∼12 km southwest of the centroid of vertical surface displacement (25.259°S, 68.483°W) based on the analysis of Remy et al. (2014). Station LLST was installed on the edifice of Lastarria volcano and is one of the highest-elevation continuously operating GPS stations in the world (5272 m), and station SOCM was intended as a far-field reference station at approximately equal longitude compared to the other stations (∼300 km east of the Nazca trench). We supplemented continuous measurements with campaign surveys, installing a total of 10 benchmarks in 2011 (Fig. 4 and Supplemental Table S1 [footnote 1]). Two campaigns were conducted in November 2011 and March 2014. All of the benchmarks were surveyed in 2011 and four benchmarks were resurveyed in 2014, with each survey consisting of a minimum of 48 consecutive hours of data collection. Table 2 lists details of the occupation histories of the sites used in this study.

We used the GIPSY/OASIS II software developed by the JPL to produce precise point positioning (PPP) solutions for the GPS data used in this study (Zumberge et al., 1997). We used reprocessed satellite orbits and clock products from JPL and transformed the daily solutions into the International Terrestrial Reference Frame 2008 (ITRF2008). These daily solutions were combined in a linear least-squares inversion to estimate velocities at each GPS site. The velocities were then referenced to the GEODVEL (a model of major tectonic plate velocities from space observations from GPS, SLR, VLBI, and DORIS over 25 years) estimate of South America plate motion (Argus et al., 2010). Further details on the processing scheme can be found in Fu and Freymueller (2012).

Our two campaign surveys occurred in March and November, which introduces the possibility of a bias due to seasonal effects. Water loading in the Amazon Basin generates a strong vertical seasonal undulation in the GPS data throughout northern South America (Davis, 2004; Bevis, 2005) as well as radial horizontal motion from the center of the loading (Fu et al., 2013). Estimates of seasonal deformation from the GRACE (Gravity and Recovery Climate Experiment by NASA) satellite capture the broad-scale regional signal but may not be able to detect more localized seasonal effects, and our sparse continuous GPS (cGPS) network makes it difficult to estimate seasonal deformation from a collection of regional stations. In fact, stations LCEN and LLST are nearly 300 km from the nearest cGPS station to the north (station CJNT) on the border with Bolivia, and 300 km from the nearest cGPS station to the south (station COPO) in Copiapo, Chile.

Instead of attempting to remove an estimated seasonal signal from the data, we focused on gauging how much uncertainty campaigns in different seasons added to our velocity estimates. To accomplish this, we ran three velocity solutions in which station LCEN only had a small number of days of data, turning it into the equivalent of a campaign site. In the first solution, LCEN had days of data corresponding to the campaign surveys in November 2011 and March 2014. In the second and third, LCEN had days of data from March 2011 and 2014 and November 2011 and 2014, respectively. The three estimates did display differences, with LCEN showing a maximum difference of 2 mm/yr between the velocity estimates in the vertical component, 1.3 mm/yr in the east horizontal component, and 1.7 mm/yr in the north horizontal component. We added these values to our formal uncertainty estimates for the campaign site velocities.

Table 3 lists the horizontal and vertical GPS velocities in the Lazufre region. The dominant linear horizontal motion is due to interseismic strain accumulation. Vertical velocities suggest a concentration of uplift around the Lazufre region compared to the far-field site (station SOCM), with the highest uplift rates apparently located southeast of Lastarria volcano. However, the highest velocities for campaign stations with sparse time series are large. For example, at station RDHL (Fig. 1) we estimate 20.8 ± 9.5 mm/yr, compared to continuous station LLST where we estimate 7.9 ± 0.8 mm/yr.

RESULTS AND MODELING

InSAR Stacks after 2011

Clearly, new InSAR data show that uplift at Lazufre has continued at least until August 2016. However, comparisons between InSAR uplift rates are complicated by different viewing geometries associated with different satellites. In order to assess a possible change in the uplift rate we used the east-west (Ux) and vertical (Uz) displacements calculated by Remy et al. (2014) from a subset of Envisat strip-map data spanning May 2006–November 2008 (Fig. 1) and then projected them into the CSK, TSX, and RS2 LOSs. We benchmarked the Ux and Uz displacement maps by comparing with our independent ALOS and Envisat ScanSAR stacks spanning 2003–2010 (Supplemental Fig. S1 [footnote 1]), which were not used to estimate the Ux and Uz displacements by Remy et al. (2014). Agreement is remarkable between the Envisat data sets even though the time period of the data sets and their LOSs are different (∼20° and 40° incidence for strip map compared to 29° and 33° incidence for ScanSAR), suggesting a relatively constant rate of uplift between 2003 and 2010. However, the mean rate from ALOS data is ∼3 mm/yr less than that observed with Envisat. Therefore, it could be argued that the rate decreased slightly between May 2010 (the end of the Envisat time series) and February 2011 (the end of the ALOS time series), but too few data points are available during that short time period to assess this robustly.

Projecting the 2006–2008 Ux and Uz components derived from Envisat into the CSK, RS2, and TSX LOSs results in rates that are 1.5 cm/yr higher than the observations made between 2014 and 2016. To assess the significance of the overestimate with respect to data noise, we used the approach of Finnegan et al. (2008) to estimate a minimum bound of InSAR stack uncertainties. This method compares InSAR stacks with different incidence angles covering the same geographic region. We used the RS2 and TSX o111 stacks and converted the mean LOS velocities into pseudo-vertical rates by dividing them by the cosine of the incidence angle. The co-registered stacks are then compared in a correlation plot. Gaussian noise of variable amplitude is added to a 1:1 correlation until the spread matches the data. We find that noise with a magnitude of 1.2 mm/yr matches the scatter between the rate maps, thus we use this figure as a proxy for the uncertainties on the LOS velocity in our interferogram stacks (Supplemental Fig. S5 [footnote 1]). This figure is about one order of magnitude smaller than the residuals of the 2014–2016 stacks minus the synthetic interferograms calculated with the 2006–2008 data (i.e., 1.5 cm/yr), hence we conclude that although the spatial footprint of the inflation signal has remained relatively unchanged, the LOS uplift rate has decreased from >3 cm/yr to 1.5 cm/yr.

We also decomposed our recent InSAR data into Ux and Uz components to constrain the ratio of maximum horizontal to maximum vertical displacement which can be diagnostic of source geometry. RS2 descending (18 April 2014–29 August 2016), TSX descending o111 (20 July 2014–06 August 2016), and CSK ascending tracks (01 August 2013–06 August 2016) were used. Maximum vertical uplift is determined to be ∼1.5 cm/yr with a Uxmax/Uzmax ratio of 0.20, somewhat lower than 0.27 obtained using data from 2006–2008 (Remy et al., 2014) and in the range that is commonly cited as diagnostic of a flat-topped source (e.g., Dieterich and Decker, 1975). We caution, however, that using this ratio as diagnostic of source geometry is only valid for a homogeneous crust (e.g., Geyer and Gottsmann, 2010) – in the Finite-Element Model section we explore the impact of a heterogeneous crust, determined by seismic studies, upon the ground deformation.

InSAR and GPS Comparison

The continuous GPS station LCEN was installed in November 2010, and data indicate a constant vertical velocity of 8.3 ± 0.3 mm/yr through April 2015 (Fig. 5). We have assumed that the vertical rate is purely due to volcanic inflation and does not contain interseismic signal. We believe this is reasonable given the prediction of <2 mm/yr uplift of interseismic coupling models (Béjar-Pizarro et al., 2013) and the measurement of negligible uplift at the reference station SOCM far from Lazufre. Importantly, there is no apparent residual trend after removing this linear vertical velocity from the LCEN time series aside from seasonal oscillations. LCEN is 12 km southwest of the inferred center of deformation, where the components of deformation from May 2006–November 2008 predict uplift of 12 mm/yr and westward motion of 9 mm/yr (Remy et al., 2014). The predicted magnitude of westward motion should clearly reduce interseismic motion at LCEN compared to regional stations (because eastward regional rates 300 km from the trench are 15–20 mm/yr in the International GNSS Service 2008 reference frame). However, this is not observed—in fact, the eastward motion of stations LLST and LCEN, 16.8 ± 0.6 and 19.2 ± 0.4 mm/yr respectively, are comparable to that of far-field reference station SOCM, 19.0 ± 0.5 mm/yr, far from the deformation field (Table 3). The predicted vertical rate at LCEN of 12 ± 1 mm/yr using Envisat data from 2006 to 2008 is much higher than the observed rate of 8.3 ± 0.3 mm/yr between November 2010 and April 2015.

Both InSAR and GPS measurements since 2010 show evidence for a decrease in uplift compared to pre-2010 rates. The question is, when did the deceleration occur? Given the variable spacing of consecutive SAR acquisitions, it is difficult to determine a precise inflection point in the uplift rate at Lazufre. Indeed, it can be seen in Figure 3 that for any 5 yr period it is difficult to distinguish an accelerating rate from a constant rate of uplift with InSAR data alone. Only considering the full data set does a clear trend emerge: Applying a double exponential model to the time series results in an inflection point in the time series in 2006 (see the Double Exponential Pressure Evolution section). Unfortunately, this inflection predates the installation of GPS, but the lack of evidence for any significant residual trend in the GPS time series since November 2010 suggests that a large deceleration took place over 5 yr between 2006 and 2011.

Elastic Okada Sill Model

Motivated by the observed change in deformation rate, we tested whether the sill model of Pearse and Lundgren (2013)—which uses distributed sill opening constrained by data from before 2010—is consistent with the InSAR data sets after 2010. Volume change for the sill is found via a linear inversion for opening on 2 × 2 km patches within a predefined sill geometry (Table 4). If descending data only are used, this model fits most data sets well with a volumetric change of 0.006 km3/yr, but residual amplitudes of ∼4 mm/yr are obtained near the maximum uplift in the RS2 stack. Including ascending TSX data changes this estimate to a volume change of 0.0068 km3/yr (6.8 × 106 m3/yr), which is approximately half of the previously estimated rate between 2003 and 2010 (Remy et al., 2014). To estimate uncertainties, we ran a Monte Carlo test in which 100 realizations of noise were added to the model results and re-inverted. The amplitude of noise was estimated from the standard deviation of values in non-deforming regions in the RS2 stack, 1.5 mm/yr (Supplemental Fig. S6 [footnote 1]). We obtained a final volume change estimate of 6.8 ± 1.25 × 106 m3/yr (Fig. 6). Maximum sill opening values are 2.5 cm/yr, and contours of opening are shown in Figure 7.

Finite-Element Model

Numerical studies have shown that upper crustal heterogeneities can significantly influence parameter estimates and modeled surface displacements from homogeneous half-space crustal deformation models (e.g., Masterlark, 2007). In order to evaluate the role of crustal heterogeneity on surface uplift at Lazufre, we constructed a finite-element model using Pylith 2.0 software (Aagaard et al., 2013). We used a rectangular mesh with a horizontal (dip = 0°) sill-like reservoir located at 10 km depth, which is the approximate depth found in previous work (Table 4). An isometric view of the mesh with sill opening and surface displacements is shown in Figure 8. All nodes along the 10-km-depth slice can be assigned an initial displacement or traction boundary condition, so that we can test rectangular sill geometries of various lengths and widths. Although most studies suggest a strike of >20°, we keep a strike of zero for mesh simplicity and note that solutions can simply be rotated. The model domain of 200 × 200 × 100 km is sufficiently large to utilize zero-displacement boundary conditions along edges and the bottom. We use hexahedral elements with 2 km sides, resulting in 500,000 elements in total. A map view depicting the extent of the finite-element model in relation to seismic tomography used to assign material properties is shown in Figure 9.

To validate the implementation of our finite-element model, we compared solutions for a point sill opening with the solution for a rectangular dislocation in an elastic half-space (Okada, 1985, 1992). The analytical solution is expected to have some deviation from the finite element model solution due to finite element model nodal displacements tapering to zero. For example, for the simplest initial condition of opening at a single node, the dislocation takes on a pyramidal shape that has a base length of twice the element size. Opening more closely approximates the analytical rectangular dislocation geometry as the mesh is refined and more nodes within the sill area are prescribed displacements. Performing a convergence test, we find that a cell size of 2 km agrees with the analytical solution to within <5% (Supplemental Fig. S10 [footnote 1]). We emphasize that a “point sill” is defined by the sill horizontal dimensions being much less than depth (e.g., ratio of <1/5). In general, finite sources with distributed opening patches in a plane result in shallower source depths (e.g., Fialko et al., 2001a, 2001b).

To add realistic heterogeneous elastic properties to the finite-element model, local seismic tomography is required. Recently, a regional (12°S–42°S) ambient noise tomography (ANT) data set became available that covers the Lazufre region (Ward et al., 2013). This analysis utilized 330 broadband stations throughout South America operating at different times over the period May 1994 to August 2012 and provides a model of seismic shear wave velocities (Vs) to a depth of 50 km (at 10 km horizontal and 1 km vertical grid postings). We present depth-averaged profiles of P-wave velocity (Vp) and S-wave velocity (Vs) for comparison with homogeneous half-space values (Fig. 10). Importantly, the ANT technique is insensitive to sharp impedance contrasts, which are better detected with receiver function analysis. A joint inversion for a Vs velocity model that combined ANT and receiver functions at Uturuncu volcano (300 km north) led to a significant changes in the three-dimensional (3-D) subsurface velocity field (Ward et al., 2014). Similarly, we suspect that the inclusion of local receiver functions at Lazufre would reduce velocities in the mid-crust and reduce the vertical extent of low-velocity zones, but this work has not yet been done.

We also utilized a local Vs velocity tomographic model centered on Lastarria volcano using seismic data from eight stations from January to March 2012, in addition to a densified local network consisting of 18 stations deployed between February and March 2008 (Spica et al., 2015). The tomographic model produced from this analysis has a smaller spatial extent compared to that of Ward et al. (2013) but higher resolution (100-m-sided cells). We present a depth-averaged profile of this model to give a first-order comparison to others discussed so far (Fig. 11).

The final step in converting tomographic models into elastic material parameters requires assuming a relationship between Vp, Vs, and density. Often, a constant Vp/Vs ratio is used, because the assumption of Poissonian material requires that Vp/Vs = 1.732. Alternatively, empirical relations derived from field studies and laboratory experiments can be used to obtain both Vp and density, for example at Mount Etna (Italy; Currenti et al., 2007). We use the empirical equations of Brocher (2005) to obtain the elastic material moduli that represent the crust in our finite-element model (Table 5).

Depth-averaged velocity models, though a simplification of the true 3-D variation, provide a simple assessment of the effect of material property variation with depth. Indeed, the tomographic model of Ward et al. (2013) contains predominantly 1-D variation with depth. For the same 10 cm of opening on a 10 × 10 km sill at 10 km depth, this heterogeneous configuration causes a subtle 1.7% increase in maximum vertical displacement and 3.5% increase in maximum radial displacement (Fig. 12). We also utilized the 3-D tomographic data set from Spica et al. (2015) in order to quantify the effect of this model on computed surface displacements. We again utilized the finite element model domain with a 10 × 10 km sill at 10 km depth and an opening of 10 cm. We find that vertical displacements are amplified by up to 6.6% and radial displacements are increased by up to 4.9% (Fig. 12). Assuming that the seismic tomography does not have significant temporal variation during the period of geodetic observations (1998–2016), existing elastic models of sill opening likely overestimate opening values by <10%.

DISCUSSION

Double Exponential Pressure Evolution

Studies to date have considered uplift at Lazufre to be due to an elastic sill-opening process, and consequently an overpressure time function identical to the functional form of the uplift time series could be invoked to explain observations. The magma recharge model presented by Le Mével et al. (2016) provides a straightforward explanation of how a double exponential trend in ground deformation time series (Fig. 3) can arise. The theory is developed for spheroidal reservoirs; however, we assume a sill can be considered a limiting geometry of a very thin ellipsoid. Essentially, the theory predicts that for a reservoir fed by a conduit with a ramping inlet pressure, the reservoir’s pressure evolution ΔP0(t) is given by: 
graphic
where ΔPi is the maximum reservoir overpressure achieved at time t* and τ is an exponential time constant that depends on magma supply dynamics. We assume that ΔPi is directly proportional to maximum ground uplift, which is at least true for the point source approximation for a penny-shaped crack (Fialko et al., 2001a) in an elastic half-space. In other words, we assume that the displacement time series in Figure 3 is directly proportional to overpressure so that we can invert the time series for the temporal terms t* and τ. We also assume that inflation began in October 1997, although there is a large uncertainty on the uplift onset due to the lack of data between this date and March 2000.

We used the nonlinear least-squares Levenberg-Marquardt algorithm to fit Equation 1 to our displacement time series, resulting in ΔPi = 48.2 ± 1.2 (scaled pressure units, discussed subsequently), τ = 9.8 ± 0.8 yr, and t* = 8.7 ± 0.4 yr. The uncertainties for these parameters were estimated by inverting synthetic data predicted by the best-fit model with synthetic noise created from a diagonal covariance matrix (Lohman and Simons, 2005) under the assumption that the data sets are spatially and temporally uncorrelated. We repeated this process for 100 synthetic noisy data sets, such that the model parameter uncertainties are the standard deviation of the ensemble of inversion results. We consider the stated uncertainties to be underestimates given the assumption of spatially and temporally uncorrelated noise.

The best-fit model (Fig. 3) requires linearly increasing overpressure over 8.7 yr (from October 1997 until August 2006) due to magma input, after which overpressure is held constant. The characteristic time scale, τ = 9.8 yr, is related to the dynamics of magma ascent assuming the shallow chamber is connected to a deeper source via a direct conduit. The characteristic time scale of ascent is related to conduit length and radius, and magma properties such as viscosity and compressibility. However, most of these parameters are poorly constrained (e.g., the conduit could extend as far as the Moho, which from receiver function studies is estimated to be 60 km deep in this region; Wölbern et al., 2009). We note that the inflection in pressure history in late 2006 could be associated with an independent observation that SO2 emissions at Lastarria volcano appeared to have peaked in 2006 (Carn et al., 2013). However, this gas analysis only considered measurements between 2004 and 2008.

The absolute value of overpressure cannot be determined independently without knowing the original reservoir size and crustal elastic moduli. The following formula relating overpressure to surface displacements for a thin sill-like reservoir holds when sill radius is much less than depth, and therefore serves as only an approximation of the absolute overpressure evolution at Lazufre: 
graphic

Considering values of shear modulus μ = 30 GPa, sill radius a = 5 km, d = 10 km, and Poissonian crust (ν = 0.25), we find that overpressure scales as 0.025 MPa/mm. Therefore, 400 mm of observed cumulative uplift can be explained by a pressure function in which overpressure grew by 1 MPa/yr between 1997 and 2006, after which point overpressure has been holding at a relatively constant value of 10 MPa/yr. This calculation illustrates the dynamics of the hypothesized pressure evolution at Lazufre; however, the stated pressure values are crudely approximated and should be considered upper bounds because for a given volumetric change, overpressure is greatly reduced for larger chambers (e.g., Fialko et al., 2001a, 2001b).

Importantly, the transition to constant overpressure in late 2006 predicts that uplift at Lazufre will decelerate until it asymptotically reaches a total cumulative displacement of 50 cm by 2040. Furthermore, if excess pressure is relieved, we expect a transition to subsidence. A recent geomorphic study identified the existence of a long-wavelength topographic dome at Lazufre (500 m height across 70 km distance) which has existed for at least 400 ka (Perkins et al., 2016). The development this topography can be explained by an average relative uplift rate of 1 mm/yr, an order of magnitude less than geodetically determined rates over the last two decades. Therefore, the active uplift episode is likely one of many occurring in the same spatial location, but separated by long periods of quiescence and possibly even periods of subsidence as seen at nearby volcanic centers such as Cerro Overo and Cerro Blanco, Chile (e.g., Henderson and Pritchard, 2013).

Comparison to Other Volcanic Systems

Similar patterns of ground uplift rate increase followed by decreasing rates have been observed for analogous deformation episodes at Three Sisters (Oregon; 1995–2010), Yellowstone (Wyoming, USA; 2004–2014), and Laguna del Maule (Chile; 2007–2014) in which one time constant describes initial uplift and another prolonged deceleration (Le Mével et al., 2015). This “double exponential” trend requires an initially increasing pressurization (Le Mével et al., 2016), and eventual deceleration can be explained either with a constant inlet pressure into a compliant elastic reservoir (e.g., Le Mével et al., 2016) or time-dependent deformation in a viscoelastic region (e.g., Newman et al., 2006).

Changes in caldera-scale inflation signals have been shown to correlate with the occurrence of seismic swarms: for example, at Yellowstone (e.g., Waite and Smith, 2002; Chang et al., 2007; Tizzani et al., 2015) and Long Valley Caldera (California, USA; e.g., Newman et al., 2001; Hill et al., 2003; Montgomery-Brown et al., 2015). In this conceptual model, caldera inflation is produced by basalt intrusion from a deep source into the brittle upper crust. When a threshold pressure is reached, a sealed layer inferred to be located within the brittle and ductile transition is breached, thus allowing fluids exsolved from the magma to be injected into the brittle shallower parts of the volcanic system where the fluid pressures are much lower. As this processes happens, the ground uplift is either reversed to subsidence or stops and is accompanied by seismic swarms with fluid-like signatures including b-values that deviate from one and stress tensors that deviate from the background stresses (e.g., Waite and Smith, 2002; Hurwitz and Lowenstern, 2014).

The uplift episode at Long Valley caldera between 1995 and 1999 also exhibited exponential growth coinciding with a large seismic swarm, followed by exponential decay. One possible mechanical model behind this deformation pattern is the intrusion of magma into a reservoir surrounded by a viscoelastic shell (Newman et al., 2001, 2006). If the crust is hot enough, it is likely that the area that surrounds the inflation source behaves as a viscous medium on time scales of decades, creating a time-dependent uplift after the intrusion of a single batch of magma.

Whether deceleration in uplift is due to a ramped pressurization history or viscoelastic relaxation process cannot be determined with currently available data. Determining the correct model could be aided by repeated microgravity studies to constrain evolving mass distributions, as has been proposed for the deformation episode at Three Sisters volcano in Oregon (Dzurisin et al., 2009). Furthermore, continuous gas or seismic monitoring could help constrain the timing of potential magmatic additions that are below the detection threshold of deformation measurements. Finally, the viscoelasticity of rocks in volcanic regions is poorly constrained, but heat flow measurements and seismic attenuation tomography could narrow in on appropriate values for the Lazufre system. In summary, with geodesy data alone it is not possible to disentangle subsurface processes responsible for the observed surface uplift rate changes at Lazufre, and it could be that a combination of magma injection, volatile exsolution, and viscoelastic rheology are responsible for the observed transient signal.

Correlation of Deformation with Other Data Sets

Unfortunately, there has not been a continuously operating seismic network at Lazufre to compare seismic activity with deformation. Seismic stations have only been deployed at Lazufre during January–March 2008 (Spica et al., 2015) and November 2011–March 2013 by the PLUTONS project (e.g., Alvizuri and Tape, 2016; Shen et al., 2016; Farrell et al., 2017), and a single permanent seismometer from OVDAS (Chilean Volcano Observatory) was installed on Lastarria volcano in December 2013. A detailed study of the shallow seismicity at Lastarria has not yet been completed, but the limited available time series will make comparison with the deformation time series difficult.

Similarly, one might expect changes in the degassing Lastarria volcano system to be correlated with changes in deformation rate, especially given the hypothesis that Lastarria degassing is the “pressure valve” for the Lazufre magmatic system (e.g., Froger et al., 2007). Gas output at Lastarria volcano is by far the greatest in the Central Andes Volcanic Zone, at times exceeding 13,000 tons/d (e.g., Tamburello et al., 2014). Temporal changes in gas composition between 2006 and 2012, specifically higher proportions of magma soluble gases H2O, SO2, and HCl and lower proportions of CO2, could reflect ascent of magma (Aguilera et al., 2012; Tamburello et al., 2014). SO2 flux at Lastarria volcano was measured at 884 tons/d with differential optical absorption spectroscopy on 27 November 2012 (Tamburello et al., 2014). On 17 March 2012, 16 March 2013, and 22 November 2014, SO2 emissions were measured at ∼800 tons/d, but large uncertainties in these measurements do not permit the conclusion of a decreasing emission rate between 2012 and 2014 (Taryn Lopez, 2016, personal commun.). More recently, the SO2 flux spanning 21 and 22 November 2014 was measured at 861 ± 250 tons/day, such that the infrequent measurements and large uncertainties do not permit concluding any recent trends. Magma ascent as evidenced by compositional gas changes could be accompanied by increasing uplift rates, which have been observed between 1998 and 2006. As the volume of gas emissions increases, reservoir pressures may be relieved, causing the decrease in uplift observed since 2006. These ideas assume an instantaneous connection between reservoir processes and surface gas emissions and can be tested in the future with dense gas sampling and integrated modeling.

CONCLUSION

We have demonstrated that deformation at Lazufre continued past 2010, but at a reduced rate. InSAR data through August 2016 show spatial uplift patterns like those observed with analysis through 2010, but maximum LOS uplift rates of up to 1.5 cm/yr. Continuous GPS time series for a point 12 km from the maximum deformation shows a clear linear vertical displacement rate of 8.3 ± 0.3 mm/yr between November 2010 and April 2015, consistent with the reduced rate measured with InSAR. Assuming the same sill-opening model of previous studies, we determine that opening rates (and therefore minimum volume changes) are approximately halved to match observations since 2010. Variations in modeled surface displacements due to subsurface heterogeneity known from current seismic tomography are <10% and mostly confined to within 10 km of the uplift source.

The deformation time series from InSAR between 1997 and 2016 is well fit by a double exponential model. In this model, magmatic reservoir overpressure increases at a constant rate over 8.7 years (from October 1997 until August 2006), after which overpressure is held constant. Importantly, this model predicts decreasing uplift over the next several decades at which point uplift should stop. Several mechanisms can be invoked to explain the rate change observed at Lazufre between 1997 and 2016, but with geodesy data alone it is not possible to disentangle the roles of time-dependent pressurization and host-rock rheology. Continued GPS measurements, combined with more frequent InSAR observations with Sentinel-1 and other satellites, as well as contemporaneous analysis of gas emission and seismic activity will help determine the physical mechanism behind the observed deceleration at Lazufre.

ACKNOWLEDGMENTS

We thank members of the PLUTONS team and Jillian Pearse for critical discussions, and Dominique Remy for sharing data files from his 2014 study. We also thank two anonymous reviewers for their insightful feedback. This work was supported by NASA grants NNX08AT02G and NNX10AN57H issued through the Science Mission Directorate’s Earth Science Division, and National Science Foundation (NSF) grant EAR-0908281, which is part of the PLUTONS project. Some of the data used are courtesy of the Committee on Earth Observation Satellites (CEOS) Volcano pilot project coordinated by Mike Poland (U.S. Geological Survey) and Simona Zoffoli (Agenzia Spaziale Italiana), and we thank the space agencies of Japan, Germany, Italy, and Canada and the European Space Agency for the data used in this paper along with the Alaska Satellite Facility for distributing some of it. COSMO-SkyMed data were provided courtesy of the Italian Space Agency (ASI), under CSK AO project ID 188. RSAT2 data were provided by the Canadian Space Agency (CSA); MacDonald, Dettwiler and Associates Ltd.; and the Science and Operational Applications Research (SOAR) and CSA-SOAR-ASI programs. GPS equipment and field support were provided by UNAVCO with support from the NSF and NASA under NSF Cooperative Agreement EAR-0735156. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Finally, we thank Sarah Doelger, Jennifer Jay, Hector Toro, Greco Ramirez, and Doug Christensen for fieldwork assistance.

1 Supplemental Material. InSAR methodology detail and additional figures; GPS tables and time series. Please visit http://doi.org/10.1130/GES01441.S1 or the full-text article on www.gsapubs.org to view the Supplemental Material.

Supplementary data