Source Model for Sabancaya Volcano Constrained by DInSAR and GNSS Surface Deformation Observation

Abstract

Sabancaya is the most active volcano of the Ampato-Sabancaya Volcanic Complex (ASVC) in southern Perú and has been erupting since 2016. The analysis of ascending and descending Sentinel-1 orbits (DInSAR) and Global Navigation Satellite System (GNSS) datasets from 2014 to 2019 imaged a radially symmetric inflating area, uplifting at a rate of 35 to 50 mm/yr and centered 5 km north of Sabancaya. The DInSAR and GNSS data were modeled independently. We inverted the DInSAR data to infer the location, depth, and volume change of the deformation source. Then, we verified the DInSAR deformation model against the results from the inversion of the GNSS data. Our modelling results suggest that the imaged inflation pattern can be explained by a source 12 to 15 km deep, with a volume change rate between 26 × 10⁶ m³/yr and 46 × 10⁶ m³/yr, located between the Sabancaya and Hualca Hualca volcano. The observed regional inflation pattern, concentration of earthquake epicenters north of the ASVC, and inferred location of the deformation source indicate that the current eruptive activity at Sabancaya is fed by a deep regional reservoir through a lateral magmatic plumbing system.

1. Introduction

Sabancaya is a 5980 m-high stratovolcano located in the Central Volcanic Zone (CVZ) of the Andes, 75 km northwest of the city of Arequipa, Perú. It is the youngest and most recently active of the three volcanoes of the Ampato-Sabancaya Volcanic Complex (ASVC), which includes Hualca Hualca to the north and Ampato to the south (Figure 1). ASVC’s volcanoes are mostly composed of andesitic and dacitic lava flows and pyroclastic deposits [1], surrounded by an extensive system of active faults and lineaments (Figure 1). Ampato is a dormant stratovolcano with no historical activity, consisting of three volcanic cones overlying an older eroded volcanic edifice [2]. The Pleistocene Hualca Hualca volcano is believed to be extinct. There is hydrothermal activity observed near Pinchollo (Hualca Hualca), but it could be of tectonic origin [3].

Sabancaya has had five historical eruptive episodes recorded since 1750 CE [4]. Between 1985 and 1997, it produced several Vulcanian eruptions with ash columns 2–3 km high and several small surges and pyroclastic flows [1,5,6]. Crustal deformation related to this eruptive cycle has not been well characterized because of a lack of data. The volcano has not been instrumented yet and very few Synthetic Aperture Radar (SAR) scenes suitable for Differential Interferometry of Synthetic Aperture Radar (DInSAR) studies were acquired by the European Space Agency’s ERS Mission over South America. However, Pritchard and Simons [7] were able to create three useful interferograms with scenes acquired between July 1992 and October 1997. These scenes showed an uplift pattern centered near Hualca Hualca, circular in shape and with a radius of ~20 km. The deformation rate reached a mean velocity of 20 mm/yr at the point of maximum displacement. Reference [7] suggests that the inflation rate was constant during the analyzed time span. Unfortunately, the poor temporal resolution of the data prevented linking the observed deformation to a specific eruptive cycle.

In 2013, Sabancaya experienced a significant increase in seismic activity. Reference [8,9] studied the crustal deformation between 2002 and 2013 by applying the DInSAR processing to data from ERS, ENVISAT, and TerraSAR-X missions. They observed several co-seismic deformation episodes and creep along the Solarpampa fault, interpreting them as tectonic in origin. The regional circular pattern previously detected in 1992 to 1997 was absent.

The present eruptive cycle of Sabancaya began in November 2016 and remains ongoing through May 2020, when this paper was accepted. It is characterized by phreatic and Vulcanian explosions associated with constant SO₂ and ash emissions, accelerated lava dome growth, and thermal anomalies [10]. Seismicity during this period is concentrated mainly north of Sabancaya, around Hualca Hualca (Figure 1), similar to previous distributions [9]. AN analysis of the Sentinel-1, TerraSAR-X, and COSMO-SkyMed missions from 2013–2018 found inflation northwest of Sabancaya as well as creep and rupture on multiple faults [11]. Modelling of the inflation source showed that the location (~7 km N of Sabancaya) and depth (~15 km) were consistent with the source inferred by [7].

In this work, we present results relevant to the present volcanic unrest (October 2014 to March 2019). The DInSAR processing of an ascending and descending Sentinel 1 dataset and Global Navigation Satellite System (GNSS) data analysis shows a radially symmetric inflation similar to the one observed for 1992–1997 by [7]. The observed deformation field is compatible with a deep source, 12–15 km below the surface and 5 km north of Sabancaya, that has possibly been feeding the current Sabancaya eruptions.

2. Deformation from DInSAR and GNSS

2.1. Synthetic Aperture Radar Data

The SAR dataset consists of 84 ascending (Path 47, Frame 1125) and 79 descending orbit (Path 25, Frames 614-646) scenes acquired by the Sentinel-1 Mission (European Space Agency, C-Band) between October 2014 and March 2019. The acquisition mode is Interferometric Wide (IW). Sub-swath 2 of the ascending orbit and sub-swath 1 of the ascending orbit cover the whole area of interest (Figure 1). They were processed using our own implementation of the DInSAR Small Baseline Subsets (SBAS) time-series approach [12], which allows a proper spatial and temporal characterization of the deformation patterns occurring within the studied area.

We treated ascending and descending orbit scenes independently and computed deformation time series and mean deformation velocity maps. We processed seven bursts of the sub-swath 2 for the ascending orbit and seven bursts of the sub-swath 1 for the descending orbit, and analyzed 254 ascending and 232 descending differential interferograms with maximum spatial and temporal baselines of 163 m and 1374 days, respectively. Multilooking factors of 10 and 40 were applied in the azimuth and range directions, respectively, giving a pixel of roughly 150 × 150 m. The topographic fringes were removed via a 30 m Digital Elevation Model from NASA Shuttle Radar Topography Mission (SRTM DEM) [13].

We detected a significant contribution of topography-related fringes in most of the interferograms, attributable to a stratified atmospheric phase component. This is particularly noticeable in the Cañon del Colca region, a 1700 m-deep canyon to the North and West of the ASVC. To compensate for the effect, which can cause a sinusoidal signal in the time series [14], we corrected each differential interferogram using the Zenith Time Delay (ZTD) corrections provided by Newcastle University through the Generic Atmospheric Correction Online Service (GACOS) [15] before a phase unwrapping and time-series calculation.

Because of the contribution from the atmosphere, the delay anomalies, the spatial averaging used to down-sample the DInSAR data, and the structure of the noise in the DInSAR data estimate, it is difficult to specify in a rigorous way the uncertainty in displacement measurements made with DInSAR. On the basis of our own experience and consultations with colleagues who specialize in DInSAR, we assigned an uncertainty to the DInSAR displacement measurements reported here of ±10 mm.

2.2. GNSS Data

The GNSS volcanic deformation monitoring network, set up in ASVC by the Volcano Observatory of the Instituto Geológico Minero y Metalúrgico (OVI-INGEMMET), has been collecting continuous data since 2014. The network consists of five continuous GNSS stations that have been progressively installed since October 2014. Three stations (SBVO, SBSE, and SBHN/SBHO) are within a 4 km radius from Sabancaya, whereas the others (SBMI and SBMU) are 9 to 12 km away from the volcano. SBHN and SBHO are the same station, but the name was changed in October 2016. Not all the stations were operative during the time span covered by the DInSAR time-series (see

Table 1. GNSS stations with their corresponding data acquisition time.

GNSS Station	Start Date	End Date
SBVO	25 September 2014	25 October 2015
SBSE	3 October 2015	31 October 2018
SBMU	24 September 2014	31 October 2018
SBMI	5 October 2016	31 October 2018
SBHO	6 October 2016	31 October 2018
SBHN	26 September 2014	24 September 2016

).

We used GAMIT/GLOBK v.10.70 [16] for processing the GNSS data. The carrier phase and other observable data from satellites visible above a 10° elevation mask were sampled at 30 s intervals. The daily position solutions were computed using precise ephemerides provided by the International GNSS Service (IGS). The initial velocity solution was computed in the South America (SOAM) reference frame defined by [17]. The position time series and preliminary velocities were adjusted to 14 IGS stations located in the Stable South American (SSA) Plate, like POVE (Porto Velo, Brazil), UNSA (Salta, Argentina), KOUR (Kourou, French Guiana), and MTV2 (Montevideo, Uruguay), among others. The uncertainties of the GNSS velocities were derived by scaling the formal error by the square root of the residual chi-square per the degrees of freedom of the solution (

Table 2. Summary of the calculated GNSS velocities per pre-eruptive and eruptive stage (before and after the beginning of the eruption of Sabancaya). The values are in millimeters per year (mm/yr).

		SBVO	SBMU	SBHN/SBHO	SBSE	SBMI
10/14–10/16	E	5 ± 2	−24 ± 2	11 ± 1	2 ± 2
	N	−21 ± 2	−7 ± 1	−10 ± 1	−16 ± 1
	U	25 ± 6	33 ± 4	34 ± 3	27 ± 7
Sabancaya Eruption (6 November 2016)
10/16–10/18	E		−28 ± 1	9 ± 1	1 ± 1	10 ± 2
	N		−7 ± 1	−11 ± 1	−15 ± 1	14 ± 1
	U		26 ± 6	29 ± 3	17 ± 5	26 ± 6

). We also applied additional corrections, considering the monument instability (random walk), instrumental error (white noise), and flicker noise [18].

To eliminate the tectonic component of the deformation velocities observed at Sabancaya, we employed the Euler pole and angular velocity of the Peruvian block. We first estimated the tectonic deformation velocities from historical data from ten GPS stations active in southern Perú between 2012 and 2014. Then, using these velocities, we calculated a Euler Pole [19] located in northern Perú at 6.81° S and 72.52° W, with an angular velocity of 0.62°/My. Finally, the tectonic deformation velocities for each GPS site of the Sabancaya monitoring network (located 1000 km away the Euler Pole) are 11 ± 9 mm/yr east and 1 ± 1 mm/yr north.

Given the stations’ availability (Table 1), the observed time span was split in two periods to improve the outcomes from processing: October 2014 to October 2016, and October 2016 to October 2018. The deformation velocity solutions for the first period include the SBMU, SBVO, SBSE, and SBHN stations, whereas for the second period we have solutions from SBMU, SBSE, SBHO, and SBMI.

2.3. Modelling

We inverted the deformation velocities for GNSS and DInSAR using the dMODELS package [20]. It implements analytical solutions for sources embedded in a homogeneous elastic half space. This approach employs pressurized cavities of simple geometry to mimic/approximate the crustal stress field produced by the actual source. None of these geometries reproduced an actual source. We employed the inversion codes for spherical, spheroidal, and sill-like sources. dMODELS implements a weighted least-squares algorithm combined with a random search grid to infer the minimum of the penalty function, the chi-square per degree of freedom ${\chi }_v^2$.

We inverted the GNSS and DInSAR data independently. This allowed us to verify our modelling results. For DInSAR, we searched for the parameters characterizing each model by jointly inverting the ascending and descending orbits’ deformation velocities. This approach allowed us to better constrain the source geometry, given that the deformation was projected in the Line-of-Sight (LOS, positive for uplift in our images) direction. In order to reduce the computational cost, we down-sampled the original data by a factor of 25 using a linear uniform subsampling method. In this case, given the low spatial gradient of the displacement, we felt that it was not necessary to employ more complex algorithms, such as gradient-based ones [21]. The inversion was performed using 17,640 measurements (pixels).

3. Results

Figure 2 displays the results obtained by DInSAR processing. The mean deformation velocity in the descending (Figure 2a) and ascending (Figure 2b) orbit reveals a deformation field centered on Hualca Hualca volcano but affecting Ampato and Sabancaya volcanoes as well. Its shape is almost circular, with a diameter of ~45 km, and it appears skewed in the opposite direction within Figure 2 because of the SAR lateral view. The flight and LOS direction are indicated in the panels of Figure 2. The LOS-projected maximum deformation is observed in the ascending orbit and reaches a cumulative displacement of 20 cm.

In the descending orbit (Figure 2a,c), the significant shift observed in the first months of 2017 is related to two earthquakes of magnitude Mw = 3.8 and Mw = 3.4 which occurred on January 15 and 30 April 2017, respectively. The time series at R3, the pixel nearest to the epicenters, shows an offset of almost 2.5 cm of LOS decrease linked to the first earthquake. The other two pixels (R1 and R2), located at the centers of the LOS-projected regional volcanic deformation patterns, show shifts linked to the second earthquake. In the latter case, the fringes suggest a LOS increase around 1 cm.

The series taken from the ascending orbit do not present such a clear interpretation due to the low data availability (no acquisition) between 2014 and mid-2016, resulting in a low temporal resolution. The earthquake signals are not as strong as in the descending orbit, probably because of the relative orientation between the deformation and the sensor flight direction. However, minor shifts coincident with the 30 April 2017 earthquake are observed.

The earthquakes’ signals are superimposed on the time series to the linear trend of the regional volcanic deformation. The regional volcanic deformation would be almost linear during the analyzed time span and is perturbed by the two rupture events shown in Figure 2c,d. The deformation mean velocity before and after the earthquakes is very similar, about 35 to 50 mm/yr. The pixels near the epicenters—R3 in Figure 2c—also show a year-long post-seismic relaxation effect followed by a progressive velocity increase.

Figure 3 and Table 2 display the mean deformation velocity (red vectors in Figure 3) computed from the GNSS data. The horizontal components, excluding station SBMI, show a radial deformation centered near Hualca Hualca between October 2016 and October 2018, affecting the whole Ampato-Sabancaya volcanic complex. SBMI is located between the epicenters of the two earthquakes in 2017, so this velocity vector is affected by co- and post-seismic displacement. The vertical components are positive, varying between 1.7 and 3.4 cm/year. These results are consistent with the deformation measured by DInSAR (Figure 3c–g).

Modelling

Table 3. Parameters of the best-fit sources for each period for both data sets for all geometries. Depth is relative to Sabancaya’s vent.

Period	Geometry	Data	χ2v	Lat °	Long °	Depth km	Radius km	δP	δV ×106 m3
10/2014–12/2016	Sphere	DInSAR	1.0	−15.7225 ± 0.0027	−71.8727°	15 ± 1	1.5	0.004 ±0.001	43 ± 5
10/2014–12/2016	Spheroid	DInSAR	6.6	−15.7407	−71.9345°	15.1	1.5	0.009	32
10/2014–12/2016	Sill	DInSAR	6.6	−15.7438	−71.8468°	15.6	10.9	0.0001	26
09/2014–11/2016	Sphere	GNSS	2.8	−15.7335 ± 0.0027	−71.8648°	12 ± 1	1.5	0.002 ±0.001	26 ± 4
09/2014–11/2016	Spheroid	GNSS	1.0	−15.7391	−71.8624°	9.3	1.5	0.481	20
11/2016–11/2018	Sphere	GNSS	25.5	−15.7340 ± 0.0027	−71.8618°	11 ± 1	1.5	0.002 ±0.001	22 ± 4
11/2016–11/2018	Sill	GNSS	17.9	−15.7306	−71.8680°	1.2	1.5	0.052	18
01/2018–4/2019	Sphere	DInSAR	1.7	−15.7159 ± 0.0027	−71.8563°	15.3 ± 0.9	1.5	0.005 ±0.001	51 ± 4

shows the inversion results using both the GNSS and the DInSAR data. To avoid the local effect of the earthquakes (e.g., post-seismic relaxation), we divided the deformation from the DInSAR time series into two periods. The first period is from October 2014 to December 2016, the second from January 2018 to March 2019. The GNSS data have also been split into two time spans linked to station availability (Table 2).

The results for both the time spans are detailed in Table 3. The independent inversion of the GNSS and DInSAR deformation data infers a deep source close to Hualca Hualca as the most probable result in all the modelled geometries.

The sill-like geometry, more appropriate for local deformation patterns produced by shallow sources, yielded parameters geologically implausible (e.g., extremely large radius or near surface depths), so we discarded it.

Of the remaining geometries, the spherical source model better fits both periods with GNSS and DInSAR data. Since the misfit is very similar, we prefer the simpler spherical source (four parameters) to the spheroidal source (seven parameters). The locations from the four best-fit inversions are within a 2 km radius area beneath Hualca Hualca. The dimensionless pressure change is well within the validity range of the linear elasticity assumption of the analytical model (0–0.1). The best fit models employing a spheroidal geometry infer a deeper source (~16 km) and produce a greater geographic dispersion between the best-fit results of the GNSS and DInSAR data.

The best fit for the independent inversion of GNSS and DInSAR data for the first period has a 12–14 km deep and 1.5 km-radius spherical source located beneath Hualca Hualca. Its volume change rate is between 26 × 10⁶ m³/yr and 43 × 10⁶ m³/yr, the dimensionless pressure change between 0.002 and 0.004, and the χ²_v between 1.0 and 2.8.

The independent inversion for the second GNSS and DInSAR deformation velocity series yields similar results (Table 3). The spherical source is located at nearly the same position and depth (12–15 km) for both data sets. Despite a much higher χ²_v for the GNSS inversion, the volume change rate is of the same order of magnitude (see Table 3). We calculated the errors for the best-fit parameters using the Monte Carlo simulation technique [22]. We added white noise to the original data set and then inverted it for the best fit solution. We repeated this process 100 times. The errors presented in Table 3 reflect the distribution of the solutions found.

Figure 3 presents the horizontal velocity vectors produced by the GNSS best-fit spherical sources (blue vectors). Note how the SBMI station shows a noticeable difference between the observed and modelled velocity. This station was affected by the 2017 seismic events, so it recorded the combined effect of the underlying regional pattern and the local earthquakes’ deformation field. Using data provided by the SBMI station results in a higher χ²_v value (~25) and a smaller fit. The modelled and observed data present a significantly better fit at the other stations.

The pattern and residuals obtained by modelling the DInSAR data are presented in Figure 4. After subtracting the modeled component, the possible seismic effects are observed in the residual from both the analyzed periods and are especially noteworthy in 2014–2017. The main seismic deformation pattern located NE of Sabancaya volcano could be associated with a post-seismic effect related to the Mw = 5.9 earthquake on 17 July 2013 [9].

4. Discussion

The analysis and interpretation of the DInSAR and GNSS deformation data allow us to better understand the Sabancaya volcanic system. In particular, our analysis shows that ground deformation affects the three volcanic edifices of the Ampato-Sabancaya volcanic complex (Figure 2 and Figure 3). Deformation patterns and best-fit sources (Figure 5) are centered more than 5 km away from the eruptive vent, close to Hualca Hualca. This lateral magmatic deformation zone is a common feature of volcanic unrest. According to a survey by [23], 24% of monitored volcanoes have had similar deformation signals.

A linear deformation trend is observed throughout the entire time span. This trend is caused by magmatic deformation disturbed by local seismicity (Figure 2). To minimize the influence of seismic deformation on our models, we split the DInSAR and GNSS velocity series into two periods—before and after the 2017 seismic events (Figure 2). The offset and long-term trends observed in the time series after the 2017 seismic events (Figure 2) are the effect of a rapid co-seismic deformation followed by a possible slower post-seismic relaxation.

The deformation source inferred by the independent inversion of the DInSAR and GNSS deformation velocities shows a deep source (12–15 km below the surface) beneath Hualca Hualca with a volume change rate of 26–43 × 10⁶ m³/yr, in agreement with the results obtained by [24] for the same period. The residuals in the DInSAR data are related to seismic deformation and provide a clear view of the faults’ movement. The model fits very well also with the deformation velocity vectors from the GNSS (Figure 3a,b). The station SMBI is affected by the earthquakes and has a combined signal from the local volcanic source and the regional seismic deformation.

The results obtained by the inversion of the GNSS data validate the DInSAR best-fit model, with both having similar locations, depths, and volume change rates. Furthermore, the similar results for both the deformation velocity series reinforce the linear trend hypothesis for the whole period (Figure 3c–g).

The best-fit source, located between 12 and 15 km beneath Hualca Hualca, reproduces the stress field and volume change of a magma reservoir that might have fed Sabancaya’s eruptions since the mid-1990s.

The deep intrusion of magma into an existing reservoir can cause the migration of magmatic fluids into the ASVC plumbing system and a change in the local heat flow [10]. The volume/pressure change started a period of volcanic unrest with inflation accompanied by seismicity. The epicenters of the earthquakes produced by fluid movements along permeable pathways and fault structures concentrate north of Hualca Hualca. Magma and hydrothermal fluids, driven by convection, can reach Sabancaya’s conduit (Figure 5).

Our results confirm the existence of a previous inflation episode in the 1990′s [7], with ground deformation centered on Hualca Hualca. In our opinion, Sabancaya is fed by a laterally extensive magmatic plumbing system (Figure 5) which allows fluids to move along a highly fractured (Figure 1) and pressured crust from the “offset” deep magma chamber to their emission point.

5. Conclusions

Continuous GNSS measurements and DInSAR observation from the ascending and descending Sentinel-1 tracks show a large inflation pattern (~45 km diameter) affecting the ASVC. The beginning date of the inflation is uncertain. The analyzed deformation velocity series have an almost linear trend for the whole period, interrupted by local seismic events. To minimize the local effect of the earthquakes on modelling, we divide both the DInSAR and GNSS data into two series: one before the 2017 earthquakes (October 2014 to December 2016) and one after the post-seismic relaxation (January 2018 to March 2019). For both datasets, the best-fit source is a deep intrusion (12 to 15 km below the surface) beneath Hualca Hualca. This is the same source inferred by [7] for the unrest experienced by Sabancaya in the mid-1990s. The regional inflation pattern centered on Hualca Hualca, the “offset” magma chamber, and the concentration of seismic activity north of ASVC suggest that a laterally extensive magmatic plumbing system may be responsible for this deformation scenario.