# The 2013–2020 seismic activity at Sabancaya Volcano (Peru): Long lasting unrest and eruption

Roger Machacca, Philippe Lesage, Hernando Tavera, Jeremy D Pesicek, Corentin Caudron, Jose L Torres, Nino Puma, Katherine Vargas, Ivonne Lazarte, Marco Rivera, et al.

► To cite this version:

Roger Machacca, Philippe Lesage, Hernando Tavera, Jeremy D Pesicek, Corentin Caudron, et al. The 2013–2020 seismic activity at Sabancaya Volcano (Peru): Long lasting unrest and eruption. Journal of Volcanology and Geothermal Research, 2023, 435, pp.107767. 10.1016/j.jvolgeores.2023.107767. hal-04273592

HAL Id: hal-04273592

https://hal.science/hal-04273592v1

Submitted on 7 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

# The 2013–2020 seismic activity at Sabancaya Volcano (Peru): Long lasting unrest and eruption

Roger Machacca<sup>a,b,c,*</sup>, Philippe Lesage<sup>c</sup>, Hernando Tavera<sup>a,b</sup>, Jeremy D. Pesicek<sup>d</sup>, Corentin Caudron<sup>e</sup>, Jose L. Torres<sup>a</sup>, Nino Puma<sup>f</sup>, Katherine Vargas<sup>a</sup>, Ivonne Lazarte<sup>a</sup>, Marco Rivera<sup>a</sup>, Alain Burgisser<sup>c</sup>

## ARTICLE INFO

### Keywords:

- Sabancaya volcano

- Long-period events

- Repeating earthquakes

- Seismic velocity change

- Crater migration

## ABSTRACT

Sabancaya volcano is the youngest and second most active volcano in Peru. It is part of the Ampato-Sabancaya volcanic complex which sits to the south of the ancient Hualca Hualca volcano and several frequently active faults, thus resulting in complex volcano-tectonic interactions. After 15 years of repose, in 2013, a series of 4 earthquakes with magnitude >4.5 occurred within 24 h, marking the beginning of a new episode of unrest. Several additional swarms of earthquakes occurred in the following years until magmatic eruptive activity started on 6 November 2016. This activity is ongoing as of this writing, with an average of 50 explosions per day. In this study, we present results of multiparametric monitoring of Sabancaya's activity observed during 2013–2020. Seismic data are used to create a one-dimensional seismic velocity model, to catalog, locate, and characterize earthquakes, to detect repeating earthquake families, and to monitor seismic velocity variations by ambient noise cross-correlation. These analyses are complemented by visual and remote sensing observations and ground deformation measurements. All monitored parameters showed significant changes on 6 November 2016, the day of eruption onset, thus dividing the eruptive activity into pre-eruptive and eruptive stages.

The unrest is characterized by high levels of seismic activity with hundreds of events detected per day. Volcano-tectonic (VT) earthquakes were dominant during the pre-eruptive period while long-period (LP) events and explosions have been most numerous since the eruption onset. Earthquake locations highlight long-lasting seismogenic zones along multiple previously active regional faults, as well as along newly identified faults. This VT seismicity is mainly distributed in a sector from the northwest to the east of the volcanic complex at distances of up to 30 km from the crater. We focus our analysis on two eruptive episodes: the eruption onset and subsequent crater migration from south to north, and the increase of lava dome extrusion rate in 2019. Both episodes are accompanied by seismic velocity decreases of up to 0.2% and are preceded by a few weeks by bursts of distal VT activity, including numerous repeating earthquakes. These repeated events were located on several remote tectonic faults (5–25 km from the vent). We suggest that these phenomena could be due to the injection of a batch of magma in the deep reservoir and/or conduit, which would generate 1) a pressure wave propagating in the hydrothermal system, triggering the bursts of seismic activity and 2) slow rising of magma by melting old material filling the conduit that eventually produced the eruptive and dome growth acceleration events.

## 1. Introduction

Sabancaya volcano (Peru) has been erupting since November 2016.

Stratovolcanoes with long-lasting eruptive cycles, such as Sabancaya, pose a significant threat to local populations for three main reasons. 1) They constantly expose local populations to variable amounts of ash that

<footnote><sup>a</sup> Instituto Geofísico del Perú, Observatorio Volcanológico del Sur, Asentamiento Humano José María Arguedas Zona A Mz. D Lt. 08, Sachaca, Arequipa, Peru</footnote>

<footnote><sup>b</sup> Pontificia Universidad Católica del Perú, Av. Universitaria s/n, San Miguel, Lima, Peru</footnote>

<footnote><sup>c</sup> Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERRE, 38000 Grenoble, France</footnote>

<footnote><sup>d</sup> U.S. Geological Survey, Volcano Disaster Assistance Program, Vancouver, WA, USA</footnote>

<footnote><sup>e</sup> Laboratoire G-Time, Department of Geosciences, Environment and Society, Université libre de Bruxelles, Belgium</footnote>

<footnote><sup>f</sup> Universidad Nacional de San Agustín de Arequipa, Arequipa, Peru</footnote>

<footnote>* Corresponding author at: Instituto Geofísico del Perú, Observatorio Volcanológico del Sur, Urb. La Marina B19, Cayma, Arequipa, Peru. E-mail address: roger.machacca@gmail.com (R. Machacca).</footnote>
can cause chronic health problems (Baxter et al., 2014), they pose a risk to the life and livelihood of local communities by damaging infrastructure and displacing populations, and they can have global impacts by disrupting business and air travel (e.g., Brown et al., 2015; Loughlin et al., 2015; Ridley et al., 2014). 2) Continuous accumulation of pyroclastic material on the flanks of the volcanic edifice promotes the formation of rain-triggered lahars (Jones et al., 2015). 3) Local authorities and populations used to consider the frequent low activity of the volcano as a normal activity level and might underestimate the potential impacts of the less frequent but larger eruptions (Mothes et al., 2015).

Many, but not all, volcanic eruptions are preceded by some form of unrest, including increased seismic activity, ground deformation, gas emissions, and thermal anomalies (e.g., Barberi et al., 1984; Newhall, 2000; Potter et al., 2015; Reath et al., 2016; White and McCausland, 2016; Pesicek et al., 2021). Long-term monitoring is the best means of recognizing these signs of unrest and determining whether they are anomalous or not, and particularly if they are indicative of an impending eruption (e.g., Phillipson et al., 2013), or changing hazards during an eruption (e.g., Sparks, 2003; Segall, 2013). However, to determine whether activity is anomalous, a baseline for background behavior must be established beforehand.

During long-lasting eruptions, eruptive dynamics can rapidly shift from low level activity to highly explosive activity that translates into substantial increase of volcanic threat (Hidalgo et al., 2015). The combined analysis of continuous geophysical signals recorded by monitoring networks can give, in near real-time, at least qualitative insights into the volcano's behavior and help to forecast eruptive events (Ripepe et al.,

2002, 2005). Monitoring and accurately interpreting the evolution of eruptive dynamics through geophysical observations ideally requires multiparametric monitoring systems including seismic, acoustic, geodetic, thermal, geochemical, infrared (IR), and visual observations.

Sabancaya was the first volcano monitored by the Instituto Geofísico del Perú (IGP). Seismological measurements date back to the 1990s; the first seismic network became operational during the eruption of 1986–1998. This first network was, however, removed after the eruption. On 23 February 2013, one day after the beginning of a new unrest episode, temporary stations were deployed. Although seismic monitoring remains the primary monitoring method, Sabancaya volcano is also now monitored by geodetic instruments, remote sensing equipment, video cameras, DOAS, multiGAS equipment, and other periodic measurements such as Self-Potential (SP).

Sabancaya volcano (15° 49.3'S, 71° 52.7'W, 5967 m) is a part of the Ampato-Sabancaya volcanic complex (ASVC) in the Central Volcanic Zone of the Andes (Fig. 1). The ASVC is situated to the southeast of Hualca Hualca volcano (HHV) (6025 m asl), a Pleistocene stratovolcano. To the north of the ASVC and HHV lies an extensive system of active faults and lineaments. This includes the Huambo-Cabanaconde Fault Zone to the northwest, which is composed of the Trigal, Solarpampa and Mojopampa normal faults and minor subparallel faults; the Ichupampa normal fault to the northeast; the Pampa Sepina fault zone and Sepina lineament to the east and northeast; and the Huanca normal fault to the southwest (Sébrier et al., 1985; Huaman-Rodrigo et al., 1993; Mering et al., 1996; Machare et al., 2003). Most of these faults were first mapped using field data and satellite images; the existence of several of them was

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Structural and locaion map of Sabancaya volcano as part of the Central Volcanic Zone (CVZ) of the Andes of southern Peru. Red rectangle in the inset maps indicates the study area shown in the enlarged map. Most of the mapped faults are obtained from the Geological and Mining Catastre (Geocatmin) system (https://geocatmin.ingemmet.gob.pe/geocatmin/). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
corroborated with earthquake swarms and interferograms (e.g., Solar-
pampa, Mojopampa, Sepina, Ichupampa), and most of these mapped
faults are available in Geocatmin system (https://geocatmin.ingemmet.
gob.pe/geocatmin/).

The seismicity at Sabancaya from 1990 to 2020 was typically shallow (within 15 km below the summit) and swarm-like in nature, i.e., concentrations of high-frequency earthquakes in time and space, rather than mainshock-aftershock sequences and the magnitudes rarely exceeded M4.5. Much of this seismicity reflects the reactivation of active fault systems related to magma intrusion (Rodríguez and Uribe, 1994; Antayhua et al., 2001, 2002; Jay et al., 2015; MacQueen et al., 2020; Boixart et al., 2020). However, the interaction of the volcanic and tec- tonic zones remains poorly understood. Multi-parametric seismic anal- ysis should help us better understand their complex relationship. In addition, a total of 67 groundwater sources have been identified sur- rounding the ASVC, 40 of which correspond to hot springs whose tem- peratures are above 20 °C (source Geocatmin system), including the Pinchollo geyser north of HHV (Fig. 1). Furthermore, preliminary results of SP measurements around ASVC (Alvarez, 2017; Puma et al., 2018) indicate that both volcanoes have active shallow hydrothermal systems, while the HHV does not, which would indicate that the thermal springs in HHV are of meteoric origin. These observations also indicate the presence of aquifers and hydrothermal systems in the surrounding region.

In this paper, we study the long-term evolution of the 2013–2020 unrest at Sabancaya volcano, including the pre-eruptive and eruptive stages. Although our analysis is largely based on seismic monitoring data, other data and observations greatly contribute to our under- standing of the volcanic activity. Using the seismic data, we estimate an optimal one-dimensional (1D) velocity model, then use it to relocate hypocenters. Next, we analyze earthquake swarms and the time occur- rence of repeaters associated with the two major eruptive events: the eruption onset followed by the crater migration and the dome growing acceleration in 2019. Lastly, we estimate apparent seismic velocity variations (AVV) using ambient noise cross-component and coda wave interferometry. In addition to these seismic analyses, 1) geodetic data were analyzed to interpret deformation patterns associated with the eruptive process, particularly during the dome extrusions in 2019, 2) visual images were processed to characterize the evolution of plumes, 3) satellite thermal data were used to identify thermal anomalies (Coppola et al., 2015), and 4) aerial images (obtained by drones and satellites) allowed us to measure the evolution of the lava dome. Finally, by inte- grating the multi-disciplinary data, we characterize the 2013–2020 eruptive episode and discuss its implications for the magma plumbing system beneath Sabancaya volcano.

## 2. Geological setting

### 2.1. *Geological history*

The Quaternary volcanic arc in southern Peru is composed of at least 16 volcanic centers and several monogenetic fields (e.g., Siebert et al., 2010; Macedo Sánchez et al., 2016), 12 of which are monitored in real time by the IGP (Machacca et al., 2021). The morphology of the ASVC, slightly affected by glacial erosion, suggests a younger age than HHV. HHV, the oldest stratovolcano consisting mainly of lava flows, was probably constructed during the Pleistocene (Gerbe and Thouret, 2004; Thouret et al., 2007). Currently, it shows a horseshoe-shaped caldera morphology, open to the north (Colca canyon) that can be linked mainly to the collapse of the northern flank (Fig. 1) and is associated with intense glacial erosion. Sabancaya and Ampato form a less dissected massif, aligned along a N50°E trend, decreasing in age from SW to NE (Fig. 1). An andesitic lava flow emplaced over Pliocene ignimbrites and forming the basement of the two stratovolcanoes places the growth of Ampato and Sabancaya during the late Pleistocene to Holocene (Gerbe and Thouret, 2004).

The eruptive chronology of the ASVC is detailed by Samaniego et al. (2016) and comprises several distinct stages. The Ampato basal edifice was built during at least two cone-building stages dated at 450–400 ka and 230–200 ka (Samaniego et al., 2016). After a period of quiescence, the Ampato upper edifice developed on the remnants of the Basal edifice about 80 to 70 ka ago. This edifice comprises several cone-building stages that successively shaped the current peaks of Ampato: the northern, southern, and central cones. The southern peak, which is the largest, was built over several effusive episodes between 40 and 20 ka (Samaniego et al., 2016). The Baylillas tephra fallout sequence and the Corinta Plinian fallout deposit testify to the intense explosive activity during this volcanic stage. The last cone-building stage of Ampato constructed the central cone from 20 to 10 ka, and probably continued until Holocene times. During the Holocene, eruptive activity migrated to the NE and built up the mostly effusive Sabancaya edifice.

The edifice of Sabancaya developed in two stages: (1) a basal edifice represented by a sequence of blocky lava flows that spilled out as far as 8 km from the summit toward the W to SE, onto the older rocks of the Ampato and HHV probably as early as the Holocene-Pleistocene boundary (Samaniego et al., 2016; Bromley et al., 2019) and (2) a young cone represented by a sequence of lava flows rests in discordancy on the lavas of the previous stage. This unit also includes a young cone, covered by pyroclastic material, with an active summit crater. These lava flows reach 4–5 km from the vent, have a thickness of 40–60 m and consist of porphyritic andesites and dacites (Samaniego et al., 2016). The northern cone contains an active semi-circular crater with a maximum diameter of 384 m, which migrated to the NW sector during this eruptive episode. Inside this crater lies a recent lava dome emplaced on 5 February 2017.

2.2. *Eruptive history*

The 1986–1998 eruptive activity ended a dormant period of about 200 years (Siebert et al., 2010). The previous recorded volcanic activity occurred during the 18th century, as suggested by Spanish chronicles (Travada and Córdova, 1752; Zamácola and Jaúregui, 1888). From the 18th century to 1986, Sabancaya was characterized by weak fumarolic activity visible in the crater. Fumarolic activity increased beginning in December 1986 and frequent pulses of steam and ash emission were reported in June 1988 (Global Volcanism Program, 1988). The activity then gradually increased until the few days of intense seismic activity preceding the crisis (VEI 2–3) that began on 29 May 1990 (Global Volcanism Program, 1990). Vulcanian explosions generated 1–5 km high ash plumes sometimes reaching heights of 7 km. Ash fall was dispersed >12 km from the summit. In September 1990, 10 to 15 seismic events per day were recorded and this rate increased to up to 50 in October 1990. These events were located about 10 km NE of the crater (Thouret et al., 1994). The eruptive activity at Sabancaya volcano consisted of alternating low and moderate magnitude (VEI 1 and 2) Vulcanian events. The bulk volume of the tephra emitted during the climactic phase of the eruption from May to October 1990 has been estimated at 0.025 km³ (Thouret et al., 1994).

The current eruptive activity ended a 15-year repose period. During 2013–2015, Sabancaya emitted white and blue gases that eventually rose to 2000 m high. Subsequently, on the 6 November 2016, a new eruptive process began. From 2016 until the present day, Sabancaya has been displaying moderate explosive activity characterized by the continuous emission of fine ash and gas rising to 4000 m above the summit of the volcano. Table 1 lists the chronology of the most important events observed at this volcano.

3. *Multi-parametric monitoring and data analysis*

In this study we use data from ground-based monitoring equipment and satellite techniques, observations from previous reports, and available photos spanning from 2013 to 2020. Although seismic methods comprise the bulk of the analysis, we begin by summarizing other
<caption>**Table 1**</caption>

Summary of the main events of eruptive and seismic activity of Sabancaya volcano between 2013 and 2020.

<table><thead><tr><td>Date</td><td>Description</td></tr></thead><tbody><tr><td>22 February 2013</td><td>A new unrest episode started with 4 earthquakes with magnitudes >4.5 in 24 h.</td></tr><tr><td>23 February 2013</td><td>The first temporary seismic station was installed, recording a high rate of seismicity in the area.</td></tr><tr><td>17 July 2013</td><td>M-5.9 earthquake occurred on the northwest flank of HHV, 18 km from the active crater of Sabancaya.</td></tr><tr><td>01 January – 14 August 2016</td><td>New fumarolic fields, located 500–1000 m to the north, northwest, and northeast of the crater, were activated, one of which generated emissions up to 500 m high, visible at 30 km.</td></tr><tr><td>15 August 2016</td><td>A M5.3 earthquake occurred on the Ichupampa fault system, ~30 km from the active crater. This is the most recent large earthquake observed during the unrest period.</td></tr><tr><td>2 November 2016</td><td>One large shallow VT event occurred.</td></tr><tr><td>6 November 2016</td><td>The first explosion with a dense ash column started at 15:40 UTC, initiating the new eruptive activity with highly recurrent Vulcanian explosions that continue to the present (January 2023). Since the eruption onset, the fumarolic fields located outside the crater have decreased in activity, although some of them remain active today.</td></tr><tr><td>24 December 2016</td><td>Apparent velocity decrease began, exceeding the limits of the seasonal fluctuations observed in years prior.</td></tr><tr><td>5 February 2017</td><td>Sentinel-2 images show new thermal anomalies at the northwest sector of the initial crater, the first evidence of the migration (160 m) from the old to the new crater.</td></tr><tr><td>February–July 2017</td><td>Sentinel-2 images show thermal anomalies inside both craters simultaneously, indicating that both craters were active during this period.</td></tr><tr><td>February 2018–April 2019</td><td>Explosive activity decreased with respect to the previous stage and remained constant until May 2020.</td></tr><tr><td>2 May 2019</td><td>Satellite images show the presence of a dome in the new crater, which remains active today.</td></tr><tr><td>12 August 2019</td><td>High-frequency event rates increased (Fig. 11b), culminating in an acceleration of dome growth on 1 October 2019.</td></tr><tr><td>2020</td><td>Activity decreased, similar to February 2018–April 2019, and remains low.</td></tr></tbody></table>

sources of monitoring data.

### 3.1. Deformation measurements

Continuous measurements of ground deformation at Sabancaya volcano started in December 2018. The GPS network consists of two Trimble NetR9 receptors with choke ring antenna, recording at an interval of 30 s, collocated with the SABA and HLCA seismic stations (Fig. 2). The time series of station SABA is presented in the supporting information (Fig. S1).

In addition, an interferogram was obtained by comparing images in the ascending pass of the Sentinel-1 Mission (relative orbit = 47) on 24 May 2015 and 27 May 2020. The result obtained with the interferometric synthetic aperture radar (InSAR) technique is displayed in Fig. S2. Additionally, we computed the temporal evolution of the relative displacement from Sentinel 1 data, on the N flank of the crater (near station HLCA) with respect to the reference point (near station PATA) between 2015 and 2020 (Fig. S1). The results suggest that during 2020 a continuous uplift process remained active north of Sabancaya volcano. This uplift is a continuation of the uplift process described by Boixart et al. (2020) for the period of 2014–2019 and MacQueen et al. (2020) for the period of 1992–2019. The deformation identified to the north of Sabancaya volcano shows at least 7 complete fringes from blue to red, which indicates a positive displacement (~20 cm), with respect to the direction of view of the satellite.

### 3.2. Video camera network

The IGP deployed three video cameras at SB1, CAJA and MRCA stations (Fig. 2) on October 2016 and April and December of 2019, located 31.8 km NE, 11.1 km SE and 9.2 km NW from the active crater, respectively, which recorded one picture every 30 s. Plume elevations from 2013 to 2016 were estimated using pictures taken during field surveys. Since October 2016, ash plume elevations have been obtained using continuous images taken by the monitoring video cameras. The maximum ash plume elevation of each event before the explosive activity (between February 2013 to 6 November 2016) was <2 km. After 6 November 2016, the maximum ash plume elevation increased, reaching >4 km during several explosions. The coloration of the plume got darker after the first explosions of November, indicating higher ash concentration (Fig. S3). The bluish coloration of some plumes is interpreted as the presence of sulfur dioxide gas (SO₂), which reacts with water in the atmosphere to form sulfuric acid (H₂SO₄). When it cools, it forms fine particles of aerosol scattering light, creating a bluish visible plume (e.g., Zelenski et al., 2015).

### 3.3. Remote sensing observations

Sabancaya volcano is also monitored by satellite-based remote-sensing on local and regional scales. Monitored parameters include SO₂ flux, hot spots in the crater, ash dispersion, topographic changes. To track the topographic changes, particularly extrusion of the lava dome in the summit crater, we use the satellite images available at Sentinel-2 Playground™ imagery (https://apps.sentinel-hub.com/sentinel-pl ayground/) that delivers one image every five days in 13 bands. Thermal emissions by hot bodies are detected in the Short-Wave InfraRed (SWIR) channel (bands 12, 8A, 4), with a spatial resolution of 20 m/pixel. We use SWIR images to detect thermal anomalies and to retrieve their temporal and spatial information (e.g., timeseries, location, size, and shape) from the lava dome into the crater.

Thermal anomalies at Sabancaya volcano are monitored by the MIROVA system (Middle InfraRed Observations of Volcanic Activity; http://www.mirovaweb.it/?action=volcanoDetails&volcano_id=354006). The system is based on the analysis of infrared data acquired by the Moderate Resolution Imaging Spectroradiometer sensor (MODIS) and uses the Middle InfraRed Radiation (MIR) recorded with 1 km² resolution to detect, locate and measure the heat radiated from volcanic activity (Coppola et al., 2015). MIROVA measures the volcanic radiative power (VRP) using the MIR-method (Wooster et al., 2003). Results in near real-time are displayed on its website (www.mirovaweb.it).

### 3.4. Seismic monitoring and analysis

#### 3.4.1. Seismic network

A network of seismometers has been developed since 2013. On 23 February 2013, two temporary stations, CAJ0 and SAB0, were deployed at approximately 11 and 3 km away from the active vent, respectively. In March, these two sites were reoccupied (now SABA and CAJA) and a third was added (PATA), all as permanent stations with telemetry. Additionally, between April and December 2013, six temporary stations were deployed (ACHM, OCHK, CBN2, LMUC, TRKN and SLL2) (Fig. 2). In 2013, a total of up to nine seismic stations were operating simultaneously. Between 2014 and 2017, the SABA, CAJA and PATA permanent stations recorded continuously, plus three temporary stations (AMPA, HLCA and OCHK). In addition, data from stations CHVY, HSAL, MISA and UBI1 of the national seismic network of IGP were used to complement the seismic network. Later, in 2018 and early 2019, 3 new telemetered seismic stations were installed (MRCA, PNCH and HLCA). Since March 2019, the IGP has had seven permanent seismic stations dedicated to monitoring the volcano, all of them being telemetered to Arequipa.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** a) Map of the permanent and temporary seismic network (blue and gray inverted triangles, respectively) around Sabancaya volcano. The yellow circles show the position of GPS stations. The green octagons display the positions of video cameras. The red star corresponds to the position of the active vent. The small black circles indicate the principal villages or reference places mentioned in the text, the red squares correspond to the position of geothermal sources. b) Activity period of seismic stations, in black broadband stations and gray short period stations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
In this study, we used continuous data recorded at the closest broadband station SABA (Guralp 40-T; 30 s – 50 Hz) during the period from 2013 to 2020 to create the seismic catalog, estimate seismic energy and compute seismic velocity changes. The other 14 stations surrounding the volcano, including three stations of the national network, were used mostly for source location.

### 3.4.2. Seismic velocity model, earthquake locations, magnitudes, and energy

In order to determine a 1D velocity model for the study area, we used the Velest package (Kissling et al., 1994; Kissling et al., 1995). This method looks for the model that produces the smallest travel time residuals for the locations of a set of events. We only used data from the 2013 campaign, when the station coverage was maximal, and we selected events according to the following criteria: (a) the events must be recorded by 5 or more stations, (b) the azimuthal gap should be <180°, and (c) the root mean square (RMS) residual should be <0.35 s. A total of 683 seismic events were selected, yielding 3784 P phases and 3594 S phases. The output velocity models (Fig. S4) were obtained by solving the coupled hypocenter-velocity problem iteratively (7 times) with Velest. The final velocity model (Table S1) was taken as the average of the 50 best models, those that exhibited the lowest misfit (RMS residual) values (< 0.13 s). In this model, for depths below sea level, velocities are consistent with the IASP91 models (Kennett and Engdahl, 1991) and with the regional S-wave velocity model obtained by ambient noise tomography studies (Ma and Clayton, 2014). The improvements on the precision of source location (RMS and horizontal and vertical error reduction) with the new model are presented in Fig. S4.

The seismic sources of the complete dataset were then located using the new velocity model and the NonLinLoc algorithm (Lomax et al., 2000). This method is based on a probabilistic approach (Tarantola and Valette, 1982) and provides a more complete description of uncertainties. Earthquakes are filtered by number of stations ≥ 4, gap ≤ 270°, horizontal error ≤ 5 km, vertical error ≤ 10 km and RMS < 0.35 s. The results are presented in Section 4.1.

To estimate the magnitudes of earthquakes at Sabancaya volcano, we developed an appropriate local magnitude scale for this region. Local magnitude scale (M_L) was defined by Richter (1935) as follows:

$$M_L = \log A - \log A_0 + S \quad (1)$$

where A is the amplitude measured in millimeters on a Wood-Anderson seismograph, S is the station correction factor, and – log A₀ is the distance correction factor. For this study, we estimate – log A₀ for Sabancaya volcano based on 20,576 amplitude measurements of 2412 seismic events recorded at 6 seismic stations and we used 17 km as the reference distance (Hutton and Boore, 1987). Results of the inversion procedure following Miao and Langston (2007) is summarized as:

$$-\log A_0 = (1.7707 \pm 0.029)\log(r_h/17) + (-0.0064 \pm 0.001)\log(r_h - 17) + 2 \quad (2)$$

where r_h is the hypocentral distance in km. The results differ from previous regional studies (Condori et al., 2017) but are similar to results in other volcanic environments (Fig. S5) with 17 km as the reference distance (e.g., Greenfield et al., 2018).

Finally, seismic energy has been calculated using the eq. (3) formulated by Johnson and Aster (2005) on the vertical (Z) component of station SABA,

$$E = \frac{2\pi r^2 \rho c S^2}{A} \int_{0}^{r} y^2(t) dt \quad (3)$$

where r is the distance from the source (3000 m), ρ is the density (2600 kg m⁻³), c is the P wave velocity (3000 ms⁻¹), A is the attenuation correction, S is the seismic site response correction (fixed at 1), y(t) is the ground particle velocity, and T is the duration of the analyzed time

window. We assumed that the source was at a shallow depth below the crater.

### 3.4.3. Seismic event classification

To classify seismic events, we used the records at station SABA, which is the closest station to the crater and thus can detect the smallest events. We adopted a terminology similar to previous studies (Lahr et al., 1994; McNutt, 1996; Chouet and Matoza, 2013; White and McCausland, 2016). Event types are based on signal characteristics, such as frequency content, type of onset, duration of the signal, and for some types, the location of the source with respect to the crater. The main types of signals observed at Sabancaya are volcano tectonic (VT), long period (LP), hybrid (HYB), tornillo (TOR), very long period (VLP), gas emission, explosion (EXP), spasmodic tremor (TRE), harmonic tremor (TRA) and lahars (LAH) (Fig. 3). VT events are characterized by clear onsets and high frequency content (up to 25 Hz) and are further grouped based on their distance to the eruptive vent. Classical VT events occur below the eruptive vent in the summit crater (usually with T_S-T_P < 1.5 s at SABA station). VTs that occur farther from the vent, generally from ~3 to 30 km or more in horizontal distance and between 2 and 15 km in depth, are termed distal VTs (dVTs) following White and McCausland (2016). At Sabancaya, this distinction by distance is used for descriptive purposes only, as some earthquakes classified as dVT may actually be proximal to the magma chamber, which is offset to the N of the crater, as discussed later. HYB are catalogued based on their waveform and frequency content (clear high frequency onset [5–15 Hz] but predominantly LP in the coda). Gas emissions are cigar-shaped emergent signals with frequency contents in the range of 4–15 Hz (Fig. 3). These events were observed mainly in the pre-eruptive period and are visually correlated with gas emission pulses. TOR events are a type of LP event mostly related to active hydrothermal systems characterized by slowly decaying coda waves (e.g., Torres et al., 1996). They are long lasting (up to several minutes), with emergent or impulsive onsets, and typically have 1–3 narrow spectral peaks (< 5 Hz.) (Fig. 3). VLP is a characteristic long-period event [3.0–1.2 s] that accompanies, in most of the cases, explosion and steam emission signals at Sabancaya volcano. TRE are sustained signals that can last from tens of seconds to hours and occur over a frequency range from 0.8 to 10 Hz, with a wide range of spectral content. TRA episodes are continuous oscillations characterized by sharp spectral peaks with multiple harmonic overtones which glide as a function of time (e.g., Benoit and McNutt, 1997; Lesage et al., 2006). The fundamental frequency ranges between ~0.9–3 Hz with typical durations of 2–7 min but occasionally longer. TOR events occur only in the pre-eruptive period, while TRA episodes only occur during the eruption. We use the terms ‘low frequency events’ (LFEs) to include LP, HYB, TOR, VLP and EXP, and ‘high frequency events’ (HFEs) to refer to dVT and VT events. EXP and LP events are distinguished by their amplitude levels and the characteristic signature of EXP, which consists of two phases (coupled events) with a time separation of 10–35 s at 3050 m to the crater.

### 3.4.4. Repeating earthquake analysis

Earthquakes with similar waveforms are often referred to as multiplets, families or repeaters, and have been identified in many tectonic and volcanic settings. For example, repeaters have been observed along the San Andreas and Calaveras Faults in California and are interpreted as repeated rupture of a fault patch (e.g., Vidale et al., 1994; Nadeau et al., 1995; Ellsworth, 1995). At volcanoes, mainly LP repeaters are commonly observed during dome building eruptive sequences (e.g., Frémont and Malone, 1987; Rowe et al., 2004; Power and Lalla, 2010; Thelen et al., 2011), during basaltic eruptions (e.g., Battaglia et al., 2003; Saccorotti et al., 2007), and during non-eruptive time periods (Lin, 2017) and have been associated with different conduit processes in various ways. However, repeaters of VT events can also occur throughout the eruptive cycle on a volcanic edifice, including on nearby tectonic faults in complex tectonic environments (e.g., White et al.,
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Waveform and spectrum for the different types of volcanic earthquakes observed at Sabancaya volcano during 2013–2020. Signals are recorded at SABA on the vertical component. Acronyms in the figure are: dVT (distal volcano tectonic), VT (volcano tectonic), LP (long period), HYB (hybrid), TOR (tornillo), VLP (very long period), EXP (explosion), TRA (harmonic tremor), TRE (spasmodic tremor) and LAH (lahar).</figcaption></figure>

2011; Massin et al., 2013; Wellik et al., 2021).

Repeaters have become increasingly important in volcano seismology and several approaches have been developed to analyze them (e. g., Matoza et al., 2015; Lengliné et al., 2016; Hotovec-Ellis and Jeffries, 2016; Shapiro et al., 2017). In this study we used the Repeating

Earthquake Detector in Python (REDPy; Hotovec-Ellis and Jeffries, 2016) program for detecting repeating earthquakes. This program detects events by applying the classic short-time average to long-time average (STA/LTA) algorithm, and cross-correlating core events (representative waveforms of each family obtained by the Reachability
functions (Ankerst et al., 1999)) with other detected events. If an event correlates above a defined threshold, it is defined as part of a family, otherwise it is labeled as an orphan event unless a match is later found. An advantage of this program is that it does not require predefined templates or training data to be provided a priori.

The families of events are computed using records from the closest station to the crater (SABA, Z component) in the band 0.8–10 Hz, and are classified based on their core event, as described in Section 3.4.3. We used LTA and STA windows of 8 and 0.8 s, respectively, and a STA/LTA ratio to turn a trigger on (off) was set to 3.0 s (2.0 s). We tested the method with different correlation coefficient (CC) values and time windows and chose CC = 0.85 and a 40 s window as our preferred values, which we then used to detect event families for the 2016–2020 period. For this period, the percentage of repeaters from all triggered events varies from a few to 5%, but there are two main periods in which the percentage of repeaters reaches 10% (sometimes even exceeding 20%) of the total triggered events (Fig. 4).

Here we focus on repeaters for the two periods that exhibit major morphological changes in the summit crater; November 2016 to September 2017 (associated with crater migration) and January to December 2019 (associated to accelerated dome growth). For July 2016–February 2017, a total of 2538 families and 16,167 repeaters were identified, while for year 2019, 2202 families and 14,394 repeaters were identified. The waveform, frequency content, and signal duration of some families are presented in Fig. 5 and discussed later in Sections 4.4 and 4.5.

### 3.4.5. Estimation of seismic velocity changes by ambient noise cross-correlation

In recent years, seismic velocity changes associated with volcanic activity have been observed at an increasing number of volcanoes (e.g., Brenguier et al., 2008; Budi-Santoso and Lesage, 2016; Bennington et al., 2018; Yates et al., 2019; Machacca-Puma et al., 2019). These tiny variations are generally estimated by coda wave interferometry using cross-correlation functions of seismic noise recorded by two stations, but can also be derived using two components of a single station (e.g., De Plaen et al., 2016; Bennington et al., 2018; Caudron et al., 2022). In this study, we used cross-component correlation functions for components of station SABA.

Data pre-processing consists of several steps. We first prepared waveforms in one-hour long segments and removed the mean and trend. As a second step, we down-sampled the signal to 50 Hz, applied spectral whitening and filtered the records in different frequency ranges (0.1–1, 0.3–1, 1–3, 3–5 and 5–8 Hz). To suppress high amplitude events, we

performed amplitude normalization in the time domain by amplitude clipping (Bensen et al., 2007) at 3 times the RMS. After that we computed the noise cross-component functions (NCFs) and stacked them over 24 h. Daily NCFs calculated with <6 h of data are discarded because of their low signal-to-noise ratio. Fig. S6a displays an example of correlogram calculated over 2019.

The velocity variations (dv/v) were estimated using the stretching method (Lobkis and Weaver, 2003; Sens-Schönfelder and Wegler, 2006), which consists of stretching or compressing the reference or the current NCF in order to maximize the correlation coefficient (CC) between the two functions in selected delay windows in the coda. In the general case, the perturbation is not uniform, and dv/v should be considered as an ‘apparent’ velocity variation (AVV). Fig. S6b displays an example of AVV resulting from stretching the causal and acausal part of NCF in the coda, with a resolution of one value per day.

Finally, in order to improve the precision, to reduce the volcanic source effect and to increase the robustness of our estimation of AVV, we calculated the velocity variation without reference following Brenguier et al. (2014) and Gómez-García et al. (2018) (see also supplementary material). In Fig. S7 we present the results of AVV in five frequency bands for SABA station. High frequency bands are affected by large VTs and tectonic earthquakes to a greater extent compared to low frequency ones (Fig. S7). For these reasons we decided to compute long term (2014–2020) AVV only in the frequency range of 0.1–1.0 Hz. For computational reasons, we split out the analysis of dv/v into segments of 14 months, with overlap of one month. Results of AVV are presented in Section 4.2.

## 4. Volcanic unrest and eruption

### 4.1. Spatio-temporal distribution of seismic activity

During the study period, Sabancaya volcano was characterized by high levels of seismic activity, located on several surrounding faults, which sometimes were activated simultaneously. In this section, we present a brief description of the spatio-temporal distribution of seismic activity based on the hypocenter locations carried out with the local seismic network and the new velocity model presented in Section 3.4.2.

The first strong seismic unrest near Sabancaya volcano was noticed when four earthquakes with magnitude 4.6, 5.2, 5.0 and 4.5 were reported within a couple of hours by the national seismological service of Peru on 22 February 2013 (Tavera et al., 2013). These earthquakes (white circles in Fig. 6a) were felt in the Colca canyon towns and were located by the national seismic network between 7 and 13 km from the

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Daily number of triggers of all types produced by REDPy (black line), and percentage of events belonging to repeaters (red line) at SABA station, during 2016–2020. The percentage of repeaters clearly increased before the onset of eruption in 2016 and the episode of dome extrusion in 2019. The top horizontal hatched bars correspond to the eruption period, and gray shade zone correspond to period of dome acceleration. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Waveforms (gray), corresponding core events (black), and spectrum for 20 repeaters recorded on the vertical component of SABA in 2019. The waveforms are band-pass filtered with 4 poles, zero phase filter and corner frequencies of 0.8 and 10 Hz and all amplitudes are normalized to the maximum. The type of volcanic earthquakes, family ID and the numbers of events in each family are indicated.</figcaption></figure>

active vent on the southeast side of Sabancaya. Hypocenter uncertainties for these first four M4+ earthquakes are on the order of 10 km. At this time there were no local seismic stations, thus they were located using stations >50 km. A day after these events, temporary stations started to be deployed around Sabancaya, but seismicity remained low without significant (>20 localized events/day) bursts. A new swarm peaked on 17 July 2013, when an earthquake with a moment magnitude (Mw) of 5.9 occurred 15 km NW from the crater, along the Mojopampa fault. The focal mechanism for this event suggests normal faulting, consistent with the trend of the Mojopampa fault and larger Huambo-Cabanaconde fault system.

In 2014, few earthquakes were located due to limited data avail-
ability (Fig. 2b). The seismic activity was located about 8 km NNE and E
of the active crater, near a place known as Hornillos, and was about 4–8
km deep (Fig. 6b). Some events occurred at distances <3 km from the
crater to the NE and with depths <0 km a.s.l.

During 2015, the seismicity was located mainly to NNE of Sabancaya volcano in the N32°E directional fault system that runs through Maca village, an area of high geological activity and mass movement (Lacroix et al., 2015). A second important source was located on the NW side of the volcano, between the Solarpampa and Mojopampa faults, similar to activity observed in 2013, but 7.5 km away from the M5.9 in 2013. Other minor sources were concentrated to the NE of the crater along the Sepina fault system and Hornillos (NNE) at depths <10 km (Fig. 6c).

In 2016 and 2017, the seismicity was concentrated along the Sepina
fault system with main foci located at 7, 13, and 23 km from the crater
(Fig. 6d, e). Another important seismogenic zone was the intersection of
Sepina fault with the Ichupampa fault system at depths <15 km. On 15
August 2016, a M5.3 strike-slip earthquake occurred at this intersection
(Fig. 6d). A new lineament system (not mapped) was identified in 2017
around Coporaque village, at 29 km NE of the crater, parallel to the
Ichupampa fault system (Fig. 6e).
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Yearly maps of seismicity around Sabancaya volcano during 2013–2020 (from a to h). The colors represent the depth of the earthquakes. Red triangles indicate the position of Sabancaya summit, white inverted triangles indicate seismic stations, star in d) indicates the pre-eruptive shallow VT location, and black lines indicate main fault systems. Focal mechanisms were obtained from the Global Centroid Moment Tensor catalog (https://www.globalcmt.org/CMTsearch.html). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
In 2018, the seismicity was concentrated in multiple foci, most of them located to the north and northeast of the volcano at depths shallower than 5 km below the surface of the volcano (Fig. 6f).

The seismicity in 2019 was located on the NW, N, NE, and E flanks of the volcano (Fig. 6g). During the second semester, strong seismic activity occurred on the east side of the volcano and preceded the increase of volcanic activity and dome growth. During this year, eight events were recorded with magnitudes greater than $M_L$4. The largest earthquake ($M_L$4.9) detected on 21 May 2019 occurred 5 km north of the crater, above the deep inflationary source, with a focal mechanism corresponding to an oblique normal fault dipping 56° to the northwest and parallel to the Sepina fault system.

In 2020, the seismic activity remained at high levels in the N, NE, and E flanks, similar to that of 2019 (Fig. 6h). The largest earthquake ($M_L$4.5) occurred 11 km to the north of the crater, inside the collapse scarp of HHV north flank. Its focal mechanism indicates an oblique normal fault dipping 56° to the north, in one of the most active fault systems.

The focal mechanisms of the large earthquakes that occurred on the northwest flank (close to the inflation source) consist of dip slip movements, sometimes with a small strike-slip component (Fig. 6a, c). The

earthquakes located NE and SE of the volcano have a strike-slip mechanism (e.g., 15 Aug. 2016) with a dip-slip component for some of them (Feb. 2013). There is a clear difference between the focal mechanisms for earthquakes located NW of the volcano and for earthquakes located at SE.

## 4.2. Apparent velocity variation observed during the eruptive stage

Fig. 7 illustrates the results of the AVV analysis for station SABA (Fig. 7a), nearby environmental data (Fig. 7b), and for several pairs of stations for the 2019 time period (Fig. 7c, d). The AVV calculated in the lowest frequency range of 0.1–1 Hz, using the vertical and east (ZE) components, present a decrease of less than –0.4% from 2014 to 2020, related to volcanic activity, especially at the beginning of the eruption and some months before the enhanced dome extrusion. This decrease exceeds the fluctuations of ±0.2% observed before the eruptive period and no clear relationship with seasonal meteorological observations can be established. No such long term AVV is detected using the vertical and north (ZN) components (Fig. 7a). Because the SABA station is located on the southern flank of Sabancaya, the north component is approximately

<figure><img src="image_10.png" /><figcaption>Fig. 7. a) Apparent velocity variations calculated for two pairs of components of SABA station. b) Temperature and precipitation at the meteorological station of Chivay (31 km NE from crater, Fig. 2) operated by The National Meteorological and Hydrological Service (Spanish acronym: SENAMHI) of Peru. c) Sensitivity kernel for velocity variations obtained from cross-component correlation functions at single stations, case SABA station, surface wave velocity c = 1.38 km/s, mean free path $l = 5$ km, T = 24 s in the coda. d) Average AVV obtained for single station (Z-E and Z-N) and AVV for pair of stations (see legend).</figcaption></figure>
radial with respect to the crater while the east component is tangential. The observation of velocity variations much larger on the tangential component than on the radial one was also made at Ubinas volcano, Peru (Machacca-Puma et al., 2019). This anisotropy is still not properly explained. We speculate that the location of a noise source below the crater, the cone shape of the andesitic volcanoes, and/or the orientation of cracks, play important roles in this anisotropy. Hotovec-Ellis et al. (2022) suggested in the same way that the geometry of ring faults must be considered to interpret AVV for data from Kilauea Volcano. The largest velocity decreases calculated with the ZE components occurred in December 2016 and, a lesser extent, in February–March 2019. The former period follows by one month the eruption onset and precedes the formation of the new crater. The latter decrease occurred a few months before the rise of extrusion rate. Both episodes might result from rock damaging, crack opening and/or pore pressure increase in the surrounding structure (Mordret et al., 2010; Carrier et al., 2015; Donaldson et al., 2019; Got et al., 2019).

### 4.3. Modifications of the crater morphology

The multiparametric monitoring system illustrates the temporal evolution of unrest parameters (Fig. 8) and reveals significant differences between the pre-eruptive and the eruptive periods. The pre-ruptive period began in February 2013 and was mainly characterized by intense seismic activity dominated by VT and dVT types (Fig. 8b, c) and by emissions of steam and volcanic gases. These emissions were generated in the crater left by the 1990s eruptive sequence and produced plumes of <2000 m above the crater level (Fig. 8g). A first thermal anomaly on Sabancaya volcano was observed in July 2014 during a period of intense phreatic activity, and occurrence of more TOR events.

The eruptive period began on 6 November 2016 and was ongoing as of this writing (January 2023). It is characterized by the occurrence of (almost daily) Vulcanian explosions, which emit volcanic ash up to 4000 m above the crater level. The level of eruptive activity varies strongly, with daily seismic energy of the explosions in the 0.01 to 90 MJ range (Fig. 8a). The AVV decreased by 0.4% in the eruptive period, especially during 2017 and 2019 (Fig. 8d). The thermal anomalies, measured by Volcanic Radiative Power, increased by almost two orders of magnitude between the pre-eruptive and the eruptive periods, reaching 100 MW for some measurements (Fig. 8f). These observations indicate the persistent emplacement of a lava body or dome inside the crater.

Fig. 9 displays a series of satellite images that illustrate the morphological changes in the summit crater of Sabancaya. The main change observed is the formation of a new crater in 2016–2017 and the migration of the active vent about 160 m to the NNW. Between January and July 2017, the two craters were simultaneously active. The time occurrences of activity in both craters, according to satellite images, are displayed in Fig. 8e. Then, beginning at the end of June 2017, the activity was concentrated in the new crater with a main vent located near its eastern edge. Some minor vents inside this crater were also observed on satellite images and could generate two or even three plumes during a single explosion sequence. On 1 October 2019, rapid lava dome growth was observed in the new crater. The extrusion rate of this dome markedly increased at the end of this year after a sequence of strong seismic activity.

HFEs seismicity was predominantly pre-eruptive and decreased in number and depth at the beginning of the eruption, while the LP activity was increasing. From the total number of HFEs (Fig. 8b), approximately 24% correspond to crater VT events, and 25% were located (Fig. 6). The dVT activity occurred in bursts, some of which reached >1000 events per day. These bursts generally correspond to seismic swarms but in some instances, they are sequences of mainshock earthquakes with magnitude >4.5 followed by aftershocks (e.g., 21 May 2019). Figs. S8 and S9 show examples of both types of seismic sequence.

In the following sections, we focus on the two periods of marked

eruptive changes: the eruption onset in 2016–2017 and the accelerated dome extrusion in 2019.

### 4.4. Onset of eruption

The onset of eruption on 6 November 2016 was preceded and followed by several types of phenomena. For this reason, we analyzed in more detail the four-month period before and after of the eruption onset. Fig. 10 shows the enlarged time series of the monitored parameters between July 2016 and February 2017.

The seismic activity in the areas surrounding the volcano suddenly increased on 15 August 2016 when several earthquakes, including a ML 5.3 earthquake, occurred on the Ichupampa fault 25 km northeast of the crater. VT activity did not decay after this event, but rather remained high and swarm-like in the area for a few months and the LP seismicity gradually increased until the eruption onset (Fig. 10b). An ML 1.6 shallow VT event occurred on 2 November at 14:28 UTC (Fig. S10), 4 days before the eruption onset.

During the first months of the eruptive period, the Vulcanian explosive activity progressively heightened and migrated from the old to the new crater, although the two craters remained active together for several months (Fig. 10d). Since December 2016, the apparent velocity variations estimated at station SABA displayed a decrease of up to 0.2% that exceeded the amplitude of the fluctuations observed during the non-eruptive period in January 2017 (Fig. 10e).

Fig. 10c shows the occurrence in time of the hundred most populated repeaters identified in this period. Each family contains 12 to 717 events which are either VT, dVT or LP. HYB events are not reported here because few of these shallow and unenergetic events are recorded at SABA station, 3 km from the crater. Similarly, no significant EXP families occurred during this period. The repeaters are mostly composed of events that are too small to be located. The duration of the repeaters is highly variable, from a few hours to several months, especially those of VT type. Most LP repeaters are long-lived, probably due to their source mechanisms. Among the ten LP families, seven started their activity before the eruption onset and remained active afterward, coincident with the activity in the old crater. The rate of appearance of new dVT repeaters sharply increased in August 2016 on two separate faults and persisted at a high level until the eruption beginning (see Fig. 10c and Section 5.4).

### 4.5. Acceleration of dome emplacement

On 5 February 2017, Sentinel 2 imagery showed for the first time a lava dome inside the new crater, without significant changes until October 2019. The evolution of this lava dome was documented by using series of satellite images and photos taken from drone overflights. Fig. 11d displays the evolution of the base area of the dome through 2019 estimated from these images. Rapid dome growth occurred from September to November 2019 when its surface area increased by a factor of 3.1 (from $18 \times 10^3$ to $57 \times 10^3$ m$^2$). This dome emplacement coincided with a marked increase in the number and energy of explosions (Fig. 11a). This episode of accelerated lava extrusion was preceded by several bursts of VT activity which were superimposed on to the almost constant background LP seismicity (Fig. 11b). A ML4.9 earthquake occurred 5 km NE of the crater on 21 May 2019, accompanied by an increase in VT event rate over a couple of weeks. Another strong swarm, which included several earthquakes with ML > 4, began on 12 August 2019–12 km E of the crater and continued until the beginning of October when the dome growth was accelerating. From February to September 2019, the apparent velocity showed a decreasing trend of about 0.2%. It returned to its initial value during the next two months (Fig. 11d).

Among the 100 most populated repeaters detected in 2019, 84 were composed of VT + dVT events, 11 of LP type and 5 of EXPs (Fig. 11c). Most LP and EXP repeaters were long-lived and were already active before the episode of accelerated dome extrusion. During this episode,
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Time series of the multiparametric observations for 2013–2020 on Sabancaya volcano. a) Daily energy of explosions. b) Daily count of VT + dVT and LP seismic events recorded at SABA station. c) Focal depth bins below sea level of the located earthquakes, large M4.5 earthquakes are plotted in white circle and pre-eruptive shallow VT in white star. d) Apparent seismic velocity variations, the shaded gray zone corresponds to the base fluctuation of AVV observed during 2014–2015. e) Periods of activity at the old and the new craters observed by Sentinel Playground imagery, the vertical lines correspond to the dates with clear visibility. f) Thermal anomaly measured by Volcanic Radiative Power. g) Plume elevation. In all panels, the horizontal hatched bars correspond to the erup- tion period.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Images showing the crater migration of Sabancaya volcano. a) Image taken on 6 June 2012; the blue line corresponds to the limit of active crater on that date. b) Image taken on 26 May 2016; the bluish zones correspond to new fumaroles fields activated during that year. c) Image taken on 2 May 2019, the blue line corresponds to the old crater, the black line corresponds to the new crater and the white line corresponds to the edge of the lava dome. d) Image of the new crater on 7 November 2019, showing the extension of the dome surface since 2 May (from dashed to solid white line). e) Image of the new crater in 2021, showing dome destruction and several vents placed where the 2019 dome was located (dashed white line). Orange hatched area corresponds to the most active vent observed during the eruptive stage. Sources: Google Earth for a,b,c and DigitalGlobe for d (courtesy of Volcano Disaster Assistance Program - VDAP). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

they continued their activity and two new repeaters, one LP and one EXP, appeared. The VT repeaters have durations much shorter than those of LP and EXP. They appeared generally during seismic bursts distal to the volcano. The most striking example is the appearance of 61 short-lasting VT and dVT repeaters during the ~45 days preceding the dome emplacement. The rate of appearance of new repeaters sharply decreased in October and November. The VT and dVT repeaters are located mainly at a few foci, 5 to 15 km from the crater and at depths of ~8 km (Fig. 12).

# 5. Discussion

## 5.1. Summary of results

The Ampato-Sabancaya volcanic complex and its surrounding region is currently associated with intense seismic activity (of volcanic and non-volcanic origin), large-scale deformation, hydrothermal activity, explosive and effusive eruptions, and gas emissions. The most recent eruptive episode started on 6 November 2016 and was preceded by a period of 3 years and 8 months of unrest characterized by the reactivation of many tectonic faults and enormous number of dVT events, the inflation of a deep reservoir and increased degassing from the crater. During the pre-eruptive and eruptive phases, hundreds of earthquakes were recorded per day including a total of 15 events with ML larger than 4.5; dVT seismicity dominated before the eruption onset with some bursts of activity including >1000 events. No significant increase of VT seismicity below the summit crater was observed in the few days

preceding the eruption onset. LP seismicity became generally dominant during the eruptive period. Explosive activity was preceded by a shallow VT event ($M_L$ = 1.6) located ~1.7 km below the volcanic edifice (Fig. S10). This event may be due to the breaking of the brittle-ductile transition boundary cap (Fournier, 2007), allowing magma to reach the crater surface. The sources of dVT events were located mainly from the NW to the East of Sabancaya at distances of up to 25 km. Very few earthquakes were located at the western and southern sides of the volcano. A significant proportion of dVT and VT events belong to families of similar waveforms, especially during some bursts of activity that preceded by a few weeks the eruption onset and the acceleration of the dome extrusion in 2019 (Fig. 4). The observations that are presented, analyzed and interpreted for this 8-year-long period are consistent with previous works on this volcanic system (Jay et al., 2015; MacQueen et al., 2020; Boixart et al., 2020; Coppola et al., 2022).

In the first sections of the paper, we presented a detailed analysis of the observations. Many questions arise from these results about the relationships between the eruptive activity, the distal VT activity, the magmatic intrusion, the tectonic faults, and the hydrothermal system of the area. In order to integrate this complex set of observations and results in a common framework, we propose in the following sections an interpretative scheme and discuss the different questions and elements of the scenario.

## 5.2. Deformation and volume of intruded magma

The deformation observed by InSAR from 2013 to 2019 extended
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Time series showing the evolution of observations between April 2016 and September 2017. a) Daily energy of explosions. b) Daily count of seismic events recorded at SABA station. c) Occurrence of repeaters observed at the vertical component of SABA station. d) Periods of activity of the old and the new crater observed in Sentinel Playground imagery (horizontal bar), the vertical line markers indicate the date with visible image. e) Apparent seismic velocity variation. f) Satellite thermal anomaly measured by VRP. In all cases the gray vertical shaded zone corresponds to the period when both craters were active, and the horizontal hatched bars correspond to the eruption period.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Time series showing the observations for 2019. a) Daily energy of explosions. b) Daily count of seismic events recorded at SABA station. c) Time occurrence of repeaters. d) Velocity variation at SABA station (gray line) and dome area (black line). In all cases the gray shaded zone corresponds to the period of dome building.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. a) Average source location of HFEs events that form part of the 2016 repeaters. b) Same for 2019 repeaters. The colored circles correspond to the HFE repeaters preceding the eruption onset (ID: 17–84) and to the dome growing (ID: 33–93), respectively. The gray circles correspond to the rest of repeaters. Dashed colored contours indicate the area of inflation between May 2015 and May 2020 (See legend).</figcaption></figure>

over an approximately circular area centered on HHV with a radius of ~20 km (Boixart et al., 2020; MacQueen et al., 2020; Fig. S2). It produced a maximum line-of-sight uplift rate of 3 to 5 cm yr⁻¹. The best models place a spherical inflating source below HHV at 5 to 7 km NNW from Sabancaya volcano at depth of 12 to 15 km below the surface (Boixart et al., 2020; MacQueen et al., 2020). During the period 2013–2019, the volume change rate was estimated as 0.026 to 0.046 km³ per year by Boixart et al. (2020) and as 0.03 to 0.05 km³ per year by MacQueen et al. (2020). These rates are approximately constant over the whole study period and no discrete injection event can be identified in the deformation time series displayed in Fig. S1. Extrapolating these values from 2013 to 2020, we obtain a total volume change ranging from 0.21 to 0.40 km³.

This geodetic estimate of intruded volume allows us to test proxy methods for estimating intruded volume based on seismicity. White and McCausland (2016) proposed an empirical relationship between intruded magma volumes and cumulative seismic moment. In the case of Sabancaya, the cumulative seismic moment depends mainly on the few earthquakes with ML > 4 that occurred in the region, especially the Mw 5.9 event on 17 July 2013. The seismic moments were estimated from the local magnitude using a relation proposed by Hanks and Kanamori (1978). The cumulative value obtained is 1.90 × 10¹⁸ Newton-meters.

Using the equation proposed by White and McCausland (2016), an intruded volume of 0.57 km³ is obtained. Recently, Meyer et al. (2021) reevaluated this relationship using a Bayesian method and more recent data. With their updated equation we obtained an estimation of intruded volume of 0.195 km³. Considering the large uncertainties associated with these proxy estimations, which can be up to 3 orders of magnitude (Meyer et al., 2021), they are generally consistent with the range of values calculated from the deformation models. All the estimations indicate an intruded magma volume between 0.2 and 0.6 km³ for the period 2013–2020.

### 5.3. Distal volcano-tectonic seismic activity

The high level of seismic activity observed in the surroundings of Sabancaya – about 227,000 HFEs detected during this 8-year-long study period, including 15 earthquakes with magnitudes equal to or larger than 4.5 – is not common in volcanic systems. Most of the events occurred between the Huambo-Cabanaconde and the Mojopampa-Hornillos faults, which are two parallel major NW-SE dipping normal

faults mapped and reported in previous studies (Jay et al., 2015; MacQueen et al., 2020). They are especially concentrated at the intersection with perpendicular secondary faults such as the intersection of the Sepina and the prolongation of Huambo-Cabanaconde faults systems at 12 km N56°E from the active crater and that of the Sepina and the Ichupampa fault at ~30 km NE from the volcano summit. The focal mechanisms of the large earthquakes, located close to the magmatic inflation source NW of the volcano, consist of dip slip movement, while earthquakes located NE and SE of the volcano have mainly strike-slip mechanisms, consistent with broad extension due to the deep volcanic inflation source and the regional tectonic framework, respectively.

The regional tectonic movements probably played a role in the accumulation of stress on the faults and in bringing them close to the rupture. However, tectonics alone cannot explain the marked increase of the seismicity rate observed in 2013–2020 with respect to the preceding decade, and a similar increase in M4.5+ seismicity was also noted prior to the previous eruptive period from 1986 to 1998 (Antayhua et al., 2001, 2002). Thus, we postulate that magma intrusion in the deep reservoir plays an important role in producing pressurization of the plumbing system, and, due to the effect of free surface, extensional and shear stresses are produced above the source (McTigue, 1987; Hautmann et al., 2014; Got et al., 2019; Hotovec-Ellis et al., 2022). In this shallow region, the Coulomb stress is progressively increased, which favors rupture on normal faults (MacQueen et al., 2020). In the case of Sabancaya volcano, the depth of the pressure source is about 13 km, which approximately corresponds to the radius of the largest zone in extension at the surface. This is similar to the width of the region where uplift was detected by INSAR (MacQueen et al., 2020; Boixart et al., 2020). Therefore, this mechanism alone cannot explain the triggering of VT activity observed at distances larger than 13 km. Deep magma intrusion can also produce heating and pressurization of the surrounding aquifers. The gas generated by the degassing deep magma reservoir rises through the old permeable conduit and the surrounding rocks and produces heating and pressurization of the shallow hydrothermal system. Fluid pressurization generates a progressive increase in pore pressure that propagates at large distance thanks to the high permeability of the faults and fractures (Talwani et al., 2007) and contributes to the triggering of VT and dVT seismicity including some large earthquakes (White and McCausland, 2016). The latter explanation is consistent with the distribution of earthquakes and the existence of aquifers and hydrothermal systems in the surrounding region.
Additional evidence for a hydrothermal system comes from Puma et al. (2018), who carried out self-potential surveys in the area and identified a large positive anomaly, indicating the upward flow of hot fluids, on the northern and eastern flanks of Sabancaya volcano. They also detected a negative anomaly, associated with the downward flow of meteoric water, on the northeastern side of HHV which is interpreted as an aquifer within an old collapse caldera. The cluster of seismicity located within the corresponding scarp ring fault is probably related to this aquifer. The numerous geothermal springs north of HHV and along the Colca valley (Fig. 1), as well as the Pinchollo and Fuye geysers, are additional evidence of the presence of aquifers.

Conversely, the low level of seismicity south of Sabancaya could be associated to the lack of hydrothermal system in this zone, existence of fewer faults, or to the hydraulic disconnection with the reservoir below HHV. This decoupling could be due to some structural barrier such as a ductile and impermeable zone under Sabancaya volcano.

## 5.4. Seismic activity associated with the main eruptive episodes

During the eruptive period, two episodes have produced significant changes in the crater morphology and in the features of the seismic activity: 1) the beginning of the eruption in November 2016 which was shortly followed by the migration of the active crater and 2) the marked increase of the extrusion rate of lava dome in 2019. Both eruptive events were preceded by reactivation of VT seismicity and especially by a strong burst of activity that occurred 2.5 and 1.5 months before each, respectively. During these two bursts, the rate of appearance of new VT and dVT repeaters sharply increased (Figs. 4, 10, 11 and 12) and decreased almost immediately after the eruption onset and the dome growth acceleration.

Some of the new repeaters remained active during periods of a few days to a few months in 2016 and up to a few weeks in 2019. In contrast, since the beginning of the eruptive period, the dominant seismic activity is composed of LP events, some of which belong to seismic repeaters. LP and EXP repeaters (Fig. S11) are much more persistent, even at the onset of the eruption or during the episode of dome extrusion acceleration. This observation suggests that the sources of these repeating events and dome extrusion probably were located at different positions, and that explosive events were not much disrupted by the changes in eruptive activity, as suggested by the presence of the active vent on the NE edge of crater (Fig. 9c, d, e). Several models have been used to explain patterns of LP families during dome growth and explosions including brittle failure of magma in the conduit (Neuberg et al., 2006; Tuffen et al., 2008), stick-slip movements of an ascending plug (Iverson et al., 2006), presence of fluids or fractures in the conduit (Chouet, 1996; Waite et al., 2008), or interaction of magma with the hydrothermal system (Matoza and Chouet, 2010). We are not, however, able to locate the LP families, and we thus cannot speculate on the source mechanism. We surmise that the persistence of the corresponding repeaters indicates that these sources are not immediately modified by the beginning of the eruption and the migration of the crater in 2016–2017 and by the acceleration of dome extrusion in 2019. They could be associated with the conduit but too deep to be affected, similar to activity described by Green and Neuberg (2006) and Thelen et al. (2011) at other volcanoes. Moreover, the long duration of the explosion repeaters in 2019 also suggests that the explosions and the lava extrusion took place in different vents as satellite imagery has shown (Fig. 9).

## 5.5. Interpretative scenario

The bursts of VT activity and the appearance of many new VT repeaters on remote tectonic faults (Fig. 12) a few weeks before the two main eruptive episodes in 2016 and 2019 suggest that these phenomena have a common origin. We cannot rule out that the temporal closeness between the repeaters and the eruptive phenomena is a simple coincidence. However, a possible alternative interpretation can be proposed.

We suggest that the bursts of seismicity and the increased volcanic activity result from processes induced by magma intrusions and occurring simultaneously in the hydrothermal system and in the plumbing system, respectively. Table 2 summarizes the proposed interpretative scenario.

A new magma injection may have started in 2012–2013 in the reservoir located at ~13 km below HHV, producing the wide and regular inflation observed by InSAR and GNSS (Boixart et al., 2020; MacQueen et al., 2020; Fig. S1). The gas generated by the batch of magma rose through a network of fractures in the closed conduit and the surrounding medium and produced the pressurization and heating of the aquifer as revealed by the thermal anomalies and the seismic unrest, including large and distal earthquakes, since 2013. The new magma ascended aseismically by slowly melting the colder material left over in the conduit from the previous eruption in the manner described by Coppola et al. (2022). The ascent rate of the column head depends on the heat flux escaping from the magma, the temperature of the old material filling the conduit, and the conduit radius. In the supporting information, we estimated the column head velocity based on the conceptual model of Coppola et al. (2022) and on their measurements of heat flux. We obtain that, at the end of the pre-eruptive period, the upward velocity could be several hundreds of meters per month. This process of magma progression in the conduit by melting old material is consistent with the low VT seismicity observed below the active crater before the eruption onset.

Around May 2016, the magma head may have encountered a shallow aquifer producing the clear increase of degassing in June reported by Kern et al. (2017) and the appearance of new fumarolic fields between January and August (Fig. 9b). The interaction of magma with water produced a strong pressure increase that propagated as a pore pressure wave through the hydrothermal system following mainly the surrounding connected faults, such as the Sepina fault. The velocity of this pressure wave depends on the permeability ($5 \times 10^{-16}$ to $5 \times 10^{-14}$ m$^2$) and the hydraulic diffusivity ($0.1$ to $10 \text{ m}^2\text{s}^{-1}$) of the fractures and on other factors (Talwani et al., 2007; Mulargia and Bizzarri, 2014, 2015). It could have reached a few km per day in the first days and progressively slowed down later (Mulargia and Bizzarri, 2014). From August to October, the pressure wave reached remote faults and triggered dVT events and repeaters at distances of up to 23 km from the volcano

<table><caption>**Table 2**
Main observations and interpretative scenario.</caption><thead><tr><th>Date</th><th>Scenario</th></tr></thead><tbody><tr><td>2012–2013</td><td>New magma injection in deep reservoir.</td></tr><tr><td>2012–2015</td><td>Large inflation due to 13-km-deep source centered 5 km N from Sabancaya summit</td></tr><tr><td></td><td>Ascent of magma by melting old material in conduit</td></tr><tr><td></td><td>Heating and pressurization of hydrothermal system by rising gases, triggering distal seismic activity</td></tr><tr><td></td><td>Magma head rises with velocity of a few hundreds of meters per month</td></tr><tr><td>May 2016</td><td>Interaction of magma column with aquifer</td></tr><tr><td></td><td>Generation of pore pressure wave and propagation through hydrothermal system</td></tr><tr><td>August–October 2016</td><td>Pressure increases trigger distal VT and repeater seismicity</td></tr><tr><td>2 November</td><td>Breaking of the ductile-brittle transition cap by magma column</td></tr><tr><td>6 November</td><td>Magma column head reaches the surface. Onset of magmatic eruption.</td></tr><tr><td>November 2016</td><td>Opening of magmatic system</td></tr><tr><td></td><td>Decrease of pore pressure in hydrothermal system</td></tr><tr><td></td><td>Reduction of VT and dVT activity</td></tr><tr><td>December 2016</td><td>Fracturing and damaging of shallow structures associated with opening of new crater, producing velocity decrease of about 0.2%</td></tr><tr><td>2019</td><td>New batch of magma rises in conduit</td></tr><tr><td></td><td>Increase of LP and explosive activity and new repeaters</td></tr><tr><td>August–September 2019</td><td>Another pore pressure wave triggers VT, dVT, and repeater activity</td></tr><tr><td>October 2019</td><td>Rate of lava dome extrusion increases</td></tr></tbody></table>
(Fig. 12a). During this period, the magma column continued to rise toward the surface with a velocity on the order of 300 m per month. Stick-slip movements of the magma on the conduit wall produced LP repeaters. The large shallow M<sub>L</sub> 1.6 VT earthquake that occurred on 2 November below the crater was probably due to the fracturing of the impermeable cap by the magma column head, a few days before the eruption began. In November, the opening of the magmatic system produced a decrease of the pore pressure in the hydrothermal system and the reduction of the VT and dVT activity. Then, the progressive opening of the new crater could be associated with the fracturing and damaging of the shallow structure of the volcano as revealed by the velocity decrease of ~0.2% detected in December.

In 2019, a new batch of magma intruded into the conduit producing an increase of the LP and explosive activity, including new repeaters, a velocity decrease of <0.2%, and the acceleration of the lava dome extrusion in October. The interaction of this intrusion with the aquifer generated another pore pressure wave that triggered, after some delay, bursts and repeaters of VT and dVT events (Fig. 12b).

## 6. Conclusions

The long-lived eruption of Sabancaya volcano is an example of complex interactions between magmatic processes, tectonic environment, and hydrothermal system. Its study required multiparametric observations provided by a permanent monitoring network operated by IGP in harsh conditions, as well as the use of state-of-the-art tools for data analysis. Sabancaya and its surrounding region were characterized by a high level of seismic activity before and during the eruptive period, composed mainly of proximal and distal volcano-tectonic events, LP events, and explosion quakes. Using an improved velocity model, we carried out precise hypocenter determinations and identified different types of seismic sequences, earthquakes with magnitudes higher than 4, and long-lived seismogenic zones at distances of up to 30 km from the volcano.

A large part of the seismic catalog belongs to the numerous families of similar waveforms, or repeaters. In particular, two bursts of seismic activity that included many new repeaters were observed a few weeks before the main eruptive events, i.e., the eruption onset in 2016 and the acceleration of lava dome extrusion in 2019. All these observations highlight the major role of the hydrothermal system in the surrounding region. Its progressive heating by magmatic gases produced increases of pore pressure and of Coulomb stress that triggered ruptures on close and remote faults. The sudden interaction of magma with the aquifer generated pore pressure waves that propagated at long distance and triggered bursts of seismicity and distal VT repeaters.

In the conduit, processes of magma convection, heat transfer, and melting of the old material may explain slow rising of magma from the deep reservoir to the surface and the lack of clear increase of VT seismic activity beneath the crater in the days preceding the eruption. The most marked velocity decrease was observed after the eruption onset and may be related to the opening of a new shallow conduit and the migration of the eruptive activity from the old crater to a new crater.

The continuous monitoring and forecasting of eruptive events are challenging at Sabancaya. Much work is still required to understand the high level and wide distribution of the different kinds of seismic activity, the role of the hydrothermal system, the propagation of pressure waves, and the triggering of distal earthquakes. Furthermore, improvements to the model of rising magma, which could take into account mechanical aspects, rock damage, the offset of the deep reservoir north from Sabancaya, and interactions with the aquifer and ductile-brittle transition are needed to better explain velocity variations, LP and explosion activity, and improve interpretations of future activity.

### CRediT authorship contribution statement

Roger Machacca: Conceptualization, Methodology, Software,

Writing – original draft. **Philippe Lesage:** Conceptualization, Methodology, Writing – review & editing. **Hernando Tavera:** Methodology, Investigation, Writing – review & editing. **Jeremy D. Pesicek:** Methodology, Investigation, Writing – review & editing. **Corentin Caudron:** Methodology, Writing – review & editing. **Jose L. Torres:** Data curation, Visualization, Writing – review & editing. **Nino Puma:** Data curation, Visualization, Writing – review & editing. **Katherine Vargas:** Data curation, Visualization, Writing – review & editing. **Ivonne Lazarte:** Data curation, Visualization, Writing – review & editing. **Marco Rivera:** Writing – review & editing. **Alain Burgisser:** Software, Methodology, Writing – review & editing.

### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The authors do not have permission to share data.

### Acknowledgements

R.M. was supported by the Franco-Peruvian doctoral school in Engineering Sciences and Geosciences (EDFPCIG; Grant 67952). Seismic campaigns during 2013 were deployed in cooperation with IRD (Dr. Bernard Valette). Thanks to V. Pinel and the IRD team for the training received in the use of NSBAS. We gratefully acknowledge data and support from the IGP team. The Volcano Disaster Assistance Program receives funding from the U.S. Agency for International Development's Bureau for Humanitarian Assistance and the U.S. Geological Survey's Volcano Hazards Program. The authors gratefully acknowledge the experts who reviewed the manuscript and provided valuable feedback: Weston Thelen, an anonymous reviewer, and Jay Wellik as USGS internal reviewer. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jvolgeores.2023.107767.

### References

Alvarez, Y.V., 2017. Estudio estructural del Sistema Hidrotermal de los volcanes Sabancaya y Hualca-Hualca mediante el método de potencial espontáneo (Tesis para optar el título profesional de Ingeniera Geofísica). Universidad Nacional de San Agustín de Arequipa, Arequipa, Perú.

Ankerst, M., Breunig, M.M., Kriegel, H.P., Sander, J., 1999. OPTICS: Ordering points to identify the clustering structure. In: ACM SIGMOD International Conference on Management of Data. ACM Press, pp. 49–60. CiteSeerX 10.1.1.129.6542.

Antayhua, Y., Tavera, H., Bernal Esquia, I., 2001. Análisis de la actividad sísmica en la región del volcán Sabancaya (Arequipa). Repositorio institucional -IGP.

Antayhua, Y., Tavera, H., Bernal Esquia, I., Palza, H., Aguilar, V., 2002. Localización hipocentrales y características de la fuente de los sismos de Maca (1991), Sepina (1992) y Cabanaconde (1998) Región Volcánica Sabancaya (Arequipa). Repositorio institucional -IGP.

Barberi, F., Corrado, G., Innocenti, F., Luongo, G., 1984. Phlegraean fields 1982–1984: Brief chronicle of a volcano emergency in a densely populated area. Bull. Volcanol. 47, 175–185.

Battaglia, J., Got, J.L., Okubo, P., 2003. Location of long-period events below Kilaua Volcano using seismic amplitudes and accurate relative relocation. J. Geophys. Res. 108 (B12), 2553. https://doi.org/10.1029/2003JB002517.

Baxter, P.J., Searl, A.S., Cowie, H.A., Jarvis, D., Horwell, C.J., 2014. Evaluating the respiratory health risks of volcanic ash at the eruption of the Soufriere Hills Volcano, Montserrat, 1995 to 2010. Geol. Soc. Lond. Mem. 39 (1), 407–425.

Bennington, N., Haney, M., Thurber, C., Zeng, X., 2018. Inferring magma dynamics at Veniaminof volcano via application of ambient noise. Geophys. Res. Lett. 45, 11,650–11,658.
Benoit, J.P., Mcnutt, S.R., 1997. New constraints on source processes of volcanic tremor at Arenal Volcano, Costa Rica, using broadband seismic data. Geophys. Res. Lett. 24 (4), 449–452.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M., Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broadband surface wave dispersion measurements. Geophys. J. Int. 169, 1239–1260.

Boixart, G., Cruz, L.F., Miranda Cruz, R., Euillades, P.A., Euillades, L.D., Battaglia, M., 2020. Source Model for Sabancaya Volcano Constrained by DINsAR and GNSS Surface Deformation Observation. Remote Sens. 12, 1852.

Brenguier, F., Shapiro, N.M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., Neressian, A., 2008. Towards forecasting volcanic eruptions using seismic noise. Nat. Geosci. 1, 126–130.

Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N.M., Briand, X., Emoto, K., Miyake, H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops. Science 345 (6192), 80–82.

Bromley, G.R.M., Thouret, J., Schimmelpfennig, I., et al., 2019. In situ cosmogenic 3He and 36Cl and radiocarbon dating of volcanic deposits refine the Pleistocene and Holocene eruption chronology of SW Peru. Bull. Volcanol. 81, 64. https://doi.org/10.1007/s00445-019-1325-6.

Brown, S.K., Sparks, R.S.J., Mee, K., Vye-Brown, C., Ilyinskaya, E., Jenkins, S.F., Loughlin, S.C., 2015. Country and regional profiles of volcanic hazard and risk. Glob. Volcanic Hazards Risk. https://doi.org/10.1017/CBO9781316276273.

Budi-Santoso, A., Lesage, P., 2016. Velocity variations associated with the large 2010 eruption of Merapi Volcano, Java, retrieved from seismic multiplets and ambient noise cross-correlation. Geophys. J. Int. 206 (1), 221–240.

Carrier, A., Got, J.-L., Peltier, A., Ferrazzini, V., Staudacher, T., Kowalski, P., Boissier, P., 2015. A damage model for volcanic edifices: Implications for edifice strength, magma pressure, and eruptive processes. J. Geophys. Res. Solid Earth 120, 567–583.

Caudron, C., Aoki, Y., Lecocq, T., De Plaen, R., Soubestre, J., Mordret, A., Seydoux, L., Terakawa, T., 2022. Hidden pressurized fluids prior to the 2014 phreatic eruption at Mt Ontake. Nat. Commun. 13, 6145. https://doi.org/10.1038/s41467-022-32252-w.

Chouet, B.A., 1996. Long-period volcano seismicity: its source and use in eruption forecasting. Nature 380, 309–316. https://doi.org/10.1038/380309a0.

Chouet, B.A., Matoza, R.S., 2013. A multi-decadal view of seismic methods for detecting precursors of magma movement and eruption. J. Volcanol. Geotherm. Res. 252, 108–175.

Condori, C., Tavera, H., Marotta, G.S., Marcelo Peres, R., George Sand, F., 2017. 2017. Calibration of the local magnitude scale (ML) for Peru. J. Seismol. 21, 987–999. https://doi.org/10.1007/s10950-017-9647-3.

Coppola, D., Laiolo, M., Cigolini, C., Delle Donne, D., Ripepe, M., 2015. Enhanced volcanic hot-spot detection using MODIS IR data: Results from the MIROVA system. In: Harris, A.J.L., De Groeve, T., Garel, F., Carn, S.A. (Eds.), Detecting, Modelling and Responding to Effusive Eruptions. Geological Society, London. https://doi.org/10.1144/SP426.5. Special Publications 426.

Coppola, D., Valade, S., Masias, P., et al., 2022. Shallow magma convection evidenced by excess degassing and thermal radiation during the dome-forming Sabancaya eruption (2012–2020). Bull. Volcanol. 84, 16. https://doi.org/10.1007/s00445-022-01523-1.

De Plaen, R.S.M., Lecocq, T., Caudron, C., Ferrazzini, V., Francis, O., 2016. Single-station monitoring of volcanoes using seismic ambient noise. Geophys. Res. Lett. 43, 8511–8518.

Donaldson, C., Winder, T., Caudron, C., White, R.S., 2019. Crustal seismic velocity responds to a magmatic intrusion and seasonal loading in Iceland's Northern Volcanic Zone. Sci. Adv. 5 (11).

Ellsworth, W.L., 1995. Characteristic earthquakes and long-term earth-quake forecasts: Implications of Central California seismicity. In: Cheng, F.Y., Sheu, M.S. (Eds.), Urban Disaster Mitigation: The Role of Science and Technology. Elsevier Science Ltd.

Fournier, R.O., 2007. Hydrothermal systems and volcano geochemistry. In: Dzurisin, D. (Ed.), Volcano Deformation. Springer, Berlin Heidelberg.

Frémont, M.J., Malone, S.D., 1987. High-precision relative locations of earthquakes at Mount St. Helens, Washington. J. Geophys. Res. 92 (B10) https://doi.org/10.1029/JB0921B10p10223, 10,223–10,236.

Gerbe, M.C., Thouret, J.C., 2004. Role of magma mixing in the petrogenesis of lavas erupted through the 1990–98 explosive activity of Nevado Sabancaya in South Peru. Bull. Volcanol. 66, 541–561.

Global Volcanism Program, 1988. Report on Sabancaya (Peru). In: McClelland, L. (Ed.), Scientific Event Alert Network Bulletin 13:6. Smithsonian Institution. https://doi.org/10.5479/si.GVP.SEAN198806-354006.

Global Volcanism Program, 1990. Report on Sabancaya (Peru). In: McClelland, L. (Ed.), Bulletin of the Global Volcanism Network, 15:5. Smithsonian Institution. https://doi.org/10.5479/si.GVP.BGVN199005-354006.

Gómez-García, C., Brenguier, F., Boué, P., Shapiro, N.M., Droznin, D.V., Droznina, S.Ya., Senyukov, S.L., Gordeev, E.I., 2018. Retrieving robust noise-based seismic velocity changes from sparse data sets:synthetic tests and application to Klyuchevskoy volcanic group (Kamchatka). Geophys. J. Int. 214 (2), 1218–1236.

Got, J.L., Amitrano, D., Stefanou, I., Brothelande, E., Peltier, A., 2019. Damage and Strain Localization around a Pressurized Shallow-Level Magma Reservoir. J. Geophys. Res. Solid Earth 124 (2), 1443–1458.

Green, D.N., Neuberg, J., 2006. Waveform classification of volcanic low-frequency earthquake swarms and its implication at Soufrière Hills Volcano, Montserrat. J. Volcanol. Geotherm. Res. 153, 51–63.

Greenfield, T., White, R.S., Winder, T., Ágústsdóttir, T., 2018. Seismicity of the Askja and Bárðarbunga volcanic systems of Iceland, 2009–2015. J. Volcanol. Geotherm. Res. 391, 1–16.

Hanks, T.C., Kanamori, H., 1978. Moment Magnitude Scale No. 12 Transactions American Geophysical Union Vol. 59. Amer Geophysical Union, 2000 Florida Ave NW, Washington, DC 20009, p. 1128.

Hautmann, S., Witham, F., Christopher, T., Cole, P., Linde, A.T., Sacks, I.S., Sparks, R.S. J., 2014. Strain field analysis on Montserrat (W.I.) as tool for assessing permeable flow paths in the magmatic system of Soufrière Hills Volcano. Geochem. Geophys. Geosyst. 15 (3), 676–690.

Hidalgo, S., Battaglia, J., Arellano, S., Steele, A., Bernard, B., Bourquin, J., et al., 2015. SO2 degassing at Tungurahua volcano (Ecuador) between 2007 and 2013: transition from continuous to episodic activity. J. Volcanol. Geotherm. Res. 298, 1–14. https://doi.org/10.1016/j.jvolgeores.2015.03.022.

Hotovec-Ellis, A.J., Jeffries, C., 2016. Near Real-time Detection, Clustering, and Analysis of Repeating Earthquakes: Application to Mount St. Helens and Redoubt Volcanoes – Invited, presented at Seismological Society of America Annual Meeting, Reno, Nevada, 20 Apr.

Hotovec-Ellis, A.J., Shiro, B.R., Shelly, D.R., Anderson, K.R., Haney, M.M., Thelen, W.A., et al., 2022. Earthquake-derived seismic velocity changes during the 2018 caldera collapse of Kilauea volcano. J. Geophys. Res. Solid Earth 127 e2021JB023234.

Huaman-Rodrigo, D., Chorowicz, J., Deffontaines, B., Guillande, R., Rudant, J., 1993. Cadre structural et risques géologiques étudiés à l'aide de l'imagerie spatiale: La région du Colca (Andes du sud pérou). Bulletin-Societegeologiquedefrance 164, 807–807.

Hutton, L.K., Boore, D.M., 1987. The ML scale in southern California. Bull. Seismol. Soc. Am. 77, 2074–2094.

Iverson, R.M., Dzurisin, D., Gardner, C.A., Gerlach, T.M., LaHusen, R.G., Lisowski, M., Major, J.J., Malone, S.D., Messerich, J.A., Moran, S.C., Pallister, J.S., Qamar, A.I., Schilling, S.P., Vallance, J.W., 2006. Dynamics of seismogenic volcanic extrusion at Mount St. Helens in 2004–05. Nature 444, 439–443.

Jay, J.A., Delgado, F.J., Torres, J.L., Pritchard, M.E., Macedo, O., Aguilar, V., 2015. Deformation and seismicity near Sabancaya volcano, southern Peru, from 2002 to 2015. Geophys. Res. Lett. 42, 2780–2788. https://doi.org/10.1002/2015GL063589.

Johnson, J., Aster, R., 2005. Relative partitioning of acoustic and seismic energy during Strombolian eruptions. J. Volcanol. Geotherm. Res. 148, 334–354.

Jones, R., Manville, V., Andrade, D., 2015. Probabilistic analysis of rain-triggered lahser initiation at Tungurahua volcano. Bull. Volcanol. 77 (8), 68.

Kennett, B., Engdahl, E., 1991. Traveltimes for global earthquake location and phase identification. Geophys. J. Int. 105 (2), 429–465. https://doi.org/10.1111/j.1365-246X.1991.tb06724.x.

Kern, C., Masias, P., Apaza, F., Reath, K.A., Platt, U., 2017. Remote measurement of high pre-eruptive water vapor emissions at Sabancaya volcano by passive differential optical absorption spectroscopy. J. Geophys. Res. Solid Earth 122, 3540–3564. https://doi.org/10.1002/2017JB014020.

Kissling, E., Ellsworth, W.L., Eberhart-Phillips, D., Kradolfer, U., 1994. Initial reference models in local earthquake tomography. J. Geophys. Res. 99.

Kissling, E., Kradolfer, U., Maurer, H., 1995. VELEST user's Guide-Short Introduction. Institute of Geophysics and Swiss Seismological Service, ETH, Zurich, p. 25.

Lacroix, P., Berthier, J., Maquerhua, E.T., 2015. Earthquake-driven acceleration of slow-moving landslides in the Colca valley, Peru, detected from Pléiades images. Remote Sens. Environ. 165, 148–158. https://doi.org/10.1016/j.rse.2015.05.010.

Lahr, J.C., Chouet, B.A., Stephens, C.D., Power, J.A., Page, R.A., 1994. Earthquake classification, location, and error analysis in a volcanic environment: implications for the magmatic system of the 1989–1990 eruptions at Redoubt Volcano, Alaska. J. Volcanol. Geotherm. Res. 62 (1–2), 137–151.

Lengliné, O., Duputel, Z., Ferrazzini, V., 2016. Uncovering the hidden signature of a magmatic recharge at Piton de la Fournie volcano using small earthquakes. Geophys. Res. Lett. 43 (9), 4255–4262.

Lesage, P., Mora, M.M., Alvarado, G.E., Pacheco, J., Metaxian, J.P., 2006. Complex behavior and source model of the tremor at Arenal volcano, Costa Rica. J. Volcanol. Geotherm. Res. 157, 49–59.

Lin, C.-H., 2017. Dynamic triggering of volcano drumbeat-like seismicity at the Taton volcano group in Taiwan. Geophys. J. Int. 210 (1), 354–359.

Lobkis, O.I., Weaver, R.L., 2003. Coda-wave interferometry in finite solids: Recovery of Pto-S conversion rates in an elastodynamic billiard. Phys. Rev. Lett. 90 (25), 254302.

Lomax, A., Virieux, J., Volant, P., Berge, C., 2000. Probabilistic earthquake location in 3D and layered models: Introduction of a Metropolis-Gibbs method and comparison with linear locations. In: Thurber, C.H., Rabinowitz, N. (Eds.), Advances in Seismic Event Location. Kluwer, Amsterdam, pp. 101–134.

Loughlin, S.C., Vye-Brown, C., Sparks, R.S.J., Brown, S.K., Barclay, J., Calder, E., et al., 2015. Chapter 1 an Introduction to Global Volcanic Hazard and Risk. Cambridge University Press.

Ma, Y., Clayton, R.W., 2014. The crust and uppermost mantle structure of Southern Peru from ambient noise and earthquake surface wave analysis. Earth Planet. Sci. Lett. 395, 61–70.

Macedo Sánchez, O., Taipe Maquerhua, E.L., Del Carpio, J.A., Ticona Paucara, J., Ramos Palomino, D.A., Puma, N., Aguilar Purhuaya, V., Machacca Puma, R., Torres, J.L., Sandoval, C., et al., 2016. Evaluación del riesgo volcánico en el sur del Perú, situación de la vigilancia actual y requerimientos de monitoreo en el futuro. Informe técnico. Informes Técnicos Instituto Geofísico del Perú (IGP), 105, pp.

Machacca, R., Del Carpio, J., Rivera, M., Tavera, H., Macedo, L., Concha, J.A., Lazarte, I. A., Centeno, R.G., Puma, N.C., Torres, J.L., Vargas, K.A., Cruz, J.E., Velarde, L., Vilca, J., Malpartida, A.R., 2021. Monitoring of active volcanoes in Peru by the Instituto Geofísico del Perú: early warning systems, communication, and information dissemination. Volcanica 4 (S1), 49–71.

Machacca-Puma, R., Lesage, P., Larose, E., Lacroix, P., Ancacci-Figueroa, R.M., 2019. Detection of pre-eruptive seismic velocity variations at an andesitic volcano using
ambient noise correlation on 3-component stations: Ubinas volcano, Peru. J. Volcanol. Geotherm. Res. 381 (83–100), 1–75.

Machare, J., Fenton, C.H., Machette, M.N., Lavenu, A., Costa, C., Dart, R.L., 2003. Database and Map of Quaternary Faults and Folds in Peru and its Offshore Region, USGS Numbered Series 2003–451: United States Geological Survey.

MacQueen, P., Delgado, F., Reath, K., Pritchard, M.E., Bagnardi, M., Milillo, P., et al., 2020. Volcano-ectonic interactions at Sabancaya volcano, Peru: Eruptions, magmatic inflation, moderate earthquakes, and fault creep. J. Geophys. Res. Solid Earth 125. https://doi.org/10.1029/2019JB019281 e2019JB019281.

Massin, F., Farrell, J., Smith, R.B., 2013. Repeating earthquakes in the Yellowstone volcanic field: Implications for rupture dynamics, ground deformation, and migration in earthquake swarms. J. Volcanol. Geotherm. Res. 257, 159–173.

Matoza, R.S., Chouet, B.A., 2010. Subevents of long-period seismicity: implications for hydrothermal dynamics during the 2004–2008 eruption of Mount St. Helens. J. Geophys. Res. 115. B12206.

Matoza, R.S., Chouet, B.A., Dawson, P.B., Shearer, P.M., Haney, M.M., Waite, G.P., Moran, S.C., Mikewell, T.D., 2015. Source mechanism of small long-period events at Mount St. Helens in July 2005 using template matching, phase-weighted stacking, and full-waveform inversion. J. Volcanol. Geotherm. Res. 120, 6351–6364.

McNutt, S.R., 1996. Seismic monitoring and Eruption forecasting of Volcanoes: A Review of the State-of-the-Art and Case Histories. In: Monitoring and Mitigation of Volcano Hazards. Springer, Berlin Heidelberg, Berlin, Heidelberg, pp. 99–146. https://doi.org/10.1007/978-3-642-80087-0_3.

McTigue, D.F., 1987. Elastic stress and deformation near a finite spherical magma body: resolution of the point source paradox. J. Geophys. Res. 92 (B12), 12931–12940.

Mering, C., Huaman-Rodrigo, D., Chorowicz, J., Deffontaines, B., Guillande, R., 1996. New data on the geodynamics of southern Peru from computerized analysis of SPOT and SAR ERS-1 images. Tectonophysics 259 (1), 153–169. https://doi.org/10.1016/0040-1951(96)00034-0.

Meyer, K., Biggs, J., Aspinall, W., 2021. A Bayesian reassessment of the relationship between seismic moment and magmatic intrusion volume during volcanic unrest. J. Volcanol. Geotherm. Res. 419, 107375.

Miao, Q., Langston, C.A., 2007. Empirical distance attenuation and the local-magnitude scale for the Central United States. Bull. Seismol. Soc. Am. 97, 2137–2151.

Mordret, A., Jolly, A.D., Duputel, Z., Fournier, N., 2010. Monitoring of phreatic eruptions using interferometry on retrieved cross-correlation function from ambient seismic noise: results from Mt. Ruapehu, New Zealand. J. Volcanol. Geotherm. Res. 191, 46–59.

Mothes, P.A., Yepes, H.A., Hall, M.L., Ramón, P.A., Steele, A.L., Ruiz, M.C., 2015. The scientific–community interface over the fifteen-year eruptive episode of Tungurahua Volcano Ecuador. J. Appl. Volcanol. 4 (1), 1–15.

Mulargia, F., Bizzarri, A., 2014. Anthropogenic triggering of large earthquakes. Nat. Sci. Rep. 4, 6100. https://doi.org/10.1038/srep06100.

Mulargia, F., Bizzarri, A., 2015. Fluid pressure waves trigger earthquakes. Geophys. J. Int. 200 (3), 1279–1283. https://doi.org/10.1093/gji/ggu469.

Nadeau, R.M., Foxall, W., McEvilly, T.V., 1995. Clustering and periodic recurrence of microearthquakes on the San Andreas Fault at Parkfield, California. Science 267, 503–507.

Neuberg, J.W., Tuffen, H., Collier, L., Green, D., Powell, T., Dingwell, D., 2006. The trigger mechanism of low-frequency earthquakes on Montserrat. J. Volcanol. Geotherm. Res. 153, 37–50.

Newhall, C.G., 2000. Volcano warnings. In: Encyclopedia of Volcanoes. Academic Press, New York, pp. 1–85.

Pesicex, J.D., Ogburn, S.E., Prejean, S.G., 2021. Indicators of volcanic eruptions revealed by global M4+ earthquakes. J. Geophys. Res. Solid Earth 126. https://doi.org/10.1029/2020JB021294 e2020JB021294.

Phillipson, G., Sobradelo, R., Gottsmann, J., 2013. Global volcanic unrest in the 21st century: an analysis of the first decade. J. Volcanol. Geotherm. Res. 264, 183–196. https://doi.org/10.1016/j.jvolgeores.2013.08.004.

Potter, S.H., Scott, B.J., Jolly, G.E., Neall, V.E., Johnston, D.M., 2015. Introducing the volcanic unrest index (VUI): a tool to quantify and communicate the intensity of volcanic unrest. Bull. Volcanol. 77 (9).

Power, J.A., Lalla, D.J., 2010. “Seismic Observations of Augustine Volcano, 1970–2007,” in the 2006 Eruption of Augustine Volcano, Alaska. J. A. Power, M. L. Coombs, and J. T. Freymueller (U.S. Geological Survey Professional Paper 1769).

Puma, N., Macedo, O., Alvarez, Y., Finizola, A., Ramos, D., 2018. Estudio estructural y del sistema hidrotermal de los volcanes Sabancaya y Hualca-Hualca mediante el método de Potencial Espontáneo. In: VIII Foro Internacional: “Los volcanes y su impacto”, Arequipa, Perú.

Reath, K.A., Ramsey, M.S., Dehn, J., Webley, P.W., 2016. Predicting eruptions from precursory activity using remote sensing data hybridization. J. Volcanol. Geotherm. Res. 321, 18–30.

Richter, C.F., 1935. An instrumental earthquake magnitude scale. Bull. Seismol. Soc. Am. 25, 1–32.

Ridley, D.A., Solomon, S., Barnes, J.E., Burlakov, V.D., Deshler, T., Dolgii, S.I., et al., 2014. Total volcanic stratospheric aerosol optical depth and implications for global climate change. Geophys. Res. Lett. 41, 7763–7769. https://doi.org/10.1002/2014GL061541.

Ripepe, M., Harris, A.J., Carniel, R., 2002. Thermal, seismic and infrasmonic evidences of variable degassing rates at Stromboli volcano. J. Volcanol. Geotherm. Res. 118, 285–297. https://doi.org/10.1016/S0377-0273(02)00298-6.

Ripepe, M., Marchetti, E., Ulivieri, G., Harris, A., Dehn, J., Burton, M., Caltabiano, T., Salerno, G., 2005. Effusive to explosive transition during the 2003 eruption of Stromboli volcano. Geology 33, 341–344. https://doi.org/10.1130/G21173.1.

Rodríguez, A., Uribe, M., 1994. Participación del Instituto Geofísico del Perú en relación con la reactivación del volcan Sabancaya, Provincia de Caylloma, Región Arequipa, Informe División de Sismología y Gravimetría, IGP-Arequipa.

Rowe, C.A., Thurber, C.H., White, R.A., 2004. Domegrowthbehavior at Soufrière Hills Volcano, Montserrat, revealed by relocation of volcanic event swarms, 1995–1996. J. Volcanol. Geotherm. Res. 134 (3), 199–221. https://doi.org/10.1016/j.jvolgeores.2004.01.008.

Saccorotti, G., Petrosino, S., Bianco, F., Castellano, M., Galluzzo, D., et al., 2007. Seismicity associated with the 2004–2006 renewed ground uplift at Campi Flegrei Caldera, Italy. Phys. Earth Planet. Inter. 165, 1424.

Samaniego, P., Rivera, M., Mariño, J., Guillou, H., Liorzou, C., Zerathe, S., et al., 2016. The eruptive chronology of the Ampato-Sabancaya volcanic complex (Siamos Peru). J. Volcanol. Geotherm. Res. 323, 110–128. https://doi.org/10.1016/j.jvolgeores.2016.04.038.

Sébrier, M., Mercier, J.L., Mégard, F., Laubacher, G., Carey-Gailhardis, E., 1985. Quaternary normal and reverse faulting and the state of stress in the Central Andes of South Peru. Tectonics 4 (7), 739–780. https://doi.org/10.1029/TC004i007p00739.

Segall, P., 2013. Volcano deformation and eruption forecasting. In: Geological Society, 380. Special Publications, London, pp. 85–106.

Sens-Schönfelder, C., Wegler, U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. Geophys. Res. Lett. 33, L21302. https://doi.org/10.1029/2006GL027797.

Shapiro, N.M., Droznin, D.V., Droznina, S.Ya., Senyukov, S.L., Gusev, A.A., Gordeev, E.I., 2017. Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer. Nat. Geosci. 10 (6), 442–445.

Siebert, L., Simkim, T., Kimberley, P., 2010. Volcanoes of the World, Third edition. Smithsonian Institution, Washington, University of California.

Sparks, R., 2003. Forecasting volcanic eruptions. Earth Planet. Sci. Lett. 210 (1–2), 1–15. https://doi.org/10.1016/S0012821X(03)00124-9.

Talwani, P., Chen, L., Gahalaut, K., 2007. Seismogenic permeability, ks. J. Geophys. Res. 112, B07309. https://doi.org/10.1029/2006JB004665.

Tarantola, A., Valette, B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion. Rev. Geophys. Space Phys. 28, 219–232.

Tavaera, H., Guardia, P., Condori, C., Fernandez, E., Arredondo, L., 2013. Sismos de la Región del volcán Sabancaya del 22 y 23 de febrero del 2013. In: Informe Técnico N° 01–2013. Instituto Geofísico del Perú.

Thelen, W., Malone, S., West, M., 2011. Multiplets: their Behavior and Utility at Dacitic and Andesitic Volcanic Centers. J. Geophys. Res. 116 (B8), B08210. https://doi.org/10.1029/2010JB007924.

Thouret, J.C., Guillande, R., Huamán, D., Bourgaud, A., Salas, G., Chorowicz, J., 1994. L'activité actuelle du Nevado Sabancaya (Sud Pérou): reconnaissance géologique et satellitaire, évaluation et cartographie des menaces volcaniques. Bull. Soc. Géol. France 1 (165), 49–63.

Thouret, J.C., Wörner, G., Gunnell, Y., Singer, B., Zhang, X., Souriot, T., 2007. Geochronologic and stratigraphic constraints on canyon incision and Miocene uplift of the Central Andes in Peru. Earth Planet. Sci. Lett. 263, 151–166. https://doi.org/10.1016/j.epsl.2007.07.023.

Torres, R.A., Gomez, D.M., Narvaez, M.L., 1996. Unusual seismic signals associated with the activity at Galeras volcano, Colombia, from July 1992 to September 1994. Ann. Geophys. 39, 299–310.

Travada, Y., Córdova, V., 1752. El Suelo de Arequipa Convertido en Cielo. In: de Odriozola, M. (Ed.), 1877. Documentos Literarios del Perú. Tomo X. Imprenta del Estado, Lima, p. 415.

Tuffen, H., Smith, R., Sammonds, P.R., 2008. Evidence for seismogenic fracture of silicic magma. Nature 453, 511–514.

Vidale, J.E., Ellsworth, W.L., Cole, A., Marone, C., 1994. Variations in rupture process with recurrence interval in a repeated small earthquake. Nature 368, 624–626. https://doi.org/10.1785/0120070026.

Waite, G.P., Chouet, B.A., Dawson, P.B., 2008. Eruption dynamics at Mount St. Helens imaged from broadband seismic waveforms: interaction of the shallow magmatic and hydrothermal systems. J. Geophys. Res. 113 https://doi.org/10.1029/2007JB005259. B02305.

Wellik, J.J., Prejean, S.G., Syahbana, D.K., 2021. Repeating Earthquakes during Multiple Phases of Unrest and Eruption at Mount Agung, Bali, Indonesia, 2017. Front. Earth Sci. 9, 653164 https://doi.org/10.3389/feart.2021.653164.

White, R.A., McCausland, W., 2016. Volcano-tectonic earthquakes: a new tool for estimating intrusive volumes and forecasting eruptions. J. Volcanol. Geotherm. Res. 309, 139–155.

White, R.S., Drew, J., Martens, H.R., Key, J., Soosalu, H., Jakobsdóttir, S.S., 2011. Dynamics of Dyke Intrusion in the Mid-crust of Iceland. Earth Planet. Sci. Lett. 304, 300–312. https://doi.org/10.1016/j.epsl.2011.02.038.

Wooster, M.J., Zhukov, B., Oertel, D., 2003. Fire radiative energy for quantitative study of biomass burning: derivation from the BIRD experimental satellite and comparison to MODIS fire products. Remote Sens. Environ. 86 (1), 83–107.

Yates, A.S., Savage, M.K., Jolly, A.D., Caudron, C., Hamling, I.J., 2019. Volcanic, coseismic, and seasonal changes detected at White Island (Whakaari) volcano, New Zealand, using seismic ambient noise. Geophys. Res. Lett. 46 https://doi.org/10.1029/2018GL080580.

Zamácola, Y., Jaúregui, J.D., 1888. Apuntes Para la Historia de Arequipa, Año de 1804. Imprenta de la Bolsa, Arequipa, p. 96.

Zelenski, M., Taran, Y., Galle, B., 2015. High emission rate of sulfuric acid from Bezmianny volcano, Kamchatka. Geophys. Res. Lett. 42, 7005–7013.