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 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 non-linear 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 (SO2), (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 km2 (Samaniego et al., 2016). The Sabancaya activity probably started 6–10 ky ago and built a volcanic edifice of 6–10 km3 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).

Fig. 1
figure 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 Detailed view of the lava dome surface on 29 October 2019. By this time the dome measured approximately 230 × 290 m

Full size image

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 km3 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 km3 per year. Notably, satellite images acquired during the eruption have shown the presence of an open crater in the summit 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 km3 yr−1, 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).

Fig. 2
figure 2

Examples of cloud-free satellite images used in this work to track the evolution of Sabancaya eruption; (a#) CNES/Airbus visible images(C) (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). 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

Full size image
Fig. 3
figure 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” 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 SO2 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

Full size image

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 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 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 (Hplume), 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 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 SO2 flux (here provided in tons day−1). The SO2 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 SO2 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)

$${\mathrm E}_{\mathrm{seismic}}=2\mathrm{\pi \, r}^2\mathrm{c}\frac1{\mathrm A}\int\mathrm S^2\mathrm U\left(\mathrm t\right)^2\mathrm{dt}$$
(1)

where r is the source-station distance, ρ is the density (2600 kg m−3), c is the P wave velocity (3000 m s−1), 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 SO2 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).

Fig. 4
figure 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) SO2 flux is measured by the ground-based DOAS 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 LANDSAT8/SENTINEL2. An increase in SO2 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)

Full size image

After the phreatic explosions in August 2014, the heat flux continued to increase, until the second half of 2015. The SO2 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 SO2 flux raising from ~ 230 tons day−1 in 2014, to ~ 518 tons day−1 in 2015 (Carn et al., 2017). Notably, the increasing degassing was characterized by periodic oscillations in the 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.

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 SO2 measurements have become more continuous and tracked this further acceleration of the unrest phase throughout the rest of the year. Notably, SO2 emissions have tripled in just a few months reaching more than ~ 2000 tons day−1 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 SO2 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. 5 e,f,g).

Fig. 5
figure 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) SO2 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

Full size image

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 SO2 flux showed a decreasing trend albeit remaining at high levels (> 2000 tons day−1), 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) × 106 m3 by December 2017 (Fig. 5c) with an average extrusion rate of ±0.04 ± 0.02 m3 s−1.

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 SO2 flux which stabilized around 2000 tons day−1 (Fig. 5e).

Notably, the overall dome growth observed during phase II was temporally interrupted by three deepening-refilling cycles (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−1 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 SO2 emissions recorded in July 2018, when the flux exceeded 6000 tons day−1, before returning to the typical levels (±1500 tons day−1) 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 m3 s−1. This recovery was characterized by a coeval increase in the SO2 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 SO2 flux (reaching 5000 tons day−1 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 m3 s−1) and was followed by a new pulse of SO2 (> 6000 tons day−1) 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 m3 s−1 and was accompanied by more than 200 explosive events per day. It is interesting to note that during this intense dome growth phase, the SO2 flux remained moderate, with values below 2000 tons day−1. On 26 October 2019, OVI-INGEMMET conducted a drone overflight that allowed to measure a lava dome volume of approximately 4.6 × 106 m3 (red cross in Fig. 5c). By this time, Sentinel1-derived measurements indicated a volume of 4.3. (±1.4) × 106 m3 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 SO2 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−1 (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−1. 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 SO2 flux increased threefold, although after reaching a peak in late August 2020 (”5600 tons day−1) 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 (Vdome), (ii) the volume of emitted ash (Vash), (iii) the volume of degassed magma (Vdegas), (iv) the volume of magma radiating an equivalent amount of thermal energy (Vtherm). 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−3 (Aguilar et al., 2019).

Fig. 6
figure 6

Magma budget of the Sabancaya eruption (all data smoothed over 3 months for clarity). (a) Cumulative volumes of extrusive (Vdome) and explosive activity (Vash). (b) Cumulative volumes of degassed magma (Vdegas) and thermally radiant magma (Vthermal). Note the strong excess of these two volumes to those related to the extrusive (Vdome) and explosive(Vash) 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 (Qdegas, Qthermal) in clear excess to erupted magma fluxes (Qdome, Qash). Note that the curve of Qdome results fragmented during period of dome stability or collapse (Qdome = 0)

Full size image

(i) Volume of extruded magma (Vdome) 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−3 (Tanaka et al., 2009), which applied to the dome of the Sabancaya would produce a DRE volume equal to 0.5 × 106 m3. Unfortunately, the density of the Sabancaya dome is unknown. For this reason, in this work, we assumed a range of densities spanning from 1500 to 2500 kg m−3 to take into account variable emplacement conditions (Fig. 6a).

(ii) Volume of ash (Vash) 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 (Hplume) above the crater is related to the mass of ash (Mash) by

$${H}_{plume}=B{{M}_{ash}}^{0.25}$$
(2)

where B is a model-derived coefficient (in m kg−4) that takes into account the specific thermal conditions of vulcanian plumes ascent (Londono et al., 2018). Here, we inverted Eq. (2) to estimate Mash from the measured Hplume and associated errors (see the “Ash plume height and color (camera network)” section). We used B = 55 m kg−4 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 (Vdegas) this is the magma volume that has degassed the measured mass of SO2 (MSO2) throughout the analyzed period. This is calculated by using the petrological method (Shinohara 2008) according to which

$${V}_{degas}=M{SO}_{2}/\left(2*({\uprho }_{DRE}*\Delta {X}_{S}\right))$$
(3)

where MSO2 is the mass (in kg) of degassed SO2 and ΔXS 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 ΔXS 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 plumbing system of Sabancaya (Gerbe & Thouret, 2004). This model was recently proposed to explain the sustained SO2 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 ΔXS 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 Aleutians 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 (Vthermal) 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}$$
(4)

where crad 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 crad equal to 1.2 (±0.6) × 107 J m−3, which is typical for an andesitic lava dome. It is important to note that the error of crad (±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 Vthermal 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 (Qdome, Qash, Qdegas, and Qthermal).

Excess degassing and excess thermal radiation

Magma budget calculations (Fig. 6) indicate that about 0.002 km3 of magma were extruded to form the “Huk” and “Iskay” lava domes, while 0.04–0.09 km3 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 (Verupt = Vdome + Vash) was 1 or 2 order of magnitude lower than the degassed magma (Vdegas ± 0.25–1.28 km3; Fig. 6b) depending on the assumed degassing scenario. This unbalance (Vdegas >  > Verupt) 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 (Vthermal = 0.04–0.13 km3) 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 (Vthermal > Verupt).

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 reach the base of the extruding lava dome (model 2).

Fig. 7
figure 7

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. 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

Conceptual models of excess degassing and excess thermal radiation of Sabancaya volcano during unrest and eruption (

Full size image

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 chamber (Girona et al., 2014). In this scenario, the rise of the magma column translates into an ever-larger SO2 degassing cell (i.e. the portion of the convective magma column above the SO2 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 (Qdegas) ranges from 2 to 10 m3 s−1 (depending on the degassing scenario considered). This value is compared with an average magma output rate (Qout = Qdome + Qash) lower than 0.08 m3 s−1. Similarly, the volumetric flux needed to explain the thermal emissions (Qthermal) ranges from 0.3 to 1 m3 s−1 (Fig. 6c), which is still 4 to 10 times Qout 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 Qthermal and Qash (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 Qdegas (mafic magma supply scenario) and the maximum limit of Qthermal almost converge for a rate of magma convection equal to 1–2 m3 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 m3 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 Qthermal and Qdegass, 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 Qdegas corresponds to the rate of magma circulation above the exsolution level of SO2 (i.e. within the full convective degassing cell), the thermally-derived magma flux (Qthermal) 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, SO2 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 SO2 degassing (from ~ 200 to ~ 2000 tons day−1) 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 m3 s−1) and collapse, accompanied by variable heat flux (0.5–25 MW), persistent but irregular degassing (700–7000 tons day−1), 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 km3 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 m3) 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 km3 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