Uturuncu volcano in southern Bolivia last erupted around 250 ka but is exhibiting signs of recent activity, including over 50 yr of surface uplift, elevated seismic activity, and fumarolic activity. We studied the spatial and temporal scales of surface deformation from 1992 to 2021 to better understand subsurface activity. We tracked Uturuncu’s recent deformation using interferometric synthetic aperture radar (InSAR) data and the global navigation satellite system (GNSS) station UTUR, located near Uturuncu’s summit. We observed a spatially coherent signal of uplift from 2014 to 2021 from Sentinel-1 A/B satellites that indicates the Altiplano-Puna magma body, located 19–24 km below ground level, and previously noted as the source of the large region of deformation, is still active. The ground is now uplifting at a rate of ~3 mm/yr compared to prior rates of ~10 mm/yr. We corroborated this waning uplift with in situ data from station UTUR. We combined the Sentinel-1 data with TerraSAR-X interferograms to constrain an ~25 km2 region of subsidence located 11 km SSW of Uturuncu, with a source depth of 2.1 km below ground level to an active period of ~2.5 yr with ~5 mm/yr subsidence. We developed a conceptual model that relates these varying depths and time scales of activity in a transcrustal magmatic system. We associate the surface uplift with pressurization from ascending gases and brines from magmatic reservoirs in the midcrust. We infer the existence of brine lenses in the shallow hydrothermal system based on low subsurface resistivity correlated with surface subsidence.

Uturuncu volcano in southern Bolivia last erupted around 250 ka but has been subject to numerous studies over the last two decades due to signs of recent activity, including over 50 yr of surface uplift (Gottsmann et al., 2018), elevated seismic activity, and fumarolic activity (e.g., Sparks et al., 2008; Jay et al., 2012; Hudson et al., 2022). We contribute to these studies by investigating recent surface deformation at Uturuncu to better determine the origin and extent of unrest.

There is growing evidence that many magmatic systems are transcrustal in nature (Cashman et al., 2017). Physical models point to the transient nature of shallow magma chambers, and geochemical data show evidence for long (>103–105 yr) residence times before more rapid mixing and eruption events (Sparks and Cashman, 2017; Costa and Morgan, 2010; Hawkesworth et al., 2004). Cashman et al. (2017) proposed a theoretical model that described vertically stacked pockets of crystal mush extending from the upper crust to the seismic Moho. This conceptual framework allows us to consider the variety of ways that magma, crystals, and volatiles can interact in the subsurface and provides a foundation for analyzing real-world examples of these systems. Uturuncu volcano is one such location that has shown evidence of activity from the lower crust to the surface (e.g., Chmielowski et al., 1999; Comeau et al., 2016; Lau et al., 2018).

Located in the Altiplano region of the southern Bolivian Andes, Uturuncu is also centrally located within the surface footprint of the Altiplano-Puna magma body (Fig. 1). The Altiplano-Puna magma body is the largest known active magma body in the continental crust at ~500,000 km3 (e.g., Chmielowski et al., 1999; Ward et al., 2013), and it is thought to have been the source of large ignimbritic flare-ups in the region over the last 10 m.y., as part of the Altiplano Puna volcanic complex (e.g., de Silva, 1989; Salisbury et al., 2011). Many studies have been conducted in recent years to investigate the subsurface structure at Uturuncu, including magnetotelluric (Comeau et al., 2015, 2016), gravity (del Potro et al., 2013; MacQueen et al., 2021; Gottsmann et al., 2022), and seismic surveys (Ward et al., 2013; Kukarina et al., 2017; Hudson et al., 2022), many of which were summarized by Pritchard et al. (2018). We sought to better understand the subsurface structure by looking at both the spatial and temporal scales of surface deformation in recent years.

Uturuncu’s recent deformation has been tracked by two main methods: interferometric synthetic aperture radar (InSAR) and a global navigation satellite system (GNSS) station, abbreviated UTUR, located on the northwest flank of Uturuncu’s edifice (67.2055°W, 22.2420°S, 5184 m). Starting with the earliest available InSAR data, European Remote Sensing (ERS-1/2) satellite acquisitions in 1992, Pritchard and Simons (2002) first noted a 70-km-diameter region of surface uplift centered on Uturuncu rising at ~1 cm/yr. Fialko and Pearse (2012) and Henderson and Pritchard (2013) confirmed the continuation of surface uplift and noted the presence of an outer ring of subsidence that was sinking at a rate of a few millimeters per year. Various models have been employed to explain the surface deformation, including a magma mush column (Gottsmann et al., 2017), an ascending diapir (Fialko and Pearse, 2012), the injection of a buoyant fluid at the brittle-ductile transition (Morand et al., 2021), and other source types summarized in Barone et al. (2019). The various models have different implications for the volcanic system but tend to agree on a source depth of 19–24 km below ground level for the deformation signal. GNSS time series through 2015 and InSAR time series through 2017 have been reported by Henderson and Pritchard (2017) and Lau et al. (2018), where both groups noted slowing surface uplift rates at Uturuncu. Additionally, Lau et al. (2018) recorded the existence of a localized zone of subsidence within the original region of uplift.

In the following sections, we present new results from 2014 to 2021 from synthetic aperture radar (SAR) satellites Sentinel-1 A/B to compare with earlier InSAR results. We observed a spatially coherent signal of surface uplift that indicates the source of deformation 19–24 km below ground level is still active, although rates have slowed. We corroborated this waning uplift with in situ data from GNSS station UTUR, extending the time series from 2015 to 2020. Furthermore, we constrained the length and rate of activity for a small region of subsidence first noted by Lau et al. (2018) using TerraSAR-X data spanning July 2012 to December 2018. We compared this deformation data set with other geophysical subsurface data and developed a conceptual model that relates these varying depths and time scales of activity to a release in pressure due to the waning deep uplift. We found that the data support the model proposed by Gottsmann et al. (2017), wherein a magma mush column with brines and gases ascends from the Altiplano-Puna magma body, and we posit these materials collect in brine lenses in the shallow subsurface.

InSAR

Sentinel-1A/B

InSAR is a widely used technique that measures relative surface displacements along radar line of sight (LOS), and it has greatly expanded volcano monitoring capabilities (e.g., Biggs and Wright, 2020). SAR satellites such as European Remote Sensing 1/2 (ERS-1/2) and Envisat (both C-band radars, ~5.6 cm) have been recording deformation at Uturuncu since as early as 1992. Henderson and Pritchard (2013) and Fialko and Pearse (2012) independently processed and analyzed strip-map InSAR data spanning May 1992–January 2011 from the ERS-1/2 and Envisat satellites. Lau et al. (2018) presented Sentinel-1 A/B data (also C-band) from 2014 to 2017. We processed available Sentinel-1 A/B data from 2014 to 2021 using four tracks: two ascending (A76, A149) and two descending (D83, D156) (Table 1). For each track, we created a network of month-long and year-long interferograms using the InSAR Scientific Computing Environment (ISCE; Rosen et al., 2012) on the cloud-computing resource hosted by Amazon Web Services (Henderson and Setiawan, 2020). We used the 30-m-resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) to correct for topographic effects in the interferograms (Farr et al., 2007). We chose a reference region in the north of each track’s data coverage due to its distance from Uturuncu and relatively low noise levels (Fig. 2).

The interferometric phase was influenced by atmospheric delays (Bevis et al., 1992) in addition to the deformation signal we sought to constrain. Turbulent tropospheric signals are frequently more spatiotemporally variable than ionospheric signals, which makes them more difficult to correct for (Rosen et al., 2018; Zebker et al., 1997). Ionospheric and stratified tropospheric signals can be estimated and removed. We used ISCE to correct for ionospheric effects, although this correction did not remove all such effects. We saw more noise retained in the ascending Sentinel-1 tracks acquired in the early evening, when ionospheric effects would be larger than the descending tracks near the magnetic equator (Liang et al., 2019).

Common ways to remove tropospheric effects are through data-based weather models or theoretical models (e.g., Bekaert et al., 2015). We tested three methods of correcting for tropospheric noise to determine the one best suited for the data: (1) An empirical linear relationship between phase and topography was removed for the whole region (Lin et al., 2010), (2) European Reanalysis v5 (ERA5) corrections were applied to the whole region (Hersbach et al., 2020), and (3) Generic Atmospheric Correction Online Service (GACOS) corrections were applied to the whole region (Yu et al., 2018a, 2018b). GACOS corrections resulted in the largest reduction in phase variance among the tested methods; therefore, GACOS corrections were used on all four tracks for processing (see Introduction and Figure S1 in Supplemental Material1 for more details).

We used the Generic InSAR Analysis Toolbox (GIAnT) and a modified small baseline subset (SBAS) method to generate a time series of ground surface movement over the 6 yr span for each pixel above a coherence threshold of 0.6 (Agram et al., 2013; Berardino et al., 2002). Uturuncu’s location in the Altiplano means that both arid conditions and low vegetation led to high temporal coherence for most regions in the interferograms (e.g., Henderson and Pritchard, 2013).

TerraSAR-X Data

To partially fill the gap between the end of the Envisat InSAR data in January 2011 and the start of Sentinel-1 data in November 2014, we created 11 interferograms with TerraSAR-X (TSX) data spanning 13 July 2012–27 October 2014 (Table 1) using the ISCE software (Rosen et al., 2012). Interferogram pairs were chosen from a set of eight SAR scenes such that perpendicular baselines were <200 m. Interferograms were processed with eight looks in azimuth and range, using 30 m SRTM DEM (Farr et al., 2007) upsampled to 10 m resolution. Interferograms were unwrapped using the snaphu-mcf algorithm (Chen and Zebker, 2001).

To enhance the detection of small deformation signals and decrease the effects of atmospheric noise in the data, we constructed a stack of seven interferograms by adding together all the interferograms spanning those dates. Before adding each interferogram to the stack, we masked areas below a coherence threshold of 0.4, and we simultaneously fit a bilinear ramp and a linear trend between elevation and the unwrapped data and removed these fitted signals. This gave us unwrapped data with reduced effects from the stratified atmosphere.

Uncertainties

To quantify uncertainty between Sentinel-1 tracks, we followed the process laid out in Finnegan et al. (2008). We compared the velocity calculations from overlapping data points in both the ascending and descending tracks. Figure S2 shows the comparisons between the velocities of the ascending (Figs. S2A and S2B) and descending (Figs. S2C and S2D) tracks along with their histograms of distance from the 1:1 line. We expected the velocity values to differ by ~10% due to the different lines of sight between tracks. To ascertain the respective uncertainties on the velocities, we fit a Gaussian to the histograms and found standard deviations of ~0.4 mm/yr for the ascending tracks and ~0.2 mm/yr for the descending tracks. We expect that these uncertainties are representative of the whole study region because the overlapping regions span a range of elevations. We used these values to quantify the precision of the Sentinel-1 time series.

For the TerraSAR-X error analysis, we did not have enough data to compare overlapping data points. Instead, we calculated an approximate noise level for the stack of seven interferograms by taking the standard deviation of all pixels in the final stack, excluding the small area of potential deformation. To determine the noise level of a single interferogram, we took the standard deviation of all pixels in the interferogram, excluding the area of potential deformation. Both the stack and the single interferogram discussed in the Transient Signals subsection of the Results section have uncertainties of ~1 mm/yr.

GNSS

There are two continuous GNSS stations installed near Uturuncu: COLO and UTUR (Henderson and Pritchard, 2017), with 8.4 and 10.25 yr of data, respectively (Fig. 1). All GNSS data were obtained from the University of Nevada at Reno Nevada Geodetic Laboratory (NGL) (http://geodesy.unr.edu; Blewitt et al., 2018). The data consist of daily acquisitions that have been referenced using the IGS14 reference frame. Henderson and Pritchard (2017) reported on the first 5.5 yr of data, from 10 April 2010 to 24 November 2015. We reproduced their analysis with the updated data sets for both stations with COLO data through November 2018 and UTUR data through July 2020, when the station ceased operation.

Model Fits

We used three different models (linear, multi-parameter, and MIDAS) to fit horizontal and vertical components of the data for UTUR and COLO. We used an additional two models (quadratic and piecewise) to fit the vertical data for UTUR. The linear model included a standard linear fit with a step on the date of the 2014 Pisagua, Chile, earthquake. The linear rate was assumed to be constant before and after the earthquake in this model. The multiparameter fit consisted of a constant linear term, a yearly seasonal term, and steps for earthquakes in the NGL catalog. MIDAS, or median interannual difference adjusted for skewness, is a fit developed by Blewitt et al. (2016) that is particularly robust at avoiding bias due to offset in the time series because it is designed to account for steps and periodic annual signals.

We also investigated if there are models that vary the uplift rate over time that better fit the data. Henderson and Pritchard (2017) reported that a quadratic fit of the vertical component from April 2010 to November 2015 had fewer functional constraints than the multiparameter fit and a similar root mean square error (RMSE). We performed this analysis again using a quadratic fit on the vertical UTUR data through 2020. We also performed a piecewise fit by dividing the vertical UTUR time series into three time periods and finding the linear rate of each period. Blewitt and Lavallée (2002) reported that fitting models to a GNSS time series shorter than 2.5 yr can lead to significant bias from seasonal signals. Because of this, each time period was over 3 yr long and started and ended in April to minimize the impact of the seasonal signal. Where possible, the time series was divided on the dates of regional earthquakes as reported in the NGL catalog to minimize the effect of coseismic movement.

Seasonality

In addition to the surface uplift rate analysis, we investigated the regional seasonal patterns with a broader network of GNSS stations that covered the same time span (Fig. 1A) to better understand the origin of the seasonal signal at UTUR. The seasonal signals clearly appear in the vertical and detrended north components (vertical shown in Fig. S3; north not shown). Based on the data from April 2010 to November 2015 and the work of Fu et al. (2013), Henderson and Pritchard (2017) hypothesized that water loading in the Amazon drives the seasonality. We tested this hypothesis with data from April 2010 to April 2018 (the longest time period for which all stations have data) at five stations (Fig. 1; Fig. S4). Continuous GNSS stations in the Altiplano are sparsely located, so these seasonal comparisons were limited to five stations. To isolate the seasonal signal, we used the NGL’s event database to locate dates of earthquakes and equipment changes and removed their corresponding phase center offsets. We detrended the data by removing a linear model to minimize any long-term tectonic or magmatic movement (Fig. S3). We found that the vertical seasonality at UTUR correlates with the seasonality of POVE, the station closest to the Amazon Basin. More details about the seasonal signal can be found in Section 3 of the Supplemental Material.

Large-Scale Deformation

Uturuncu’s surface movement in the satellite LOS from 1992 through 2021 is presented in Figure 3 for a point near the UTUR GNSS station. This figure includes the previously processed ERS-1/2 and Envisat data, along with updated Sentinel-1 and GNSS data. The inset emphasizes the historic surface uplift rates in comparison with the current, lower uplift rates. As previously reported in Henderson and Pritchard (2013), the peak rate from 1992 to 2005 was ~10 mm/yr. Since 2005, the surface uplift rate has been tapering off, most noticeably in the more recent Sentinel-1 data. The peak LOS velocities for the Sentinel-1 tracks (which span 2014–2021) are ~2.5 mm/yr. Given the 0.4 mm/yr and 0.2 mm/yr errors for the ascending and descending tracks, this region of surface uplift is above the level of noise.

Linear velocity maps for each track are presented in Figure 2. The faint signal of surface uplift centered on Uturuncu is spatially coherent in the descending tracks, while it is visible but less obvious in the ascending tracks. The ascending acquisitions were taken at dusk, so we would expect to have an increased impact of ionospheric effects in our low-latitude study region (Liang et al., 2019). Furthermore, we note increased phase variance from December through March, the Altiplano’s rainy season (Garreaud et al., 2003). The previously identified ring of subsidence is no longer discernible from noise in our results, as Lau et al. (2018) also noted.

Time series for COLO and UTUR are presented in Figure 4, and bestfit rates are presented in Table 2. The obtained surface uplift rates at UTUR range from 3.2 mm/yr to 3.6 mm/yr depending on the type of curve fitting used (Table 2). We disregarded the quadratic fit because the quadratic term contributed very little to the fit with the additional ~5 yr of data (0.008t2 + 3.59t + 4.03, RMSE = 4.5 mm; t is in days). The fit with the smallest RMSE for all data was the multiparameter fit, with an RMSE of 4.0 mm (0.4 mm/yr) for the vertical component at UTUR. Both the multiparameter and MIDAS fits agreed with the Sentinel-1 uplift rates, within error.

The models of interest (linear, multiparameter, and MIDAS) all assumed constant rates. We split the vertical UTUR time series into three time periods and found the linear rate of each to determine if there were any obvious changes in surface uplift rate. With this, we found a surface uplift rate of 2.8 mm/yr from 2010 to 2014, which increased to 4.7 mm/yr from 2014 to 2017 and then decreased to 3.2 mm/yr from 2017 to 2020 (Fig. 4). These differences in rates are larger than the error bounds for each fit. The GNSS models are further discussed in the GNSS Signal subsection of the Discussion section.

Transient Signal

The Sentinel-1 time series of deformation indicates that the small area of subsidence ~11 km SSW of Uturuncu (Figs. 2 and 5), which was first discussed by Lau et al. (2018), stopped deforming in early 2017. The deformation was already ongoing in October 2014 when the Sentinel-1 time series from Lau et al. (2018) began, so we used TerraSAR-X (TSX) data dating back to 2012 to attempt to better constrain the start date. We created a single TSX interferogram spanning 27 October 2014–12 December 2018 (Fig. 5A), which clearly shows the small subsidence area (dashed circle) south of Uturuncu also visible in Sentinel-1 velocity maps spanning this time period (Fig. 2). The noise level of this interferogram is ~1 mm/yr, and the maximum subsidence reaches ~5 mm/yr. No individual TSX interferograms prior to 27 October 2014 showed the subsidence source. The interferogram stack in Figure 5B spans 13 July 2012–27 October 2014 and shows a small area of possible subsidence in the northwest part of the known subsidence area (magenta star in Fig. 5B).

The deformation rate in this area, however, is comparable with deformation rates in other stable areas of the map and was not consistently observed in our analysis of the interferograms, appearing and disappearing depending on the corrections applied or plotting method. Based on this, we conclude that any deformation in this area, if present, is below the noise level of ~1 mm/yr between 13 July 2012 and 27 October 2014 in the stacked TSX data. The subsidence south of Uturuncu likely began after October 2014, although there are some alternative possibilities considered in the discussion. Using the Sentinel-1 time series, we determined a maximum rate of subsidence of ~5 mm/yr from October 2014 to January 2017, i.e., greater than the uncertainties of 0.4 mm/yr and 0.2 mm/yr determined for the ascending and descending tracks of Sentinel-1 data.

Historic Uplift with Moat of Subsidence

The initial signal that compelled researchers to study Uturuncu in detail was a large, 70-km-diameter surface uplift of ~1 cm/yr, surrounded by a moat of subsidence of a few millimeters per year (Pritchard and Simons, 2002; Fialko and Pearse, 2012). This deformation pattern persisted from May 1992 to January 2011 (e.g., Henderson and Pritchard, 2017; Lau et al., 2018). Several theories have been proposed to explain this perplexing signal, including an ascending diapir (Fialko and Pearse, 2012), buoyant fluid stored at the brittle-ductile transition (Morand et al., 2021), a magma mush column (Gottsmann et al., 2017), and various depths of pressurized magma chambers (e.g., Pritchard and Simons, 2002; Hickey et al., 2013; Barone et al., 2019). While each of these theories explains the spatial extent of the surface uplift, they have different implications for the expected deformation over time. With the updated deformation time series and velocity maps, we can better inform our conceptual model of the subsurface dynamics.

From 2014 to present, there has been a spatially coherent uplift signal within the footprint of the original 70-km-diameter uplift. The surface uplift rate has dropped to 2–3 mm/yr, and the descending tracks agree within error with the GNSS station UTUR (Table 3). The Sentinel-1 data have been referenced to a common reference point within the study region, while the GNSS station uses the IGS14 reference frame, so we do not expect their agreement to be exact. The ring of subsidence surrounding the central uplift has either ceased deforming or is under the noise threshold that we reasonably expect to see in these data. Our results indicate a slowing of deformation. Gottsmann et al. (2017) proposed that the deformation is a result of episodic mush reorganization. If we are nearing the end of an episode of mush reorganization, the volume of fluids and gases percolating to the surface will lessen, leading to a decrease in the pressurization rate. Over time, we expect the surface deformation would reverse with stress relaxation. The current slowing deformation may indicate an upcoming pause and reversal; further monitoring is necessary to track this interpretation. We additionally compared the current (2014–present) radial velocity profiles from the center of the uplift to the 1992–2011 velocities published in Henderson and Pritchard (2017) (see Fig. S5 herein). The current radial velocity profiles are less steep than they have been in the past, indicating a change in the pressurization rate, and potentially revealing a decrease in subsurface fluid migration.

Small-Spatial-Scale Deformation

Lau et al. (2018) first noted a shallow region of subsidence ~11 km SSW of Uturuncu. With the TSX data, we were able to put a lower bound of 2.5 yr on the duration of subsidence, with a velocity of 5 mm/yr from October 2014 to January 2017. The subsidence end date was determined using the time series in Figure 6B. The deformation could have potentially begun at lower rates before October 2014, but the limited availability of TSX data makes detection more unlikely. Additionally, the 70-km-diameter uplift signal could have concealed the smaller subsidence signal in months before October 2014. Therefore, 2.5 yr is a lower bound for the signal duration. Short-lived signals have been noted at Uturuncu before; Barone et al. (2019) proposed the existence of an active inflating source at ~4.5 km below sea level from August 2006 through February 2007. The existence of these transient signals of deformation indicates activity at multiple spatial and temporal scales.

We performed a forward Mogi model (Mogi, 1958) for the source of subsidence and found that the source depth of 2.12 km below ground level proposed by Lau et al. (2018) was consistent with our signal, and our changes in volume are comparable within error (this study, dV = 1.80 × 10−4 km3/yr; Lau et al. [2018], dV = 1.95 × 10−4 km3/yr) (see Conclusions section and Fig. S6 in Supplemental Material for more details). We compared the source depth of subsidence (2.12 km below surface) to a resistivity model developed by Comeau et al. (2016) from magnetotelluric data collected from 2011 to 2013. We determined that the area of subsidence corresponds to a low-resistivity zone estimated at ~3 Ω·m (Fig. 7 herein; Comeau et al., 2016). In a volcanic environment such as Uturuncu, low resistivity typically corresponds to partial melt, hydrothermal alteration, or brines (e.g., Evans, 2012; Comeau et al., 2015, 2016). Uturuncu has a history of primarily dacitic eruptions, and the observed low resistivities in this depth range with the petrological constraints on the dacitic melt would require ~100% melt (Pritchard et al., 2018). This is geologically improbable given the amount of time since the last eruption (~250 k.y.).

The presence of brines has been hypothesized for the shallow (2–4.5 km below ground level) region and could be the source of the deformation activity (Pritchard et al., 2018; Comeau et al., 2016). Afanasyev et al. (2018) described the ability of high-salinity brines to form above degassing magma bodies beneath dormant volcanic systems, which may persist for hundreds of thousands of years. According to their model, once degassing ends, the brines will spread laterally and sink back through the formational conduit (Afanasyev et al., 2018). We propose the existence of numerous disconnected lenses forming the low-resistivity anomalies in Figure 7. The ~25 km2 region of subsidence noted from October 2014 to January 2017 was likely a singular lens collapse within a much larger plexus of brine lenses because of the small size of the deformation anomaly compared to the much larger resistivity anomaly, as shown in Figures 6 and 7. It is possible that we do not see additional transient signals either because they were active before 2014, when we had less frequent InSAR acquisitions, or because they were hidden by the larger signal. There are no obvious surface features such as fumaroles that have been noted near the subsidence feature. We expect that the brine lens collapsed diffusely, so we do not anticipate related surface uplift due to brine movement. If the surface uplift rate remains low, then the number of brine lens collapses should increase with time, so we expect more small-scale signals from brine lens collapses in the future.

Conceptual Model of Subsurface Dynamics

The two leading theories for explaining the observed surface deformation at Uturuncu volcano are an ascending magmatic diapir (Fialko and Pearse, 2012) and magma mush reorganization (Gottsmann et al., 2017). There is still spatially coherent uplift of ~2–3 mm/yr within the original 70-km-diameter region of surface uplift, but the surrounding moat of subsidence has either stopped or decreased to below the noise threshold of the Sentinel-1 data. Lau et al. (2018) argued that the slowing deformation rate remains consistent with an ascending magmatic diapir and explained that short-term fluctuations in uplift rate could be related to “local magma intrusions or migration of volatiles” (p. 46). This is not necessarily inconsistent with the Gottsmann et al. (2017) model. As Lau et al. (2018) noted, the main difference in explaining the short-term fluctuations lies in the relative importance of volatiles versus magmas. We argue for the magma mush reorganization model here primarily due to geomorphic evidence showing no evidence of the long-term uplift typically associated with an ascending diapir at Uturuncu (Perkins et al., 2016). The presence of brines also supports the important role of volatile flux in feeding them (Afanasyev et al., 2018).

The magma mush reorganization model is consistent with the available geodetic data. Gottsmann et al. (2017) proposed that, in time, this system would allow for a reversal of deformation patterns, which would explain the lack of long-term deformation in the region (Perkins et al., 2016). A reversal of deformation is not necessarily expected with the ascending diapir model and would be further evidence to determine the preferred model for the region. Continued monitoring using both the InSAR and available GNSS data is necessary to track the waning signal and determine the long-term processes of which we may be observing the effects.

Based on the available InSAR data, there is evidence for multiple levels of deformation. In this study, we show that the 70-km-diameter uplift that started ~50 yr ago has been waning over the past 10 yr. We also show that the localized region 11 km SSW of Uturuncu underwent a short period of subsidence as the larger surface uplift lessened. We propose that these deformation patterns can be connected by extending the magma mush model proposed by Gottsmann et al. (2017) to explain the large-scale “sombrero” signal.

Importantly, the model by Gottsmann et al. (2017) proposes that fluids and gas are exsolving from a magmatic column connected to the Altiplano-Puna magma body, and they form a shallow brine-rich hydrothermal reservoir (Fig. 7). Now that we see evidence for slowing uplift, the proposed conduit could be shutting off, meaning a lack of volatile influx into the shallow hydrothermal system. We propose that this led to a localized brine lens collapse 11 km SSW of Uturuncu. We indicate the conceptualized fluid and melt movement with arrows on Figure 7 and distinguish between the movement happening on a decadal (or longer) time scale (blue and gray arrows) versus a subdecadal time scale (purple arrows). The waning input from the deep system predicts that there will be a deeper subsidence signal similar in footprint to the 70-km-diameter historic uplift according to Gottsmann et al.’s model. In addition, more shallow subsidence features similar to the inferred subsidence due to brine lens collapse should be observed in the future.

GNSS Signal

Past work indicates that the surface uplift rate at Uturuncu fluctuated over the time period studied (Henderson and Pritchard, 2017). It is important to investigate the functional forms of these temporal variations to understand their underlying causes. Of the five methods of fitting the GNSS time series, the multiparameter fit performed best of the time-invariant fits (linear, MIDAS, multiparameter), and the piecewise linear fit performed best of the time-varying fits (quadratic, piecewise) (Table 2). In accounting for both the seasonal signal and time steps due to events such as earthquakes, the multiparameter fit minimized the RMSE. The quadratic fit proposed in Henderson and Pritchard (2017) no longer explains the updated data set through 2020. It is more likely that the shorter time period used in the initial analysis (through November 2015) appeared more quadratic in nature. The piecewise linear function shows changes in uplift rate in mid-2014 and mid-2017 (Fig. 3). One interpretation of the fluctuating uplift rates is that it is in response to the driving force of uplift. The acceleration and deceleration of surface uplift are consistent with the viscoelastic relaxation of the crust in response to an ascending diapir but can also be indicative of pulses of exsolving volatiles. Both the ascending diapir and magma mush models allow for time-varying rates of deformation. Even though the geodetic data do not rule out the ascending diapir theory, we believe that the literature analyzing geomorphology at Uturuncu provides more evidence for the magma mush reorganization theory (Perkins et al., 2016).

The improved spatial and temporal resolution of Sentinel-1 data has allowed us to discern short-duration, small signals, such as the subsidence 11 km SSW of Uturuncu first noted by Lau et al. (2018). This subsidence corresponds spatially with a low-resistivity zone, likely representing brines in the shallow subsurface. The timing of the subsidence corresponds with the waning uplift signal at Uturuncu, giving us reason to believe the two signals are connected. The lessening pressure and decreased influx of volatiles could lead to a shallow brine lens collapse, which would result in a signal such as the one we observed. Additional data sets such as the resistivity models (Comeau et al., 2015, 2016), seismic surveys (Ward et al., 2013; Kukarina et al., 2017; Hudson et al., 2022), gravity surveys (MacQueen et al., 2021), and degassing data will contribute to a more complete understanding of the subsurface dynamics. Continued monitoring is necessary to investigate the time series for additional transient signals of deformation to form a more complete understanding of Uturuncu’s transcrustal magmatic system.

The GNSS station UTUR covers an essential period of time from early 2011 to late 2014, when InSAR data are lacking for this region (although there is limited coverage for part of this interval from TerraSAR-X). These continuous measurements on Uturuncu’s edifice give a more detailed look at the story of surface uplift. While various fits to the vertical time series give uplift rates of ~3 mm/yr from 2010 to 2020, a closer look by dividing the time series into segments shows that the uplift rate could be varying more than initially thought. Pulses of fluids migrating toward the surface could lead to changes in the surface uplift rate because of fluxes in pressurization rate. Decreased pressurization rates can lead to brine lens collapses in the shallow subsurface. This study shows the importance of monitoring volcanic systems that have no recent record of eruption. We can gain insight into what could be the end of a volcano’s life cycle, better understand the changing geophysical signals that are related to noneruptive activity, and track the evolution of a dynamic transcrustal magma system.

1Supplemental Material. Additional details on InSAR processing and analysis, GNSS analysis, and Mogi model. Please visit https://doi.org/10.1130/GEOS.S.21703889 to access the supplemental material and contact editing@geosociety.org with any questions.
Science Editor: David E. Fastovsky
Guest Associate Editor: Gary Michelfelder

We would like to thank the Observatorio San Calixto, La Paz, Bolivia (especially Ruben Tintaya and Gonzalo Fernandez), Faustino Ticona of the Escuela Militar de Ingeniería, La Paz, Bolivia, and the Servicio Nacional de Áreas Protegidas de Bolivia for their assistance with the GNSS time-series data. Elizabeth Eiden, Matthew Pritchard, and Patricia MacQueen were partly supported by National Science Foundation (NSF) grant EAR 1757495, Scott Henderson was partly supported by NSF grant EAR 0908281, and Patricia MacQueen was partly supported by National Aeronautics and Space Administration (NASA) Future Investigators in NASA Earth and Space Science and Technology (FINESST) award number 80NSSC19K1339. We thank Felipe Aron, an anonymous reviewer, and editors Gary Michelfelder and David Fastovsky for their helpful comments which greatly improved this work.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.