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

HAL
open science

# Shallow magma convection evidenced by excess degassing and thermal radiation during the dome-forming Sabancaya eruption (2012–2020)

Diego Coppola, S Valade, P Masias, M Laiolo, F Massimetti, A Campus, R Aguilar, R Anccasi, F Apaza, B Ccallata, et al.

► To cite this version:

Diego Coppola, S Valade, P Masias, M Laiolo, F Massimetti, et al. Shallow magma convection evidenced by excess degassing and thermal radiation during the dome-forming Sabancaya eruption (2012–2020). Bulletin of Volcanology, 2022, 84 (16), 10.1007/s00445-022-01523-1. hal-03543369

HAL Id: hal-03543369

https://hal.univ-reunion.fr/hal-03543369

Submitted on 26 Jan 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Shallow magma convection evidenced by excess degassing and thermal radiation during the dome-forming Sabancaya eruption (2012–2020)

D. Coppola<fnref n="1" />,<fnref n="2" /> · S. Valade<fnref n="3" /> · P. Masias<fnref n="4" /> · M. Laiolo<fnref n="1" />,<fnref n="2" /> · F. Massimetti<fnref n="1" /> · A. Campus<fnref n="1" /> · R. Aguilar<fnref n="4" /> · R. Anccasi<fnref n="4" /> · F. Apaza<fnref n="4" /> · B. Ccallata<fnref n="4" /> · C. Cigolini<fnref n="1" />,<fnref n="2" /> · L. F. Cruz<fnref n="4" />,<fnref n="5" /> · A. Finizola<fnref n="6" />,<fnref n="7" /> · K. Gonzales<fnref n="4" /> · O. Macedo<fnref n="5" /> · R. Miranda<fnref n="4" /> · M. Ortega<fnref n="4" /> · R. Paxi<fnref n="4" /> · E. Taipe<fnref n="4" /> · D. Valdivia<fnref n="4" />

Received: 31 August 2021 / Accepted: 30 December 2021

© The Author(s) 2022

## Abstract

We used a large set of satellite- (visible, infrared, and radar images from Planetscope, MODIS, VIIRS, Sentinel2, Landsat 8, and Sentinel 1) and ground-based data (optical images, SO₂ flux, shallow seismicity) to describe and characterize the activity of the Sabancaya volcano during the unrest and eruption phases that occurred between 2012 and 2020. The unrest phase (2012–2016) was characterized by increasing gas and thermal flux, sourced by a convective magma column rising along with the remnants of a buried plug still permeable to fluid flow. Conversely, a new conduit, adjacent to the previous one, fed the eruptive phase (2016–2020) which was instead characterized by a discontinuous extrusive activity, with phases of dome growth (at rates from 0.04 to 0.75 m<fnref n="3" /> s⁻<fnref n="1" />) and collapse. The extrusive activity was accompanied by fluctuating thermal anomalies (0.5–25 MW), by irregular SO₂ degassing (700–7000 tons day⁻<fnref n="1" />), and by variable explosive activity (4–100 events d⁻<fnref n="1" />) producing repeated vulcanian ash plumes (500–5000 m above the crater). Magma budget calculation during the eruptive phase indicates a large excess of degassing, with the volume of degassed magma (0.25–1.28 km<fnref n="3" />) much higher than the volume of erupted magma (<0.01 km<fnref n="3" />). Similarly, the thermal energy radiated by the eruption was much higher than that sourced by the dome itself, an unbalance that, by analogy with the degassing, we define as “excess thermal radiation”. Both of these unbalances are consistent with the presence of shallow magma convection that fed the extrusive and explosive activity of the Sabancaya dome.

**Keywords** Sabancaya · Dome-forming-eruption · Excess degassing · Excess radiation · Shallow magma convection

## Introduction

Lava-dome-forming eruptions are one of the most spectacular volcanic manifestations, as well as among the most dangerous (Newhall and Melson, 1983; Calder et al., 2005). These eruptions may have a wide spectrum of characteristics

Editorial responsibility: Y. Moussallam

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

D. Coppola
diego.coppola@unito.it

<footnote n="1">1. 1Dipartimento Di Scienze Della Terra (DST), Università Degli Studi Di Torino, Torino, Italy</footnote>

<footnote n="2">2. 2Centro Interdipartimentale Sui Rischi Naturali in Ambiente Montano E Collinare (NATRISK), Università Degli Studi Di Torino, Torino, Italy</footnote>

<footnote n="3">3. 3Departamento de Vulcanología, Instituto de Geofísica, Universidad Nacional Autónoma de México, Mexico City, Mexico</footnote>

<footnote n="4">4. 4Observatorio Vulcanológico del INGEMMET (OVI), Arequipa, Peru</footnote>

<footnote n="5">5. 5Universidad Nacional de San Agustín de Arequipa, Arequipa, Peru</footnote>

<footnote n="6">6. 6Laboratoire GéoSciences Réunion, 6Université de La Réunion, 97744 Saint Denis, France</footnote>

<footnote n="7">7. 7Université de Paris, Institut de Physique du Globe de Paris, CNRS, 75005 Paris, France</footnote>
ranging from the gentle effusion of a viscous lava body (Eichelberger et al., 1986) to catastrophic collapses causing Pyroclastic Density Currents (PDCs), explosive eruptions, and sudden volcanic blasts (e.g. Nakada et al. 1999; Vallance et al. 2008; Calder et al. 2002; Carn et al. 2004; Ogburn et al., 2015). They can trigger cascade risks, including lahars or tsunamis, and very often can suddenly change the eruptive style from effusive to explosive (Calder et al., 2002; Boudon, 2015; Pararas-Carayannis 1992; Plank et al., 2019). For these reasons, the dynamics of lava dome extrusion have been extensively studied (Calder et al., 2005 and reference therein). Yet, several processes remain unclear, most likely because they are also dependent on external forcing (such as heavy rains or large earthquakes; Elsworth et al., 2004; Walter et al., 2008), volcano-tectonic or local conditions (Walter et al., 2015; Zorn et al., 2019), or nonlinear processes affecting the magmatic plumbing system (Melnik and Sparks, 1999). The latter include fluctuations in the chamber pressure, cyclicity in magma discharge rate and/or degassing, as well as the mutual interactions of these parameters with the changing properties of magma (i.e. viscosity, crystallinity, permeability) along the conduit (Sparks 1997; Melnik and Sparks, 1999; Melnik and Sparks; 2005; Costa et al., 2007).

It has been widely recognized that many dome-forming eruptions are also characterized by excess degassing (Andres et al., 1991; Wallace, 2001), implying that the lava dome is connected to the chamber through (i) a gas-permeable volcanic conduit, or (ii) a convecting magma column (Shinohara 2008 and references therein). These different degassing mechanisms can strongly influence the parameters typically monitored on a volcano, such as seismicity, thermal activity, deformation, and have an obvious control on the surface activity and its potential fluctuations.

A fundamental challenge of volcanic monitoring is, therefore, to better understand the processes at the origin of the monitored parameters and recognize mechanisms or short-to intermediate-term conditions that can lead to lava dome failure and/or transition between eruptive styles (Pallister and McNutt, 2015; Moussallam et al., 2021). For this purpose, the correlation between acquired data and eruptive/degassing processes is fundamental, although not always trivial because it is dependent on both the number of monitored parameters and the conceptual model used to interpret these data (Pallister and McNutt, 2015; Reath et al., 2020).

In many high-altitude, dome-forming volcanoes, it is very difficult to obtain continuous and homogeneous information about the eruptive activity in progress, especially because of the harsh environmental conditions which make the installation of permanent monitoring networks arduous and dangerous (Pallister and McNutt, 2015). In these cases, satellite and ground data very often complement each other as they offer a comprehensive, synoptic vision of the ongoing eruptive

phenomena (Furtney et al., 2018; Laiolo et al., 2019; Reath et al., 2019a, b).

In this work, we present a unique, multiparametric dataset, consisting of satellite- and ground-based measurements, that allowed us to characterize the precursory and eruptive activity of Sabancaya volcano (Peru), in the period 2012–2020. Specifically, the large set of acquired data was used to analyze the shallow processes associated with this eruption, which include (i) thermal emission, (ii) gas emission (SO₂), (iii) phases of dome growth/collapse, (iv) height of the plume, and (v) shallow seismicity associated with explosive events.

This dataset is finally used to calculate the magma volumes involved in different processes at the origin of the observed activity (i.e. extrusion, degassing, thermal and explosive activity) which in turn allow constraining the shallow magma budget of the Sabancaya eruption.

## Sabancaya volcano

Sabancaya (15°47'S; 71°51'W) is a 12-ky-old Peruvian stratovolcano, located in the Central Andean Zone, a volcanic zone originated by the oceanic-continental collision and the consequent subduction of the Nazca Plate under the South American Plate. The 5976-m high Sabancaya belongs to the Ampato-Sabancaya volcanic complex and is constituted by a succession of lava flows and pyroclastic covering an approximate area of 70 km² (Samaniego et al., 2016). The Sabancaya activity probably started 6–10 ky ago and built a volcanic edifice of 6–10 km³ in volume. The top of the edifice, covered by ice and snow during a large part of the year, hosts an active crater of approximately 350 m in diameter (Fig. 1), responsible for the Holocene eruptive activity characterizing the volcano. This last is characterized by low-to-moderate magnitude (VEI 1–3) vulcanian or phreatomagmatic activity, interspersed by periods of repose lasting decades to centuries (Global Volcanism Program, 2013; Juvigné et al., 1998; 2008; Samaniego et al., 2016), with erupted products petrochemically ranging from andesite to dacite (Gerbe and Thouret, 2004; Samaniego et al., 2016).

In 1988, after about 200 years of rest, Sabancaya entered a new eruptive phase, which lasted until 1997 (Global Volcanism Program, 1988; Gerbe and Thouret, 2004). This eruption was characterized by mild vulcanian explosions (VEI 1–2) accompanied by small eruptive columns (5–7-km height; Gerbe and Thouret, 2004). The eruption reached a climax in May–October 1990 during which approximately 0.025 km³ of tephra were deposited (Gerbe and Thouret, 2004). In the following years, a gradual decrease of the activity was observed with a magma production rate of 0.001–0.01 km³ per year. Notably, satellite images acquired during the eruption have shown the presence of an open crater in the summit

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 1** a Geological setting of Sabancaya volcano, and location of the monitoring network used in this work; b An eruptive column 4500 m height imaged by the SCOP camera on 24 August 2018. c</figcaption></figure>

Detailed view of the lava dome surface on 29 October 2019. By this time the dome measured approximately 230×290 m

area, but no evidence of a lava dome was observed throughout the 1988–1997 eruption (Gerbe and Thouret, 2004).

After 15 years of quiescence, the reactivation of Sabancaya became clear in late 2012, when increased fumarolic activity was accompanied by frequent seismic swarms and

deformation (Global Volcanism Program, 2013; Jay et al., 2015; Moussallam et al., 2017; Kern et al., 2017; Boixart et al., 2020; MacQueen et al., 2020). On November 6, 2016, a series of vulcanian explosions marks the beginning of a new eruptive phase of Sabancaya (MacQueen et al., 2020),
still ongoing at the time of this writing. The sequence of
November 2016 explosions was followed by the extrusion
of the first lava dome in more than 300 years.

Petrological and geophysical studies (Gerbe & Thouret, 2004; Jay et al., 2015; Boixart et al., 2020; MacQueen et al., 2020) suggest that the main magmatic reservoir of Sabancaya is off-centered with respect to the summit, and is located about 7 km NE (beneath Hualca Hualca) at ~13 km depth. This reservoir is thought to be supplied by an andesitic mafic magma which causes it to inflate at a rate of 0.02–0.05 km³ yr⁻¹, and destabilize regional faults (Boixart et al., 2020; MacQueen et al., 2020). Gerbe & Thouret (2004) also suggest that during the eruptive phases a shallow magma chamber (located under the Sabancaya at 6 km depth), is temporarily fed by continuous or intermittent arrival of magma from the deeper reservoir. Hence, it is supposed that heat and magmatic gases from this shallow magma chamber travel along pre-existing permeable pathways in the volcanic conduit of Sabancaya, causing thermal anomalies and elevated degassing from the summit crater (Moussallam et al., 2017; Jay et al., 2015; MacQueen et al., 2020).

Dataset

To characterize the Sabancaya eruption between 2012 and
2020, we used eight different datasets which are described
below and grouped into (i) satellite-based and (ii) ground-
based measurements.

Satellite-based dataset

Lava dome area from high-resolution visible images
(PLANETSCOPE, CNES/Airbus)

Constellations of small satellites named “CubeSats” are opening new perspectives in the earth observation panorama, due to their both high spatial resolution and almost daily revisit frequency (Cigna et al., 2020). In this work, we used PlanetScope imagery (on “Doves” CubeSat satellites; https://earth.esa.int/eogateway/missions/planetscope), acquiring visible images with a spatial resolution of ~3 m and a temporal resolution of ~1–72 h (Adelghi et al. 2019). Access to the online Planet web platform (https://www.planet.com/explorer) allowed us to observe the evolving dome features, manually draw its perimeter, and calculate the dome area (Fig. 3c) for cloud-free images where a change was appreciable. Similarly, two CNES/Airbus images acquired on 27 May 2016 and 17 May 2019, and made available via the Google Earth platform©, were used to visualize in greater detail (0.5 m res.) the morphology of the pre- and sin-eruption crater (Fig. 2a).

Crater Depth and lava dome volume (SENTINEL 1)

SENTINEL-1 (S1) Synthetic Radar Aperture (SAR) intensity images were processed by the volcano monitoring system MOUNTS (Valade et al. 2019, www mounts-project.com) to analyze the dome growth. Images with VV polarization from the ascending orbit were processed (flight path direction N°-12, relative orbit #47); the radar was therefore looking at Sabancaya crater with an ENE direction, i.e., with a Line-of-Sight (LOS) direction of N°78 (Fig. 2c). The temporal resolution of the acquisitions on this orbit is 1 image every 12 days, while the spatial resolution of the images is ~14.1 m in the azimuth direction (i.e., along the flight path), and 2.3 m in range direction (i.e., along the LOS). Crater depth estimates were obtained from SAR images in radar geometry by measuring the length of the shadow cast by the western crater wall and multiplying it by the cosine of the radar incidence angle (Wadge et al., 2011). The shadow length decreases as the dome grows, thereby providing the crater depth (Fig. 3d), an indirect measure of the dome height. The combination of the area of the dome (“Lava dome area from high-resolution visible images (PLANETSCOPE, CNES/Airbus)” section) and its height (from Sentinel 1 data) allowed us to estimate the lava dome volume, by assuming an inverted truncated cone geometry having a minor base diameter of 75±25 m (consistent with the initial area of the dome measured on 5 February 2017; see next chapter).

Hot area and spatial distribution of thermal
anomalies (SENTINEL 2 and LANDSAT 8)

SENTINEL-2 MultiSpectral Instrument and LANDSAT-8 Operational Land Instrument (hereafter S2 and L8) InfraRed images were used to track thermal anomalies of the Sabancaya dome (Fig. 2b). These multispectral sensors, launched in 2015/2017 (S2A&S2B) by ESA and in 2013 (L8) by NASA/USGS, offer a high spatial resolution in InfraRed wavelengths (20 and 30 m/pixel, respectively), with a total revisit time from a few days to weeks, depending on the latitude of the target. Here, we used the Top-Of-Atmosphere (TOA) reflectances for both data products, resampled for L8 data to 20 m spatial resolution, thus having the same geometric grid for both imageries. A recently automated hybrid hotspot detection algorithm, based on fixed ratios in SWIR bands (0.8 to 2.2 µm) and contextual statistical analysis (Massimetti et al., 2020), has been applied on Sabancaya S2-L8 dataset. The algorithm detects “hot” pixels in each image, from which we derive (i) the Hot Area exposed inside the crater, simply multiplying the number of hot pixels by their dimension (Fig. 3c); (ii) the evolution of the Thermal Profile (Fig. 3b), generated by stacking together the thermal

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 2** Examples of cloud-free satellite images used in this work to track the evolution of Sabancaya eruption; (a#) CNES/Airbus visible images⁽¹⁾ (through Google Earth Platform; b#) left side False Color RGB Sentinel 2 obtained with band combination 4–3–2 (visible bands at 10 m resolution) fused with infrared bands 12–11 (20 m resolution); right side vertical thermal profile obtained by summing the thermal index (TI) of hot pixels present in each row (Massimetti et al., 2020).</figcaption></figure>

The “hot” colormap is used in Fig. 3b to represent the intensity of the anomaly along with each stacked thermal profile. (c#) Radar image Sentinel 1 (radar geometry arbitrarily deformed for graphical convenience). Upper and lower panels display the typical configuration of the summit crater during the unrest (#1) and eruptive (#2) phases, respectively

profiles (oriented North–South) of each image (Fig. 2b) (see Laiolo et al., 2019; Massimetti et al., 2020).

## Volcanic Radiative Power (MODIS)

MIROVA (Middle InfraRed Observation of Volcanic Activity; www.mirovaweb.it) is a fully automatic hot-spot volcano detection system based on the analysis of Moderate Resolution Imaging Spectroradiometer (MODIS) data (Coppola et al., 2016). MODIS sensor has mounted onboard the Terra and Aqua NASA's satellites that were launched on February 2000 and May 2002 and covering their polar orbit with a temporal resolution acquiring about 4 images per day (2 nighttime and 2 daytime) at equator latitude with a spatial resolution of 1 km in the infrared bands. The MIROVA algorithm allows to detect, locate and quantify volcanic hotspots, measuring the Volcanic Radiative Power (VRP) that is heat flux radiated by hot volcanic features (> 200 °C Coppola et al., 2016). This approach provides VRP measurements with an error of ±30% (Wooster et al., 2003) and is used here to reconstruct the time series of VRP recorded at Sabancaya between 2012 and 2020. The same algorithm has been also applied to the VIIRS (Visible Infrared Imaging

Radiometer Suite) data, a multispectral sensor mounted onboard the Suomi-NPP (since January 2012), and on the NOAA-20 (since January 2018) satellites. In this specific case, we analyzed the night-time data acquired by the two infrared imaging bands of VVIRS, having a spatial resolution of 375 m (Campus et al., 2021). The better spatial resolution of VIIRS allows to detect thermal anomalies much lower than those detectable by MODIS (i.e. VRP < 1 MW), but suffer from saturation problems for more intense thermal anomalies (i.e. VRP > 10 MW per pixel). Overall, the combined analysis of MODIS and VIIRS data between January 2012 and December 2020 results in the detection of 2343 alerts with an average of 0.72 alerts per day (Fig. 3a).

## Ground-based dataset (INGEMMET)

### Ash plume height and color (camera network)

The network of optical cameras of the Sabancaya volcano was installed in 2013 and is composed of 4 optical cameras (SCOP, SAMP, SMUC, and SHUA; Fig. 1a) that operated in different periods. These stations are distributed around the volcano at distances between 4.3 and 28.7 km, and
Satellite-based

Ground-based

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Multiparametric plot of monitoring data collected during the pre-eruptive (unrest) and eruptive phases at Sabancaya, subdivided into a satellite- and ground-based parameters. (a) Volcanic Radiative Power (VRP) retrieved from MODIS (blue) and VIIRS (cyan) data; (b) Stacked thermal profiles (aligned N-S) along the summit area (See Fig. 2 b). Colors from red to bright yellow represent the increasing intensity of heat emission: Pre-eruptive crater rim is shown on the left; (c) “Total” lava dome area (from Planet Scope) and “Hot”</figcaption></figure>

lava dome area (from Sentinel 2); (d) Crater floor depth derived from Sentinel 1; (e) Height of the ash plume. Light gray and dark gray ash plumes are plotted using corresponding colors; (f) Daily SO₂ flux measured by the DOAS network; (g) Explosion rate (explosion per day) as recorded by the analysis of seismic data. The vertical dashed line corresponds to the beginning of the eruptive phase that occurred on 6 November 2016. See the text for a detailed explanation of each dataset

elevations comprised between 3590 and 5180 m (above sea level; asl). Volcanic emissions have been classified according to colors retrieved from the optical images, being yellowish, bluish, whitish, reddish, light gray, and dark gray. In this work, we used only the data relating to light or dark

gray plumes, which are related to the emission of volcanic ash plumes. The height of the ash plumes above the crater rim (H<sub>plume</sub>), is calculated by using a graphical scale that has been overlaid on the photographs of each camera, based on specific reference points (Fig. 1b). This approach provides

<figure><img src="image_9.png" /></figure>
height estimates with an error of ~10%. The time series of ash-plume height and color is presented in Fig. 3e.

## Sulfur dioxide flux (DOAS network)

Gas monitoring is performed by 3 DOAS scanners that are part of the Network for Observation of Volcanic and Atmospheric Change, NOVAC (Galle et al., 2010). This equipment is installed between 3 and 5 km from the Sabancaya crater (SAD1, SAD3, and SAD5; Fig. 1a) scans the sky looking for volcanic plumes passing overhead, and measures SO₂ flux (here provided in tons day⁻¹). The SO₂ flux measured in the 3 stations is evaluated according to the wind direction and the coverage of each piece of equipment (Ilanko et al., 2020). Only the valid readings of each station are considered and the one with the best coverage of the broadcast is used to provide the best estimate of daily SO₂ flux (Fig. 3f).

## Number and energy of explosions (seismic network)

The seismic network from Sabancaya volcano was installed gradually from 2009 until now. During the analyzed period, the network was composed of 18 broadband seismic stations (Fig. 1a) that operated in different periods. Thirteen stations were working with GURALP CMG – 6TD and the other five with SILICON sensor and MINIMUS digitizers. In the special case of seismic signals associated with explosions, the analysis, classification, and determination of its seismic energy, was performed using mainly the SB07 seismic data, because this station is situated only 2.6 km SE of the active vent and worked continuously during the analyzed period (Fig. 1a). The number of daily explosions is shown in Fig. 3g. To calculate the seismic energy, we used the equation proposed by Johnson and Aster (2005)

$$E_{\text{seismic}} = 2\pi r^2 c \frac{1}{A} \int S^2 U(t)^2 dt \quad (1)$$

where r is the source-station distance, ρ is the density (2600 kg m⁻³), c is the P wave velocity (3000 m s⁻¹), A is the attenuation correction, S is the seismic site response correction and U(t) is the particle velocity. For this calculation, we considered the data recorded by the closest seismic station (SB07), and we assumed the source was fixed below the active vent so that A and S were fixed at 1.

## Multiparametric characterization of the 2014–2020 activity of Sabancaya

The analyzed dataset (Fig. 3) is subdivided into two distinct main periods corresponding to (i) the unrest period (2014–6 November 2016) and (ii) the eruptive period (6 November

2016–31 December 2020). The key aspects characterizing these two phases are presented separately in the next sections.

## The unrest period (2012–2016)

The unrest of Sabancaya seems to have started many years before the November 2016 eruption, probably as early as 2003 when weak, but slowly increasing thermal anomalies were firstly detected by ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) inside the Sabancaya crater (Macedo et al., 2014; Jay et al., 2015; Reath et al., 2019b). In January 2012, thermal anomalies inside the Sabancaya crater started to be detected also by VIIRS (Fig. 3a), with growing VRP reaching about 0.3 MW by the end of the year. Clearer signs of reactivation began to manifest on 22–23 February 2013 when three earthquakes, with magnitudes between 4.6 and 5.2 ML, were recorded at shallow depth (~8 km) approximately 6 km northeast of the crater (Machacca et al., 2018). The swarm was followed by a gradual increase in fumarole activity, and SO₂ degassing from the summit crater (IGP-OVS, Annual Report 2014; Reath et al., 2019b; MacQueen et al., 2020) as well as by increasing VRP rising for the first time above 1 MW on 7 August 2014 (Fig. 4b). This increase was followed on 9 and 25 August 2014, with the first phreatic explosions of Sabancaya (IGP-OVS, Annual Report 2014). Notably, in the same period, MODIS detected its first thermal anomalies (Fig. 4a) to confirm the growing thermal flux observed by VIIRS. On 26 August 2014, just 1 day after the second phreatic explosion, also the Landsat 8 SWIR data detected the first hotspot inside the crater and suggests the enlargement of high-temperature zones (>200 °C) inside the crater (Massimetti et al., 2019) (Fig. 4c). Notably, the Landsat 8 and the following Sentinel 2 images allow us to appreciate that these precursory thermal anomalies were not associated with the development of punctual thermal anomalies (i.e. new hot cracks), but rather they were sourced by an almost uniform circular hot area (~150 m in diameter), located on the southern part of the crater (Fig. 2b and 4b).

After the phreatic explosions in August 2014, the heat flux continued to increase, until the second half of 2015. The SO₂ measurements acquired in 2014 and 2015, although sporadic, confirmed an increasing degassing trend (Fig. 4d) which was also corroborated by the satellite measurements retrieved from IASI (Infrared Atmospheric Sounding Interferometer; Moussallam et al., 2017) and OMI (Ozone Monitoring Instrument; Carn et al., 2017). In particular, OMI data indicated an annual average of SO₂ flux raising from ~230 tons day⁻¹ in 2014, to ~518 tons day⁻¹ in 2015 (Carn et al., 2017). Notably, the increasing degassing was characterized by periodic oscillations in the
<figure><img src="image_10.png" /><figcaption>**Fig. 4** Evolution of the thermal and gas emissions during the 2012–2016 unrest phase of Sabancaya. (a) Examples of Landsat 8 and Sentinel 2 images (band combination as in Fig. 2) show the location of thermal anomalies on the southern portion of the summit crater. Image timings are indicated by red circles in the plot (c). Note the appearance of a fumarolic zone on the northern crater rim in July and October 2016 images; (b) Time series of VRP (weekly average) obtained from the combined analysis of MODIS and VIIRS (see Fig. 3a); (c) “hot” area inside the crater as measured by Sentinel2 and Landsat8. (d) SO₂ flux is measured by the ground-based DOAS</figcaption></figure>

stations (this work), and other resources. Note that the first phreatic explosions that occurred in the summer of 2014 (black stars) were preceded by a relative maximum of VRP and roughly coincide with the first high-temperature thermal anomalies detected by LAND-SAT8/SENTINEL2. An increase in SO₂ degassing rate is observed since 2012 (annual averages for 2012 and 2015 are from Carn et al., 2017) but accelerates in April 2016. This acceleration is followed by the appearance of new fumaroles on the northern flank since July 2016 (horizontal arrows in b and c)

gas composition (periodicity of few minutes) which were interpreted by Moussallam et al. (2017) as cyclic input of deeper, gas-rich magma to shallow depth, likely due to vigorous magma convection within the conduit. Based

on the gas composition measured in November 2015, the same authors suggest a minimum equilibrium temperature of 730° C compatible with the presence of a dacitic magma composition at depth.

<figure><img src="image_11.png" /></figure>
From the second half of 2015, there has been a decrease in the VIIRS heat flux, coupled with a reduction in the hot area detected by Sentinel2 and Landsat 8. The decreasing trend reversed again in January 2016 when a new phase of increasing thermal flux is recorded (Fig. 4b). Since early 2016 ground-based SO₂ measurements have become more continuous and tracked this further acceleration of the unrest phase throughout the rest of the year. Notably, SO₂ emissions have tripled in just a few months reaching more than ~2000 tons day⁻¹ in mid-2016 (Fig. 4b). In May 2016, very high water vapor emissions were measured by Kern et al. (2017) and were interpreted as the first evidence of boiling off of the subsurface aquifer. During the same period, two new fumarolic areas appeared outside the crater (Kern et al., 2017; Global Volcanism Program, 2016), on the northern flank of Sabancaya, (Fig. 4a). The first fumarolic area, having an NW-SE trend, was located just outside the northern crater rim at about 5700 m a.s.l. (Fig. 4a). The second area was observed farther outside (~1 km) the crater, on the northeastern flank of the edifice at about 5600 m a.s.l. (Supplementary Figure S1 and Fig. 12 in Kern et al., 2017). These new fumarolic zones had temperatures comprised between 71 and 91 °C (Global Volcanism Program, 2016), not sufficient to be detected as “hotspots” by MODIS or VIIRS, neither by Sentinel2 (Fig. 4a, b). However, on the high-resolution Sentinel 2 images, it is still possible to appreciate the enlargement of these new degassing zones as evidenced by the bright areas in Fig. 3c and Supplementary Fig. S1. The northward expansion of these thermal anomalies was accompanied and followed by a further intensification of intra-crater degassing, with temperature anomaly and SO₂ flux reaching their respective maxima by the end of October 2016 (Fig. 4b and 4d). This very high thermal and gas emissions preluded the beginning of the eruptive phase that occurred on November 6.

## The eruptive period (2016–2020)

The eruptive phase began on 6 November 2016, and, at the time of this work, is still ongoing. We hereafter present the changes in the parameters monitored during this complex eruption that we subdivided into 6 phases. As a whole, this eruptive period was characterized by an initial vent-clearing stage, by three main episodes of dome growth (phase II, IV, and VI), one period of general stability (phase III), and one collapsing phase (V) (Fig. 5c). All these phases were accompanied by the occurrence of recurrent, discrete vulcanian explosions (VEI 1–2), with a general decreasing height (Fig. 5d). High thermal anomalies, fluctuating degassing and highly variable explosive seismicity accompanied the evolution of the lava dome and outlines the occurrence of complex eruptive dynamics (Fig. 5e,f,g).

## Phase I beginning of explosive activity (November-mid December 2016)

On 6 November 2016, a series of vulcanian explosions marked the beginning of the eruptive period. Phase 1 lasted a few weeks and consisted of frequent vulcanian explosions producing eruptive columns that gradually increased in height reaching up to 4.5 km above the crater by mid-December 2016 (Fig. 5d). On the other hand, the SO₂ flux showed a decreasing trend albeit remaining at high levels (> 2000 tons day⁻¹), comparable to the late stages of the pre-eruptive period (Fig. 4d and 5d). The number of explosions in this phase was the highest of the whole observed sequence (Fig. 5h), although the energies involved were still low (Fig. 5g). There was no evidence of the presence of lava domes in this purely explosive phase (Fig. 5a,b,c).

## Phase II Beginning of dome growth (December 2016–January 2018)

In mid-December 2016, thermal anomalies started to increase rapidly, suggesting the arrival of hot magma at the surface (Fig. 5a). However, evidence of a new lava dome appeared only in late December, when the Sentinel 1 satellite image of 26 December 2016 showed the rise of the crater floor of several tens of meters (Fig. 3d). This was followed on February 5, 2017, by the first Sentinel 2 image showing the thermal fingerprint of a new hot lava dome, of about 85 m in diameter, located in the northern sector of the crater, adjacent to the hot area detected during the unrest phase (Fig. 4a). Since then, the new lava dome, named “Huk” (“one” in the Quechua language; Ramos et al., 2019), began to grow inside the northern portion of the crater, reaching a volume of approximately 1.4 (±0.4)×10⁶ m³ by December 2017 (Fig. 5c) with an average extrusion rate of ±0.04 ± 0.02 m³ s⁻¹.

Phase II was characterized by evident thermal anomalies which accompanied the emplacement and gradual expansion of the Huk lava dome inside the crater (Fig. 5b). According to Sentinel 2 data, the heat was sourced by the whole dome surface, showing some fluctuation, eventually affected by a seasonal cloud covering (Fig. 5a,b). This initial extrusive phase was accompanied by intense explosive activity, characterized by evident variations in the number and energy of explosions (Fig. 5f, g), and by the emission of dark-gray ash plumes still reaching a height of 3–4 km above the crater. Since March 2017, the explosive activity slightly decreased and evolved into the emission of light-gray colored ash-plumes (Fig. 5g). This evolution was accompanied by a gentle decrease of the SO₂ flux which stabilized around 2000 tons day⁻¹ (Fig. 5e).

Notably, the overall dome growth observed during phase II was temporally interrupted by three deepening-refilling cycles
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>
◄ **Fig. 5** Multiparametric time series of the Sabancaya eruptive phase. (a) Volcanic Radiative Power (weekly average) obtained by the combined analysis of MODIS and VIIRS; (b) Thermal profiles (aligned N-S) along the summit area (see Fig. 3). The white dashed line envelops the maximum extension of the thermal anomaly (moving average over 1 month). The white arrows indicate the occurrence of pyroclastic flows outside the summit crater.; (c) Net lava dome volume derived by assuming an inverted truncated cone geometry (see the “Crater Depth and lava dome volume (SENTINEL 1)” section). Black arrows indicate the occurrence of deepening-refilling cycles “d-r”; Red cross indicates the estimated dome volume from a drone overflight on 26 October 2019. (d) SO₂ flux measured by the DOAS network; (e) Height of the ash plume above the crater. Different colors refer to the emission light gray and dark gray ash plume; (f) Explosion rate (red) and associated weekly seismic energy (black) released by explosions (g). The vertical dashed line corresponds to the beginning of the eruption on 6 November 2016

(hereby referred to as d-r cycles), characterized by the rapid lowering and raising of the crater floor (Fig. 3c). These cycles occurred on 15 February, 5 July, and 20 December 2017, and lasted 3, 7, and 4 weeks (Fig. 5c), respectively. Each episode was accompanied by a relative increase in the number of explosions (Fig. 5f) and was characterized by the deepening of the crater floor by 17–32 m. The deepening stages occurred at rates between 0.9 and 2.7 m d⁻¹ with pre-collapse levels recovered within 2 weeks (Fig. 3d).

## Phase III stable explosive activity, stable dome volume (January 2018–March 2019)

A marked reduction in the number of explosions occurred in January 2018 (Fig. 5). This low and stable explosive activity characterizes the entire phase III and was accompanied by apparent stability of the lava dome which stopped growing. However, the pause in the extrusive activity was still characterized by high and persistent thermal anomalies (Fig. 5a), likely sourced by the entire surface of the dome (Fig. 5b), eventually displacing toward the south (on 15 April 2018 an isolated thermal anomaly shows the presence of hot material in the southern part of the crater not occupied by the dome; Fig. 5b). The only notable event that occurred during the stable phase remains a sharp increase of the SO₂ emissions recorded in July 2018, when the flux exceeded 6000 tons day⁻¹, before returning to the typical levels (±1500 tons day⁻¹) that also characterized the previous phase (Fig. 5d). Surprisingly this degassing pulse did not perturb the explosive activity (which remained characterized by light-gray ash emissions) nor the apparent calm in the dome evolution, nor was it accompanied by a complimentary variation of the thermal flux.

## Phase IV acceleration in dome growth (March 2019–October 2019)

Phase IV started in March 2019 and corresponds to the reactivation and intensification of both the extrusive and the

explosive processes. Based on the measured parameters, this period has been divided into 4 distinct sub-phases. During phase IVa (13 March–11 May 2019) the “Huk” dome began to grow anew with a rate of 0.14 ±0.05 m³ s⁻¹. This recovery was characterized by a coeval increase in the SO₂ flux (Fig. 4e) as well as in the energy of the explosions (Fig. 5f, g). The increased explosive activity culminated on 11 May 2019 with the detection of hot material outside the crater, on the northern flank of the volcano (Fig. 5b and Supplementary Figure S2a). This was the first time that a thermal anomaly was detected outside the crater during the eruptive phase, and marks the transition into phase IVb (12 May–31 July 2019). This latter was characterized by a stable lava dome, low explosive activity, but highly irregular degassing. Two consecutive pulses in SO₂ flux (reaching 5000 tons day⁻¹ each; Fig. 5e) were recorded and accompanied the high thermal flux sourced by the dome surface (Fig. 5a,b). The pulsating degassing preluded a new episode of the dome’s growth which characterized the following phase IVc (1 August-25 September 2019). This phase of renewed extrusive activity occurred at a rate similar to phase IVa (0.15 ±0.05 m³ s⁻¹) and was followed by a new pulse of SO₂ (>6000 tons day⁻¹) and by explosive activity becoming gradually more frequent (Fig. 5e,f,g). Starting from late August 2019, the ash plumes became dark gray, similar to the first months of the eruption. On 16 September 2019, Landsat 8 recorded for the second time a thermal anomaly outside the crater, likely associated with the hot deposit extending

"500 m north of the crater (Fig. 5f and Supplementary Figure S2b). From this point onwards, begins phase IVd (25 September–23 October 2019) that sees the most evident acceleration of the extrusive and explosive activity (Fig. 5). In this phase, the dome grew at a rate of 0.75 ±0.25 m³ s⁻¹ and was accompanied by more than 200 explosive events per day. It is interesting to note that during this intense dome growth phase, the SO₂ flux remained moderate, with values below 2000 tons day⁻¹. On 26 October 2019, OVI-INGEM-MET conducted a drone overflight that allowed to measure a lava dome volume of approximately 4.6×10⁶ m³ (red cross in Fig. 5c). By this time, Sentinel1-derived measurements indicated a volume of 4.3. (±1.4)×10⁶ m³ which represents the maximum volume reached by the “Huk” lava dome.

## Phase V slow dome/crater collapse (November 2019–September 2020)

Starting from November 2019, there was a clear reversal of the trend that characterizes the previous phase. The number and energy of the explosions dropped rapidly to settle on particularly low values never recorded since the eruption onset (Fig. 5g,f). The SO₂ flux and plume height remained low (Fig. 5e), and were accompanied by the slow subsidence of the dome surface which deepened at an average velocity
of 0.22 m d⁻¹ (about an order of magnitude less than the d-r cycles). During this sinking phase, the thermal anomalies dropped to a minimum (<10 MW; Fig. 5a) and were spatially distributed along a ring structure surrounding the crater center (Fig. 5b and Supplementary Figure S3). The slow subsidence was temporally interrupted by a new d-r cycle characterized by a rapid deepening of "50 m that occurred at a speed of "3.4 m s⁻¹. Also during this waning phase, the occurrence of the d-r cycle was coeval with a surge of thermal and explosive activity (Fig. 5a,f,g).

## Phase VI (September 2020–December 2020) new phase of dome growth

A new lava dome, named “Iskay” (“two” in the Quechua language; Ortega et al., 2021) started to grow in the north-eastern sector of the collapsed crater in September 2020 (Fig. 5 and Supplementary Fig. S3 and S4). Unfortunately, this new dome is poorly visible by S1 images (Fig. 5c). The reason for this is that the dome grew very close to the northern rim of the crater, and the row where the depth is calculated does not intersect this new dome. However, the new extrusive phase was accompanied by the sharp resumption of thermal activity with VRP > 25 MW and the thermal anomalies slightly shifted to the north where the new dome was growing (Fig. 5b). At the beginning of this new phase, the SO₂ flux increased threefold, although after reaching a peak in late August 2020 ("5600 tons day⁻¹) the degassing trend started to decline. The energy and number of explosions recorded during Phase VI indicate that the extrusion

of the new dome was accompanied by renewed explosive activity and by the emission of dark-gray ash plumes.

## Shallow magma budget

The multiparametric dataset has been used to quantify the shallow magma budget (Harris and Stevenson, 1997a, b) by estimating magma volumes involved in the different processes responsible for the observed activity (Fig. 6a, b). More specifically, we calculate (i) the volume of the extruded lava dome ($V_{dome}$), (ii) the volume of emitted ash ($V_{ash}$), (iii) the volume of degassed magma ($V_{degas}$), (iv) the volume of magma radiating an equivalent amount of thermal energy ($V_{therm}$). Below, we define these 4 types of magma volumes in more detail, by specifying the assumptions and the method used to calculate each of them. All the volumes were recalculated for Dense Rock Equivalent (DRE) value by assuming a DRE magma density ($ρ_{DRE}$) of ~2700 kg m⁻³ (Aguilar et al., 2019).

(i) *Volume of extruded magma ($V_{dome}$)* The bulk lava dome volume is derived from the lava dome area and crater floor depth, as explained in the “Crater depth and lava dome volume (SENTINEL 1)” section. Previous studies on lava domes have found extremely variable lava dome densities as a function of magma composition and voids fraction (Harnett et al., 2019). For example, muon tomography studies over lava domes, show densities as low as ~1000 kg m⁻³ (Tanaka et al., 2009), which applied to the dome of the Sabancaya would produce a DRE volume equal to 0.5 × 10⁶ m³. Unfortunately, the density of the

**Fig. 6** Magma budget of the Sabancaya eruption (all data smoothed over 3 months for clarity). (a) Cumulative volumes of extrusive ($V_{dome}$) and explosive activity ($V_{ash}$). (b) Cumulative volumes of degassed magma ($V_{degas}$) and thermally radiant magma ($V_{therm}$). Note the strong excess of these two volumes to those related to the extrusive ($V_{dome}$) and explosive($V_{ash}$) activity alone. (c) Timeseries of the relevant magma fluxes were obtained by the magma budget calculations. The excess degassing and excess thermal radiation appear as magma fluxes ($Q_{degas}$, $Q_{therm}$) in clear excess to erupted magma fluxes ($Q_{dome}$, $Q_{ash}$). Note that the curve of $Q_{dome}$ results fragmented during period of dome stability or collapse ($Q_{dome}=0$)

<figure><img src="image_14.png" /></figure>
Sabancaya dome is unknown. For this reason, in this work, we assumed a range of densities spanning from 1500 to 2500 kg m⁻³ to take into account variable emplacement conditions (Fig. 6a).

(ii) Volume of ash ($V_{ash}$) to estimate the volume of ash emitted we follow the approach of Londono et al., (2018) which applied the Druitt et al. (2002) and Bonadonna et al. (2002) models to vulcanian explosions of Nevado del Ruiz volcano. Accordingly, the height of an ash plume ($H_{plume}$) above the crater is related to the mass of ash ($M_{ash}$) by

$$H_{plume} = BM_{ash}^{0.25} \qquad (2)$$

where B is a model-derived coefficient (in m kg⁻⁴) that takes into account the specific thermal conditions of vulcanian plumes ascent (Londono et al., 2018). Here, we inverted Eq. (2) to estimate $M_{ash}$ from the measured $H_{plume}$ and associated errors (see the “Ash plume height and color (camera network)” section). We used $B = 55$ m kg⁻⁴ which has been applied to model volcanic plumes at Nevado del Ruiz volcanoes, having magma compositions and eruptive styles similar to Sabancaya (Londono et al., 2018). The mass of ash is retrieved for each explosion detected by the camera network (Fig. 1b), then cumulated over weekly time windows and converted into DRE volumes.

(iii) Volume of degassed magma ($V_{degas}$) this is the magma volume that has degassed the measured mass of SO₂ (MSO₂) throughout the analyzed period. This is calculated by using the petrological method (Shinohara 2008) according to which

$$V_{degas} = MSO_2/(2 * (\rho_{DRE} * \Delta X_S)) \qquad (3)$$

where MSO₂ is the mass (in kg) of degassed SO₂ and ΔX_S is the S volatile loss. Having no direct estimates on the initial content of S, we assumed 2 distinct scenarios. In the first case (Scenario 1) we assume ΔX_S equal to 500 ppm, here considered as an upper boundary for andesitic compositions (Shinohara 2008). In the second case (Scenario 2), we assume that gas emissions were supplied by an S-rich mafic melt (basalt to basaltic andesite) mixing with the andesitic melt in the shallow plumb-ing system of Sabancaya (Gerbe & Thouret, 2004). This model was recently proposed to explain the sustained SO₂ emissions at the andesitic Mt. Cleveland volcano, where olivine-hosted melt inclusions in andesitic magmas were found to contain up to 2300 ppm S (Werner, et al., 2017, 2020). For this scenario, we thus assume ΔX_S equal to 2500 ppm to represent the initial concentrations in the S-rich mafic melt. However, we must point out that the subduction setting is extremely different between the Aleu-tians and South America, so that this high sulphur content is to be considered as an extreme end-member scenario, not based on any petrological evidence.

(iv) Volume of radiating magma ($V_{thermal}$) this is the volume of magma retrieved from the thermal approach, and related to the measured radiant energy. According to Coppola et al., 2013, during an effusive eruption, this volume is related to the volcanic radiant energy (VRE) through

$$V_{thermal} = \frac{VRE}{crad} \qquad (4)$$

where $c_{rad}$ is an empirical coefficient that takes into account the effective rheology of the emplacing lava body. For the Sabancaya lava dome, we used an average silica content of 62.5 wt.% (Ortega et al., 2021) and we calculate (following Coppola et al., 2013) a $c_{rad}$ equal to 1.2 (±0.6) × 10⁷ J m⁻³, which is typical for an andesitic lava dome. It is important to note that the error of $c_{rad}$ (±50%) expresses the uncertainty in using the silica content as a simple proxy of the rheological and insulation conditions of the erupted lava (Coppola et al., 2013). This results in an upper and lower limit of $V_{thermal}$ which likely includes the real emplacement conditions for the Sabancaya lava dome.

All volumes (upper and lower estimates) were calculated on a 7-day basis (Fig. 6) to provide the weekly fluxes and associated errors ($Q_{dome}$, $Q_{ash}$, $Q_{degas}$, and $Q_{thermal}$).

## Excess degassing and excess thermal radiation

Magma budget calculations (Fig. 6) indicate that about 0.002 km³ of magma were extruded to form the “Huk” and “Iskay” lava domes, while 0.04–0.09 km³ of ash were produced by the explosive activity according to Eq. (1) (Fig. 6a). Although it is unclear whether a non-juvenile, previously degassed magma may have contributed to the total erupted volume to some degree, it is evident that the total erupted magma volume ($V_{erupt} = V_{dome} + V_{ash}$) was 1 or 2 order of magnitude lower than the degassed magma ($V_{degas} ± 0.25–1.28$ km³; Fig. 6b) depending on the assumed degassing scenario. This unbalance ($V_{degas} > V_{erupt}$) is typical of open vent volcanoes and suggests that not only the unrest phase of Sabancaya was characterized by degassing of unerupted magma (Moussallam et al., 2017; Jay et al., 2015; MacQueen et al., 2020), but also that the following eruptive activity was characterized by a strong excess of magma degassing (Fig. 6c).

Similarly, the “radiant magma volume”, derived from the thermal approach ($V_{thermal} = 0.04–0.13$ km³) is higher (1 order of magnitude) than the erupted magma volume (Fig. 6b) suggesting that most of the detected heat (VRP) was sourced by additional process(es) other than the simple magma extrusion (Harris and Stevenson 1997a, b). By analogy with gas emissions, we thus introduce the notion of
"excess thermal radiation", to describe the thermal unbalance that occurs when there is more magma releasing heat (and cooling down) than the magma erupted at the surface ($V_{thermal} > V_{erupt}$).

Two models have been invoked to explain the excess degassing in subduction volcanoes (Shinohara, 2008). The first model (permeability controlled flow model; Edmonds et al., 2003) assumes the gas flow derives from mafic magma at depth and ascends along a permeable conduit to the surface. In this case, bubble-magma separation occurs at depth (i.e. within a convective magma chamber) with gas emissions controlled by the rate of magma supply and surface thermal emissions also affected by the length of the permeable conduit (Stevenson, 1993). The second model (convection controlled flow model; Shinohara, 2008) assumes the gas flow is continuously supplied by a convective magma column that rises at shallower levels, within the volcanic conduit (Stevenson and Blake, 1998). Eventually, when the top of the convecting magma reaches the surface, a lava dome is extruded and the gas and heat

flow becomes mainly controlled (i) by the convective rate inside the conduit and (ii) by the permeability of the lava. From a thermal point of view, these two models differ only in the depth at which the convection occurs (magma chamber or volcanic conduit), and therefore in the length of the permeable plug that the gas must pass through (Stevenson, 1993).

At Sabancaya, these two models seem to represent two end members of a continuous process (Fig. 7), in which convection occurs at decreasing depths, starting from a crustal magma chamber, eventually fed by mafic magma (model 1) up to shallow levels in the conduit where andesitic magma reaches the base of the extruding lava dome (model 2).

The outgassing and the low thermal anomalies recorded during the unrest phase were effectively sourced by the southern part of the crater floor, presumably occupied by the remains of the buried conduit feeding the 1988–1997 eruption (Fig. 4a). This structure was still permeable to the flow of gas and heat released by the rising magma and may have a role in attenuating the pressurization of the shallow magma

<figure><img src="image_15.png" /><figcaption>**Fig. 7** Conceptual models of excess degassing and excess thermal radiation of Sabancaya volcano during unrest and eruption (adapted from Shinohara, 2008). (a) During the initial phase of unrest degassing of unerupted mafic magma occurred through a permeable fluid flow along the old conduit of Sabancaya. The rise of magma into the crust (b) caused the growth of the convective degassing cell and the decrease of the height of the permeable plug. The two processes promoted increasing gas and heat emissions recorded at the surface.</figcaption></figure>

At more shallow levels (c), the resistance of the old plug causes the ascending magma to migrate northwards and interact with the hydrothermal system and the opening of new fumaroles. (b) During the eruption a convective magma column feeds the dome's activity at the surface, causing a large quantity of magma to degass and cool at very shallow levels. This produces high thermal anomalies both sourced by the extrusive/explosive activity at the surface, as well as by the very shallow magma circulation below the lava dome
chamber (Girona et al., 2014). In this scenario, the rise of the magma column translates into an ever-larger SO₂ degassing cell (i.e. the portion of the convective magma column above the SO₂ exolution level) as well as into a decreasing height of the permeable conduit that the gas must pass through, losing heat (Fig. 7). We thus suggest that the ascent of magma in the unrest phase of the Sabancaya caused the increasing thermal- and gas-derived magma flux characterizing (Fig. 4 and Fig. 6c) since the convective magma cell becomes bigger and closer to the surface (Fig. 7a,b).

Nonetheless, the appearance of new fumaroles on the north flank during the period of increasing degassing (June 2016) can be interpreted as a signal that an ascending magmatic intrusion had intercepted the hydrothermal system (Kern et al., 2017) and progressively moved toward the north to a position adjacent to the old permeable plug (Fig. 4a). This transition suggests that while the upper part of the old plug was found to be permeable to gas flow, it still formed a physical barrier to the magma ascent which had to open a new pathway further north (Fig. 7c).

Our calculations (Fig. 6c) suggest that once the eruption has started the magma flux required to sustain the measured gas emissions ($Q_{\text{degas}}$) ranges from 2 to $10 \text{ m}^3 \text{ s}^{-1}$ (depending on the degassing scenario considered). This value is compared with an average magma output rate ($Q_{\text{out}} = Q_{\text{dome}} + Q_{\text{ash}}$) lower than $0.08 \text{ m}^3 \text{ s}^{-1}$. Similarly, the volumetric flux needed to explain the thermal emissions ($Q_{\text{thermal}}$) ranges from $0.3$ to $1 \text{ m}^3 \text{ s}^{-1}$ (Fig. 6c), which is still 4 to 10 times $Q_{\text{out}}$ and again suggests the presence of unerupted magma degassing and cooling at a very shallow level.

During the eruptive activity, the excess degassing and thermal radiation were both sourced by the dome surface (Fig. 2, 3) which was continuously disrupted by the explosive activity. A portion of the excess thermal radiation was probably associated with the explosive activity, which repeatedly exposed the hot core to the atmosphere. This is supported by the good correlation between the VRP and the surface explosive activity (Fig. 5a and 5g) and, indirectly, from the temporally consistent fluctuation between $Q_{\text{thermal}}$ and $Q_{\text{ash}}$ (Fig. 6). However, to keep the thermal emission high, after repeated explosions of the lava dome, a continuous supply of heat is necessary, which again can be explained by a convective magmatic column, located under a highly permeable dome, and modulating the activity at the surface.

A further argument in favor of shallow magmatic convection during the eruption is provided by the subsidence/collapse periods of the lava dome. Two different types of dome subsidence/collapse were observed during the whole eruptive period. The first type corresponds to the sudden deepening-refilling cycles (d-r cycles) observed during phases II and V. These events lasted for a few weeks at most (rapidly

recovering their pre-deepening conditions) and were accompanied by a relative increase in explosive activity (number of events per day). However, they were not associated with any other variations in the magma feeding system (degassing, thermal, explosion energy; Fig. 5). We thus suggest that these events tracked the occurrence of temporary, local, and/or extremely superficial disrupting dynamics, involving exclusively a crater sector or specific portions of the lava dome. On the contrary, the subsidence recorded during phase V was slower, involved the entire crater sector, and was accompanied by the lowering of explosive and thermal activity to their minimum values. This phase of subsidence occurred just after the extrusive pulse of phase IV which led the Huk lava dome to reach its maximum volume and height (Fig. 3 and 5). Once the extrusive pulse of phase IV was exhausted, the magma column was no longer able to maintain itself at the same level (pressure) and probably began to collapse causing the gradual subsidence of the dome (see Supplementary Figure S4). Therefore, rather than being associated with superficial processes, as the d-r cycles, this slow subsidence seems to be linked to the lowering of the shallow magma column that sustained the active lava dome. During the collapsing period, we observed the formation of annular thermal anomalies (see Supplementary Figure S3) compatible with degassing hot rings surrounding the crater's center. Similar 'drain back' mechanisms and annular thermal features have been observed at other andesite domes (Oppenheimer et al. 1993; Matthews et al., 1997; Werner et al., 2017) and support the inference of a low viscosity convective magma column at very shallow levels, whose collapse can be determined by the loss of pressure inside the conduit (Matthews et al., 1997).

All these arguments support the idea that a very shallow magmatic convection is at the origin of both excess degassing and excess thermal radiation. It is worth noting that the minimum limit of $Q_{\text{degas}}$ (mafic magma supply scenario) and the maximum limit of $Q_{\text{thermal}}$ almost converge for a rate of magma convection equal to $1-2 \text{ m}^3 \text{ s}^{-1}$ (Fig. 6c). This value is consistent with the continuous deformation measured under the volcano, whose inversion suggests sub-crustal volume's change in the order of $0.6-1.5 \text{ m}^3 \text{ s}^{-1}$ (Boixart et al., 2020; MacQueen et al., 2020).

On the other hand, it seems if we take into consideration the full variability of values of both $Q_{\text{thermal}}$ and $Q_{\text{degas}}$, they provide different convective rates. These differences may be associated with limitations for the heat emissions to reach the surface (i.e. part of the heat is going to be released through the conduit walls and will never reach the surface). However, we do not rule out that the fluxes obtained from gas and thermal data might correspond to the rate of magma transport (convection) in different portions (cells) of the magmatic system. While $Q_{\text{degas}}$ corresponds to the rate of magma circulation above
the exsolution level of SO₂ (i.e. within the full convective degassing cell), the thermally-derived magma flux (Q<sub>thermal</sub>) would correspond to a smaller, shallower portion of the degassing system. This zone includes the permeable lava dome and the upper portion of the convective magma column, whose persistent heat loss is compensated by the continuous supply of magma from depth. Whether and to what extent these fluxes are linked to the efficiency of magma transport at various levels of the plumbing system is an intriguing topic that deserves further study. The continuous acquisition of multi-parameter data (satellite- and ground-based) at Sabancaya and other dome-forming, open-vent volcanoes will possibly clarify these aspects.

## Conclusions

We used a multiparametric dataset of satellite- (visible, infrared, and radar images) and ground-based (optical camera, SO₂ flux, explosion rate, and energy) measurements to characterize the activity of Sabancaya between 2012 and 2020. Thermal, degassing, extrusive and explosive activity has been quantified and allowed us to subdivide the analyzed period into a phase of unrest (2012–2016) and a dome-forming eruptive phase (2016–2020); with the latter still ongoing at the time of writing (September 2021).

During the unrest phase, the rise of a convective magma column caused increasing SO₂ degassing (from ~200 to ~2000 tons day⁻¹) and thermal emission (from ~0.04 to ~1.5 MW) at the surface, favored by the presence of an old conduit permeable to gas flow. Conversely, the complex eruptive phase that followed was characterized by episodes of dome growth (at rates from 0.04 to 0.75 m³ s⁻¹) and collapse, accompanied by variable heat flux (0.5–25 MW), persistent but irregular degassing (700–7000 tons day⁻¹), and weak to moderate seismic energy of explosions associated with repeated vulcanian ash plumes (500–5000 m above the crater).

Overall, the eruptive activity was characterized by a strong excess of gas and thermal radiation both indicative of an efficient magmatic convection process, operating up to the shallowest levels of the column and causing the extrusive and explosive activity of the overlaying lava dome. Based on magma budget calculations, we estimated that between 2016 and 2020, at least 0.25–1.28 km³ of magma degassed and at shallow depth, before being cycled back to the deeper magma chamber. A fraction of this degassed magma volume (0.04–0.13 k m³) reached the top of the convective magmatic column, losing heat and producing an excess amount of heat compared to that simply radiated by the dome. Finally, less than 0.001 km³ have erupted through a combination of extrusive and explosive

activity. Our results indicate that, like many other dome-forming eruptions, the eruptive activity of Sabancaya is characterized by an evident excess of gas and thermal release. We suggest that these two types of excesses (degassing and radiation) are indicative of magma convection (possibly at different levels of the plumbing system) and play a fundamental role in determining the mid-to-long-term eruptive dynamics that are observed on the surface.

## Supplementary information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00445-022-01523-1.

**Acknowledgements** We thank the Editor, Christoph Kern, and an anonymous reviewer for the comments made on the first version of the manuscript and for the constructive suggestions they gave to the interpretation of the data. MODIS and VIIRS data were provided by the LANCE-MODIS system (https://lance.modaps.eosdis.nasa.gov/) and LAADS-DAAC system (https://ladsweb.modaps.eosdis.nasa.gov/), respectively. Sentinel data (1 and 2) were provided by Copernicus Open Access Hub (https://scihub.copernicus.eu/). Planet Scope and CNES / Airbus images are distributed by © 2021 Planet Labs Inc and © 2021 Google Earth platforms, respectively.

**Funding** This work was partially funded by the Italian Ministry of University and Research (MUR). Part of this research was funded thanks to the PAPIIT project IA102221. https://ladsweb.modaps.eosdis.nasa.gov/

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Aguilar, R.C., 2019. Exposure-based risk assessment of tephra-fall associated with long-lasting Vulcanian eruptions at Sabancaya volcano, Peru. These de Specialisation in the assessment and management of geological and climate related risk - CERG-C 2018, University of Geneva.

Aldeghi C, Escobar-Wolf G (2019) Volcano monitoring from space USING High-cadence PLANET Cubesat IMAGES applied to Fuego Volcano. Guatemala Remote Sensing 11:2151. https://doi.org/10.3390/rs11182151

Andres RJ, Rose WI, Kyle PR, DeSilva S, Francis P, Gardeweg M, Moreno Roa H (1991) Excessive sulfur dioxide emissions from Chilean volcanoes. J Volcanol Geoth Res 46:323–329. https://doi.org/10.1016/0377-0273(91)90091-d
Boixart G, Cruz LF, Miranda Cruz R, Euillades PA, Euillades LD, Battaglia M (2020) Source model for Sabancaya volcano constrained by dinsar and GNSS surface Deformation Observation. Remote Sensing 12:1852. https://doi.org/10.3390/rs12111852

Bonadonna C, Macedonio G, Sparks RS (2002) Numerical modelling of TEPHRA fallout associated with dome collapses and Volcanian explosions application to hazard assessment on Montserrat. Geological Society, London, Memoirs 21:517–537. https://doi.org/10.1144/gsl.mem.2002.021.01.23

Boudon, G., Balcone-Boissard, H., Villemant, B., Morgan, D.J., 2015. What factors control superficial lava dome explosivity? Sci Rep 5. https://doi.org/10.1038/srep14551

Calder ES, Luckett R, Sparks RS, Voight B (2002) Mechanisms of lava dome instability and generation of rockfalls and pyroclastic flows at Soufrière Hills Volcano, Montserrat. Geological Society, London, Memoirs 21:173–190. https://doi.org/10.1144/gsl.mem.2002.021.01.08

Calder, E.S., Lavallée, Y., Kendrick, J.E., Bernstein, M., 2015. Lava dome eruptions. The Encyclopedia of Volcanoes 343–362. https://doi.org/10.1016/B978-0-12-385938-9.00018-3

Campus A., Coppola D., (2021). Monitoring volcanic thermal emission using VIIRS. Abstract Volume, 4a Conferenza "A. Rittmann" Giovani Ricercatori - Online Conference, April 6th-9th 2021. C. Montagna, E., Nicotra and G.Risica eds, Misc. INGV 59 1–88. https://doi.org/10.13127/misc/59

Carn SA, Watts RB, Thompson G, Norton GE (2004) Anatomy of a lava DOME collapse the 20 March 2000 event at Soufrière Hills Volcano, Montserrat. J Volcanol Geoth Res 131:241–264. https://doi.org/10.1016/S0377-0273(03)00364-0

Carn, S.A., Fioletov, V.E., McLinden, C.A., Li, C., Krotkov, N.A., 2017. A decade of global volcanic SO2 emissions measured from space. Sci Rep 7. https://doi.org/10.1038/srep44095

Cigna F, Tapete D, Lu Z (2020) Remote sensing of volcanic processes and risk. Remote Sensing 12:2567. https://doi.org/10.3390/rs12162567

Coppola D, Laiolo M, Piscopo D, Cigolini C (2013) Rheological control on the radiant density of active lava flows and domes. J Volcanol Geoth Res 249:39–48. https://doi.org/10.1016/j.jvolgeores.2012.09.005

Coppola D, Laiolo M, Cigolini C, Donne DD, Ripepe M (2015) Enhanced volcanic hot-spot detection using Modis IR data results from the Mirova system. Geological Society, London, Special Publications 426:181–205. https://doi.org/10.1144/sp426.5

Coppola, D., Laiolo, M., Cigolini, C., Massimetti, F., Delle Donne, D., Ripepe, M., Arias, H., Barsotti, S., Parra, C.B., Centeno, R.G., Cevuard, S., Chigna, G., Chun, C., Garaebiti, E., Gonzales, D., Griswold, J., Juarez, J., Lara, L.E., López, C.M., Macedo, O., Mahinda, C., Ogburn, S., Prambada, O., Ramon, P., Ramos, D., Peltier, A., Saunders, S., de Zeeuw-van Dalfsen, E., Varley, N., William, R., 2020. Thermal remote sensing for global volcano monitoring experiences from the Mirova system. Front Earth Sci 7. https://doi.org/10.3389/feart.2019.00362

Costa, A., Melnik, O., Sparks, R.S., Voight, B., 2007. Control of magma flow in dykes on cyclic lava dome extrusion. Geophys Res Lett 34. https://doi.org/10.1029/2006GL027466

Druitt TH, Young SR, Baptie B, Bonadonna C, Calder ES, Clarke AB, Cole PD, Harford CL, Herd RA, Luckett R, Ryan G, Voight B (2002) Episodes of cyclic VULCANIAN explosive activity with fountain collapse at Soufrière Hills Volcano, Montserrat. Geological Society, London, Memoirs 21:281–306. https://doi.org/10.1144/gsl.mem.2002.021.01.13

Edmonds M, Oppenheimer C, Pyle DM, Herd RA, Thompson G (2003) SO2 emissions from Soufrière Hills volcano and their relationship to conduit permeability, hydrothermal interaction and degassing regime. J Volcanol Geotherm Res 124:23–43. https://doi.org/10.1016/S0377-0273(03)00041-6

Eichelberger JC, Carrigan CR, Westrich HR, Price RH (1986) Non-explosive silicic volcanism. Nature 323:598–602. https://doi.org/10.1038/323598a0

Elsworth D, Voight B, Thompson G, Young SR (2004) Thermal-hydrologic mechanism for rainfall-triggered collapse of lava domes. Geology 32:969. https://doi.org/10.1130/g20730.1

Furtney MA, Pritchard ME, Biggs J, Carn SA, Ebmeier SK, Jay JA, McCormick Kilbride BT, Reath KA (2018) Synthesizing multi-sensor, multi-satellite, multi-decadal datasets for global volcano monitoring. J Volcanol Geoth Res 365:38–56. https://doi.org/10.1016/j.jvolgeores.2018.10.002

Galle, B., Johansson, M., Rivera, C., Zhang, Y., Kihlman, M., Kern, C., Lehmann, T., Platt, U., Arellano, S., Hidalgo, S, 2010. Network for Observation of Volcanic and Atmospheric Change (NOVAC)—a global network for volcanic gas monitoring network layout and instrument description. Journal Geophysical Research 115D05304. https://doi.org/10.1029/2009JD011823

Gerbe, M.-C., Thouret, J.-C., 2004. Role of magma mixing in The Petrogenesis of tephra erupted during The 1990–98 explosive activity of Nevado SABANCAYA, southern Peru [WWW Document]. Bulletin of Volc

Girona T, Costa F, Newhall C, and B. Taisne, (2014) On depressurization of volcanic magma reservoirs by passive degassing. J Geophys Res Solid Earth 119:8667–8687. https://doi.org/10.1002/2014JB011368

Harnett CE, Kendrick JE, Lamur A, Thomas ME, Stinton A, Wallace PA, Utley JEP, Murphy W, Neuberg J and Lavallée Y (2019) Evolution of mechanical properties of lava dome rocks across the 1995–2010 eruption of Soufrière Hills Volcano, Montserrat. Front. Earth Sci. 77. https://doi.org/10.3389/feart.2019.00007

Harris AJL, Stevenson D (1997a) Magma budget and steady-state activity of Volcano and Stromboli. Geophys Res Lett 24(9):1043–1046. https://doi.org/10.1029/97GL00861

Harris AJL, Stevenson D (1997b) Thermal observations of degassing open conduits and fumaroles at Stromboli and Volcano using remotely sensed data. J Volcanol Geotherm Res 76(3-4):175–198. https://doi.org/10.1016/S0377-0273(96)00097-2

IGP-OVS, 2014. Annual Report, Sabancaya https://www.igp.gob.pe/servicios/centro-vulcanologico-nacional/sites/default/files/2020/reporte_sabancaya_201407.pdf

Ilanko, T., Pering, T.D., Wilkes, T.C., Apaza, F., Kern, C., Moreno, A., De Angelis, S., Layana, S., Rojas, F., Vasconez, F., McGonigle, A.J.S, 2019. Degassing at Sabancaya volcano measured by UV cameras and the NOVAC network. Volcanica 2, 239–252. https://doi.org/10.30909/vol.02.02.239252

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

Jay JA, Delgado FJ, Torres JL, Pritchard ME, 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 JB, Aster RC (2005) Relative partitioning of acoustic and seismic energy during Strombolian eruptions. J Volcanol Geoth Res 148:334–354. https://doi.org/10.1016/j.jvolgeores.2005.05.002

Juvigné E, Thouret JCL, Gilot E, Leclercq L, Gourgaud A (1998) L'activité du volcan Nevado Sabancaya (Pérou) Au Cours DE L'HOLOCÈNE [Activity of Nevado Sabancaya Volcano (peru) throughout the holocene]. Quaternaire 9:45–51. https://doi.org/10.3406/quate.1998.2105

Juvigné, E., Thouret, J.C., Loutsch, I., Lamadon, S., Frechen, M., Fontugne, M., Rivera, M., Dávila, J., Mariño, J., 2008. Retombées volcaniques Dans DES tourbières Et lacs autour Du massif des Nevados ampato ET Sabancaya (PÉROU méridional, Andes
Centrales). Quaternaire 157–173. https://doi.org/10.4000/quate
rnaire.3362

Kern C, Masias P, Apaza F, Reath K, Platt U (2017) Remote measure-
ment of high preeruptive water vapor emissions at Sabancaya vol-
cano by passive differential optical absorption spectroscopy. Jour-
nal Geophysical Research Solid Earth 122:3540–3564. https://doi.
org/10.1002/2017JB014020

Laiolo M, Ripepe M, Cigolini C, Coppola D, Della Schiava M, Genco R, Innocenti L, Lacanna G, Marchetti E, Massimetti F, Silengo MC (2019) Space- and ground-based geophysical data tracking of magma migration in shallow feeding system of Mount Etna Volcano. Remote Sensing 11:1182. https://doi.org/10.3390/rs11101182

Londono, J.M., Galvis, B., 2018. Seismic data, photographic images and physical modeling of volcanic plumes as a tool for monitoring the activity of Nevado del Ruiz volcano, Colombia. Frontiers in Earth Science 6. https://doi.org/10.3389/feart.2018.00162

Macedo, O. (2018). Observaciones de campo en la zona de Pinchollo al norte del volcán Sabancaya, Informe Técnico, Instituto Geofísico del Perú.

Machacca Puma, R., Ortega, M. A., Gonzáles, R., Miranda Cruz, R., Apaza Choquehuayta, F. E., Lazarte Zerpa, I. A., et al. (2018). Monitoreo Multi-paramétrico del volcán Sabancaya y evolución del proceso eruptivo 2016–2017. In Libro de resúmenes VIII Foro Internacional los Volcanes y su Impacto, Arequipa, 26 y 27 de abril del 2018, pp. 105–108. Arequipa, Peru, Instituto Geológico Minero y Metalúrgico -INGEMMET.

MacQueen, P., Delgado, F., Reath, K.A., Pritchard, M.E., Bagnardi, M., Milillo, P., Lundgren, P.R., Macedo, O., Aguilar, V., Ortega, M., Anccasi, R., Lazarte Zerpa, I.A., Miranda, R., 2020. Volcano-
tectonic interactions at Sabancaya volcano, peru eruptions, mag-
matic inflation, moderate earthquakes, and fault creep. https://doi.
org/10.1002/essoar.10501561.1

Massimetti F, Coppola D, Laiolo M, Valade S, Cigolini C, Ripepe M (2020) Volcanic hot-spot detection Using Sentinel-2 A comparison with MODIS-MIROVA thermal data series. Remote Sensing 12:820. https://doi.org/10.3390/rs12050820

Matthews SJ, Gardeweg MC, Sparks RS (1997) The 1984 to 1996 cyclic activity of Lascar Volcano, northern Chile cycles of dome growth, dome subsidence, degassing and explosive eruptions. Bull Volcanol 59:72–82. https://doi.org/10.1007/s004450050176

Melnik, O., Sparks, R.S., 2005. Controls on conduit magma flow dynamics during lava dome building eruptions. Journal of Geophysical Research Solid Earth 110. https://doi.org/10.1029/2004j b003183

Melnik O, Sparks RS (1999) Nonlinear dynamics of lava dome extru-
sion. Nature 402:37–41. https://doi.org/10.1038/46950

Moussallam Y, Tamburello G, Peters N, Apaza F, Schipper CI, Curtis A, Aiuppa A, Masias P, Boichu M, Bauduin S, Barnie T, Bani P, Giudice G, Moussallam M (2017) Volcanic gas emissions and degassing dynamics at Ubinas and Sabancaya volcanoes; implications for the volatile budget of the central volcanic zone. J Volcanol Geoth Res 343:181–191. https://doi.org/10.1016/j.jvolg eores.2017.06.027

Moussallam, Y., Barnie, T., Amigo, A., Kelfoun, K., Flores, F., Franco, L., Cardona, C., Cordova, L., Toloza, V., 2021. Monitoring and forecasting hazards from a slow growing lava dome using aerial imagery, tri-stereo Pleiades-1A/B imagery and PDC numerical sim-ulation. Earth Planet. Sci. Lett. 564, 116906 ISSN 0012-821X10.1016/j.epsl.2021.11690

Nakada S, Shimizu H, Ohta K (1999) Overview of the 1990–1995 eruption at Unzen Volcano. J Volcanol Geoth Res 89:1–22. https://doi.org/10.1016/s0377-0273(98)00118-8

Newhall CG, Melson WG (1983) Explosive activity associated with the growth of volcanic domes. J Volcanol Geoth Res 17:111–131. https://doi.org/10.1016/0377-0273(83)90064-1

Ogburn, S.E., Loughlin, S.C., Calder, E.S., 2015. The association of lava dome growth with major explosive activity (VEI ≥ 4) Dome-Haz, a global dataset. Bull Volcanol 77. https://doi.org/10.1007/s00445-015-0919-x

Oppenheimer C, Francis PW, Rothery DA, Carlton RWT, Glaze LS (1993) Infrared image analysis of volcanic thermal features: Lás-car Volcano, Chile, 1984–1992. J Geophys Res 98(B3):4269–4286. https://doi.org/10.1029/92JB02134

Ortega, M., et al., 2021. Informe técnico anual Vigilancia del volcán Sabancaya, periodo 2021. Instituto Geológico, Minero y Metalúrgico. Informe técnico N° A7136. 125p.

Pallister, J., McNutt, S.R., 2015. Synthesis of volcano monitoring. The Encyclopedia of Volcanoes 1151–1171. https://doi.org/10.1016/ b978-0-12-385938-9.00066-3

Pararas-Carayannis G (1992) The tsunami generated from the erup-
tion of the volcano of Santorin in the Bronze Age. Nat Hazards
5:115–123. https://doi.org/10.1007/bf00127000

Plank S, Walter TR, Martinis S, Cesca S (2019) Growth and collapse of a Littoral lava Dome during the 2018/19 eruption of Kadovar Volcano, Papua New guinea, analyzed by multi-sensor satellite imagery. J Volcanol Geoth Res 388:106704. https://doi.org/10.1016/j.jvolgeores.2019.106704

Ramos, D., et al. 2019. Informe técnico anual Vigilancia del volcán Sabancaya, periodo 2019. Instituto Geológico, Minero y Metalúrgico. Informe técnico N° A6998. 99p.

Reath K, Pritchard M, Poland M, Delgado F, Carn S, Coppola D, Andrews B, Ebmeier SK, Rumpf E, Henderson S, Baker S, Lundgren P, Wright R, Biggs J, Lopez T, Wauthier C, Moruzzi S, Alcott A, Wessels R, Griswold J, Ogburn S, Loughlin S, Meyer F, Vaughan G, Bagnardi M (2019a) Thermal, deformation, and degassing remote sensing time series (CE 2000–2017) at the 47 most active volcanoes in Latin AMERICA implications for volcanic systems. Journal of Geophysical Research Solid Earth 124:195–218. https://doi.org/10.1029/2018jb016199

Reath K, Pritchard ME, Moruzzi S, Alcott A, Coppola D, Pieri D (2019b) The AVTOD (ASTER VOLCANIC thermal Output database) Latin America archive. J Volcanol Geoth Res 376:62–74. https://doi.org/10.1016/j.jvolgeores.2019.03.019

Reath, K., Pritchard, M., Biggs, J., Andrews, B., Ebmeier, S.K., Bagnardi, M., Girona, T., Lundgren, P., Lopez, T., Poland, M., 2020. Using conceptual models to RELATE multiparameter satellite data to subsurface volcanic processes in Latin America. Geochem Geophys Geosys 21. https://doi.org/10.1029/2019g c008494

Global Volcanism Program, 2017. Report on Sabancaya (Peru) (Craf-
ford, A.E., and Venzke, E., eds.). Bulletin of the Global Volcan-
ism Network, 425. Smithsonian Institution. https://doi.org/10.
5479/si.GVP.BGVN201705-354006

Samaniego P, Rivera M, Mariño J, Guillou H, Liorzou C, Zerathe S, Delgado R, Valderrama P, Scao V (2016) The eruptive chronology of the Ampato-Sabancaya volcanic COMPLEX (Southern Peru). J Volcanol Geoth Res 323:110–128. https://doi.org/10.1016/j.jvolgeores.2016.04.038

Shinohara, H., 2008. Excess degassing from volcanoes and its role on eruptive and intrusive activity. Rev Geophys 46. https://doi.org/10.1029/2007rg000244

Sparks RSJ (1997) Causes and consequences of pressurisation in lava dome eruptions. Earth Planet Sci Lett 150:177–189. https://doi.org/10.1016/s0012-821x(97)00109-x

Stevenson DS (1993) Physical models of fumarolic flows. J Volcanol Geoth Res 57(3-4):139–156. https://doi.org/10.1016/0377-0273(93)90009-G

Stevenson DS, Blake S (1998) Modelling the dynamics and thermody-
namics of volcanic degassing. Bull Volcanol 60:307–317. https://
doi.org/10.1007/s004450050234
Tanaka, H. K. M., T.Uchida, M. Tanaka, M. Takeo, J. Oikawa, T. Ohminato, Y. Aoki, E. Koyama, and H. Tsuji (2009), Detecting a mass change inside avolcano by cosmic-ray muon radiography (muography) first results from measurements at Asama volcano, Japan, Geophys.Res. Lett., 36, L17302, https://doi.org/10.1029/2009GL039448

Valade, S., Ley, A., Massimetti, F., D'Hondt, O., Laiolo, M., Coppola, D., Loibl, D., Hellwich, O., Walter, T.R., 2019. Towards global volcano monitoring using multisensor sentinel missions and artificial intelligence The mounts monitoring system. Remote Sensing 11, 1528. https://doi.org/10.3390/rs11131528

Vallance, J.W., Schneider, D.J., Schilling, S.P., 2008. Growth of the 2004–2006 lava-dome complex at Mount St. Helens, Washington. Professional Paper 169–208. https://doi.org/10.3133/pp17509

Global Volcanism Program, 2013. Volcanoes of the World, v. 4.10.1 (29 Jun 2021). Venzke, E (ed.). Smithsonian Institution. Downloaded 27 Jul 2021. https://doi.org/10.5479/si.GVP.VOTW4-2013.

Wadge G, Cole P, Stinton A, Komorowski J-C, Stewart R, Toombs AC, Legendre Y (2011) Rapid topographic change measured BY high-resolution satellite radar at Soufriere Hills Volcano, Montserrat, 2008–2010. J Volcanol Geoth Res 199:142–152. https://doi.org/10.1016/j.jvolgeores.2010.10.011

Wallace PJ (2001) Volcanic so 2 emissions and the abundance and distribution of exsolved gas in magma bodies. J Volcanol Geoth Res 108:85–106. https://doi.org/10.1016/s0377-0273(00)00279-1

Walter TR, Subandriyo J, Kirbani S, Bathke H, Suryanto W, Aisyah N, Darmawan H, Jousset P, Luehr B-G, Dahm T (2015) Volcanotectonic control of Merapi's lava dome splitting The November 2013 Fracture observed from high Resolution TerraSAR-X data. Tectonophysics 639:23–33. https://doi.org/10.1016/j.tecto.2014.11.007

Walter, T.R., Wang, R., Zimmer, M., Grosser, H., Lühr, B., Ratdomopurbo, A., 2007. Volcanic activity influenced by tectonic earthquakes static and dynamic stress triggering at Mt. Merapi. Geophysical Research Letters 34. https://doi.org/10.1029/2006g01028710

Werner C, Kern C, Coppola D, Lyons JJ, Kelly PJ, Wallace KL, Schneider DJ, Wessels RL (2017) Magmatic degassing, lava dome extrusion, and explosions from Mount Cleveland volcano, Alaska, 2011–2015 insight into the continuous nature of volcanic activity over multi-year timescales. J Volcanol Geoth Res 337:98–110. https://doi.org/10.1016/j.jvolgeores.2017.03.001

Werner C, Rasmussen DJ, Plank T, Kelly PJ, Kern C, Lopez T et al (2020) Linking subsurface to surface using gas emission and melt inclusion data at Mount Cleveland volcano, Alaska. Geochemistry, Geophysics, Geosystems 21(7). https://doi.org/10.1029/2019G0008882

Wooster M (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:83–107. https://doi.org/10.1016/s0034-4257(03)00070-1

Zorn, E.U., Le Corvec, N., Varley, N.R., Salzer, J.T., Walter, T.R., Navarro-Ochoa, C., Vargas-Bracamontes, D.M., Thiele, S.T., Arámbula Mendoza, R., 2019. Load stress controls on directional lava dome growth AT VOLCÁN de colima, mexico. Front Earth Sci 7. https://doi.org/10.3389/feart.2019.00084