Cyclic activity of the Fuego de Colima volcano (Mexico): insights from satellite thermal data and nonlinear models

The Fuego de Colima volcano (Mexico) shows a complex eruptive behavior, with periods of rapid and slow lava dome growth punctuated by explosive activity. We reconstructed the weekly discharge rate average between 1998 and 2018 by means of satellite thermal data integrated with published discharge rate data. By using spectral and wavelet analysis, we found a multiyear long-term, multi-month intermediate-term, and multi-week short-term cyclic behavior during the period of the investigated eruptive activity like that of many other dome-forming volcanoes. We use numerical modeling in order to investigate the nonlinear cyclic eruptive behavior considering a magma feeding system composed of a dual or a single magma chamber connected to the surface through an elastic dyke developing into a cylinder conduit in the shallowest part. We investigated cases in which the periodicity is controlled by (i) the coupled deep–shallow magma reservoirs, (ii) the single shallow chamber, and (iii) the elastic shallow dyke when it is fed by a fixed influx rate or constant pressure. Due to the limitations of the current modeling approach, there is no single configuration that can reproduce all the periodicities on the three different timescales. The model outputs indicate that the observed multiyear periodicity (1.5–2.5 years) can be described by the fluctuations controlled by a shallow magma chamber with a volume of 20–50 km³ coupled with a deep reservoir of ca. 500 km³, connected through a deep elastic dyke. The multi-month periodicity (ca. 5–10 months) appears to be controlled by the shallow magma chamber for the same range of volumes. The short-term multi-week periodicity (ca. 2.5–5 weeks) can be reproduced considering a fixed influx rate or constant pressure at the base of the shallower dyke. This work provides new insights on the nonlinear cyclic behavior of Fuego de Colima and a general framework for comprehension of the eruptive behavior of andesitic volcanoes.

1 Introduction

Lava-dome-forming eruptions are relatively long-lived events, lasting from several months to several decades (e.g., Merapi, Indonesia, Siswowidjoyo et al., 1995; Kelut, Indonesia, De Bèliaz et al., 2012; Fuego de Colima, Mexico, Lamb et al., 2014; Santiaguito, Guatemala, Harris et al., 2003), and are usually punctuated by dome collapses and explosive (Vulcanian) episodes. Discharge rates can change widely over a range of timescales, reflecting the physical mechanisms involved in the transfer of magma to the Earth's surface (Melnik et al., 2009; Odbert and Wadge, 2009).

Volcanoes, including Santiaguito (Guatemala; Harris et al., 2003), Mt. St. Helens (USA; Swanson and Holcomb, 1990), and Soufrière Hills, Montserrat (Voight et al., 1998; Loughlin et al., 2010; Wadge et al., 2010; Nicholson et al., 2011), have shown periodic dome growth. Periodicity in eruption parameters can be complex, showing systematic or nonsystematic temporal changes as the eruption progresses (Denlinger and Hoblitt, 1999; Costa et al., 2007a; Melnik et al., 2009; Bernstein et al., 2013; Wolpert et al., 2016), and can be characterized by short-, intermediate-, and long-term wavelengths (Costa et al., 2007a, 2012, 2013; Melnik et al., 2009; Melnik and Costa, 2014). Short- and intermediate-term periodicities (hours or weeks) are generally explained by the upper conduit pressurization related to the nonlinear ascent of magma flow (Denlinger and Hoblitt, 1999; Melnik and Sparks, 1999; Voight et al., 1999; Wylie et al., 1999; Ozerov et al., 2003; Lensky et al., 2004; Costa et al., 2007a, b, 2012; Kozono and Koyaguchi, 2009, 2012). This is because the lower part of the dyke conduit can act as a capacitor that allows magma to be stored temporarily and released during the more intense phase of discharge (Costa et al., 2007a, b, 2012, 2013; Melnik et al., 2009). The long-term periodicity, with timescales from several months to decades (Voight et al., 2000; Belousov et al., 2002; Sparks and Young, 2002; Wadge et al., 2006), is usually controlled by pressure variations in magma reservoirs (Barmin et al., 2002; Costa et al., 2007b; Melnik et al., 2009; Melnik and Costa, 2014).

Since historical times, the Fuego de Colima volcano (Mexico; Fig. 1a) has been characterized by decade-lasting cycles of dome growth alternating with Vulcanian explosions and ending with sub-Plinian eruptions. The last two occurred in 1818 and 1913 (Luhr, 2002; Saucedo et al., 2005; Norini et al., 2010; Heap et al., 2014; Massaro et al., 2018b). The most recent cycle started after the 1913 eruption, and it is characterized by lava domes extruded with minor seismicity at high magma temperatures (960–1020 ^∘C) (Savov et al., 2008). As for other dome eruptions (Sparks, 1997), dome growth at Fuego de Colima can be explained by complex nonlinear pressure variations during magma ascent from magma reservoirs (e.g., Melnik and Costa, 2014), cooling, crystallization, degassing (e.g., Melnik and Sparks, 1999; Lensky et al., 2004; Nakanishi and Koyaguchi, 2008; Kozono and Koyaguchi, 2012), and upper conduit geometric configurations characterized by multiple pathways (e.g., Lavallee et al., 2012; Reubi et al., 2015).

Two magma chambers located at different depths characterize the feeding system of the Fuego de Colima volcano (Fig. 1b), with roofs located at ca. 6 (shallow magma chamber) and ca. 15 km (deep magma chamber) of depth, as indicated by petrographic studies (Macìas et al., 2017) and geophysical data (Spica et al., 2017). The purpose of this study is to investigate the existence of patterns of fluctuations in discharge rates during the 1998–2018 eruption activity at the Fuego de Colima volcano. The available geological, geophysical, and petrological data for this recent activity provide a remarkable opportunity to improve the characterization and our understanding about the physical processes underlying the cyclic extrusion of lava domes. In particular, we used thermal remote sensing data along with published effusion rates for reconstructing the oscillatory magma discharge rate behavior of effusive activity at Colima.

The availability of satellite thermal images in the last decade has strengthened the use of thermal data for observing volcanic activity (Ramsey and Harris, 2013), especially in studying the relationships with lava discharge rates (Coppola et al., 2009; Harris et al., 2010; Garel et al., 2012). Coppola et al. (2013) propose that the radiant density of effusive–extrusive activity can be used to estimate lava discharge rates and erupted volumes by means of empirical relationships based on the SiO₂ content of the erupted lava. Although still under debate, the so-called “thermal approach” (Dragoni and Tallarico, 2009) offers a good way to monitor volcanic activity, especially when direct observations are limited or absent. Here we focus our attention on the dynamics of fluctuations in the magma discharge rate at different timescales at the Fuego de Colima volcano during 1998–2018. By using time series analytical techniques (i.e., Fourier and wavelet analysis) we have identified three fundamental periodicities in subsets of the time series: (i) long term (ca. 1.5–2.5 years), (ii) intermediate term (ca. 5–10 months), and (iii) short term (ca. 2.5–5 weeks), similar to those observed at many lava dome eruptions (Costa et al., 2012; Melnik and Costa, 2014; Christopher et al., 2015). These periodicities were compared with numerical simulations provided by the model of Melnik and Sparks (2005) as generalized by Costa et al. (2007a) to account for the presence of a shallow dyke, as well as with Melnik and Costa (2014) for describing the control of a coupled dual chamber system.

Numerical modeling of the different parts of the plumbing system can successfully reproduce the first-order cyclic behavior of Fuego de Colima during the 1998–2018 eruption activity. Our results highlight the fact that the dual magma chamber dynamics control the long-term periodicity evident during 2002–2006 and 2013–2016, while the single magma chamber dynamics are more effective to explain the intermediate-term periodicity in the same periods. Finally, the shallow dyke dynamics regulate the multi-week cycles observed during 2002–2006 and 2011–2016.

The present work is divided in five main sections. The first describes the historical activity of the Fuego de Colima, with particular attention to the recent period, from 1998 to 2018. The second section describes the methods applied to the dataset composed of satellite thermal data integrated with published data. The third section is dedicated to the input and target data used for numerical simulations. The fourth section presents the results obtained by the spectral and wavelet analyses. This latter allows us to establish significance levels for the wavelet power spectrum. The periodicities observed in this spectrum were compared to the results obtained by numerical simulations. The last section contains a discussion on the eruptive behavior that occurred at Fuego de Colima during 1998–2018, providing new insights from observational data and nonlinear models.

Figure 1 (a) Digital elevation model of the Colima Volcanic Complex (NC: Nevado de Colima volcano; FC: Fuego de Colima volcano) and Colima Rift with the main tectonic and volcano–tectonic structures (modified from Norini et al., 2010). In the inset, the location of the Colima Volcanic Complex (CVC) within the Trans-Mexican Volcanic Belt (TMVB) is shown in the frame of the subduction-type geodynamic setting of Central America. (b) Schematic view of the conduit feeding system framework used for numerical simulations (modified after Melnik and Costa, 2014).

2 The historical activity of the Fuego de Colima volcano

Since historical times Fuego de Colima has represented the most active volcano in Mexico, posing a serious threat to all surrounding populations (Cortés et al., 2005; Gavilanes-Ruiz et al., 2009; Bonasia et al., 2011; Roverato et al., 2011). The earliest accounts of the volcano activity can be found in Historia Antigua de Mexico (Clavijero, 1780), wherein the destructive effects of its explosive activity are carefully described (Brèton-Gonzalez et al., 2002).

The historical activity of Fuego de Colima was described and interpreted by several authors (Luhr and Carmichael, 1980; Medina-Martinez, 1983; De la Cruz-Reyna, 1993; Brèton-Gonzalez et al., 2002; Luhr, 2002). The Fuego de Colima has shown a transitional eruptive behavior spanning from effusive to explosive activity, dominated by dome growth and Vulcanian eruptions. Occasionally sub-Plinian events occurred (1576, 1606, 1690, 1818, and 1913), indicating a recurrence time of approximately 100 years (De la Cruz-Reyna, 1993; Luhr, 2002; Saucedo et al., 2005; Gavilanes-Ruiz et al., 2009; Massaro et al., 2018b). The sub-Plinian event that occurred in 1913 (Saucedo et al., 2010) is the largest historical eruption and it has been used as a benchmark for volcanic hazard studies (Martin del Pozzo et al., 1995; Saucedo et al., 2005; Bonasia et al., 2011).

The 1998–2018 eruptive activity

The 1998–2018 period is the only period of post-1913 activity for which there is sufficiently available information to explore the cyclic activity of Fuego de Colima.

Different periods of effusion (domes and lava flows) punctuated by Vulcanian eruptions and dome collapses characterized the volcano activity between 1998 and 2018 (Savov et al., 2008; Varley et al., 2010a; Hutchison et al., 2013; Mueller et al., 2013; Zobin et al., 2015; GVP, 2017). The duration of extrusive activity and the magma discharge rate varied through time, generally divided into five eruptive phases up to 2015: (I) 1998–1999; (II) 2001–2003; (III) 2004–2005; (IV) 2007–2011; (V) 2013–2015 (Zobin et al., 2015; Aràmbula-Mendoza et al., 2018).

The first dome extrusion of the 1998–1999 phase started in November 1998 and quickly filled the 1994 explosion crater, forming lava flows that descended the southern flanks of the Fuego de Colima cone during most of 1999 (<5 m³ s⁻¹ on average for Mueller et al., 2013; 4.11 m³ s⁻¹ on average for Reubi et al., 2013). At the beginning, this dome grew rapidly (ca. 4.4 m³ s⁻¹), reaching a volume of ca. 3.8×10⁵ m³ in 24 h. During this period the effusion rate reached a peak value around 30 m³ s⁻¹ (Navarro-Ochoa et al., 2002; Zobin et al., 2005; Reubi et al., 2015) and showed a cyclic damped behavior soon after.

During 1999–2001 a series of explosions destroyed the dome and excavated a large apical crater (Brèton-Gonzalez et al., 2002). A slow outpouring of lava (<1 m³ s⁻¹ for Mueller et al., 2013, 0.17 m³ s⁻¹ for Reubi et al., 2013, 2015) resumed in May 2001 and continued for 22 months. In February 2002, the lava dome overflowed the crater rims, producing lava flows. During this eruptive phase, the magma extruded from three separate vents with only minor explosive activity at a rate of ca. 0.9 m³ s⁻¹ (GVP, 2002). Vulcanian explosions dismantled the dome during July and August 2003 (GVP, 2003).

In September 2004, low-frequency seismic swarms heralded the onset of the new effusive phase (Varley et al., 2010a; Aràmbula-Mendoza et al., 2011; Lavallee et al., 2012) with a small increase in the average discharge rate of 0.6 m³ s⁻¹ (Reubi et al., 2013, 2015). The lava dome building occurred from the end of September until the beginning of November, with a magma effusion rate up to 7.5 m³ s⁻¹ in October (Zobin et al., 2008, 2015). The effusive activity was accompanied and followed by intermittent Vulcanian explosions. The explosive activity diminished in intensity during December 2004–January 2005.

From February to September 2005, effusion and large explosions occurred. In the following months, small, short-lived domes were observed, with an estimated effusion rate between 1.2 and 4.6 m³ s⁻¹ (Varley et al., 2010b; Reubi et al., 2015). In May and June, the explosive activity produced pyroclastic density currents reaching distances up to 5.4 km from the volcano summit (Varley et al., 2010a).

In February 2007, a new lava dome began to grow and explosions were reported in the period between January 2009 and March 2011. The 2007–2011 period of dome extrusion represents the slowest growth rate in the recent history of Fuego de Colima. Hutchison et al. (2013) calculated a mean effusion rate of ca. 0.02 m³ s⁻¹ from 2007 to 2010 using digital photographic data, in good accordance with Zobin et al. (2015), who reported extrusion rates of 0.03 m³ s⁻¹ during 2007. Mueller et al. (2013) estimated the magma extrusion rate between 0.008±0.003 and 0.02±0.007 m³ s⁻¹ during 2010, which dropped down to 0.008±0.003 m³ s⁻¹ again in March 2011. On 21 June 2011 an explosion heralded the cessation of dome growth and marked the end of the effusive period. After 1.5 years of rest, in January 2013 a sequence of explosions cored out the 2011 dome and generated pyroclastic density currents that reached distances of up to 2.8 km from the summit (GVP, 2013). From March to October, the calculated discharge rate was in the range of 0.1–0.2 m³ s⁻¹ (Reyes-Dávila et al., 2016).

Successively, the middle to low explosive activity took place up to February–March 2014 until a new pulse of magma was observed in July, with an approximate rate of 1–2 m³ s⁻¹ (Aràmbula-Mendoza et al., 2018). On 11 January 2015, a new lava dome was observed inside the crater (Thiele et al., 2017) and its growth continued until July, with an effusion rate of ca. 0.27 m³ s⁻¹ (Zobin et al., 2015). Between 10 and 11 July 2015 the recent dome was destroyed by the most intense activity since the 1913 eruption (Capra et al., 2016; Reyes-Dávila et al., 2016).

In the 2013–2015 period, the average extrusion rate was ca. 0.2 m³ s⁻¹ (Thiele et al., 2017), with peak values >10 m³ s⁻¹ (Varley, 2015). After that, the eruptive activity ceased until January 2016 when daily ash plumes started to occur along with active lava flows and explosions. In early July a new dome began to grow, overtopping the crater rim. A large explosion was recorded on 10 July 2016, followed by daily and multiple daily ash plumes up to the end of year. Multiple flows descended from the lava dome during September–December. In 2017 frequent strong explosions and ash emissions were recorded until March. Through June decreasing seismicity and minor landslides were reported with no evidence of effusive activity or new dome growth (GVP, 2017).

Here we provide a more systematic overview of the 1998–2018 eruption activity obtained by satellite thermal data along with some published data, explained in the following section.

3 Methods

We analyzed the thermal energy spectrum of the Fuego de Colima volcano available from March 2000 to October 2018 from the Middle Infrared Observation of Volcanic Activity (MIROVA) hot-spot detection system (Coppola et al., 2016). The period 1998–1999 was integrated using published discharge rates (Navarro-Ochoa et al., 2002; Zobin et al., 2005). The MIROVA system is based on the near-real-time (NRT) analysis of the MODerate resolution Imaging Spectroradiometer (MODIS) data, distributed by the LANCE-MODIS data system (http://modis.gsfc.nasa.gov/, last access: 22 August 2019).

The thermal emission from an object is attenuated by the atmosphere resulting from absorption by gases and scattering by particles. The MIROVA system focuses on the middle infrared region (MIR), which shows the lowest attenuation levels, to better detect and analyze thermal radiation emitted from volcanic sources. While the standard MODIS forward processing delivers Aqua and Terra images within 7–8 h of real time, LANCE-MODIS allows for the creation of MIROVA radiant flux time series within 1–4 h from the satellite overpass (http://www.mirovaweb.it, last access: 22 August 2019). This thermal data collection was converted into lava discharge rate estimates and integrated with some published data in order to reconstruct the weekly mean discharge rate spectrum from 1998 to 2018 (Fig. 2a).

Figure 2 Dataset of the averaged discharge rates of Fuego de Colima during 1998–2018, derived from the MIROVA thermal data (black points) and published data (blue crosses) (Navarro-Ochoa et al., 2002; Zobin et al., 2005; Reubi et al., 2013; Mueller et al., 2013; Varley, 2015; Reyes-Dávila et al., 2016; Thiele et al., 2017; GVP, 2002, 2017). Values >0.1 (m³ s⁻¹) are considered to be “high” (dark blue area) and values <0.1 (m³ s⁻¹) are “low” discharge rates (light blue area). The 0.01 (m³ s⁻¹) value is the threshold under which the MIROVA system does not provide reliable data (blue line); (a) weekly average discharge rates. The boxes contain symbols of volcanological observations reported in the literature; (b) monthly average discharge rates; (c) yearly average discharge rates.

In this work, we refer to Coppola et al. (2013), who describe the relationship between the heat lost by lava thermal radiance variations and discharge rates by means of a unique empirical parameter. They compared the energy radiated during several distinct eruptions to the erupted lava volumes (m³). The relationship between the volcanic radiated energy (VRE) and the erupted volume was defined by introducing the concept of radiant density (c_rad; J m⁻³).

This parameter is analyzed as a function of the SiO₂ content and the bulk rheological properties of the related lava bodies. It is strongly controlled by the characteristic thickness of the active lavas at the time of a satellite overpass, whereas the effects of variable degrees of insulation, morphology, and topographic conditions produce only a limited range of variability (±50 %) (Coppola et al., 2013). For the Fuego de Colima we used a value of $c_{\mathrm{rad}}=\mathrm{3.90}\times {\mathrm{10}}^{\mathrm{7}}$ (J m⁻³) for a SiO₂ content of 59.6 % (Savov et al., 2008; Coppola et al., 2013). We obtained the cumulative volumes of effusion per year (from 2000 to 2018) considering the ratio between the average VRE estimations and c_rad.

It is important to stress that the instrumental limit of the MIROVA system is not able to detect thermal anomalies below 0.5–1 MW. Since we used a radiant density c_rad of 3.90×10⁷ J m⁻³, the minimum reliable value of the discharge rate is 0.01 m³ s⁻¹ (Coppola et al., 2013). As reported by Coppola et al. (2016), the thermal data obtained from MIROVA are not corrected due to the presence and/or attenuation of clouds. For this reason, the estimates of effusion rates and volumes are to be considered minimum estimates. Because the 2002–2006 and 2013–2016 intervals are the most active in the analyzed period, we firstly applied the Fourier analysis to the monthly average discharge rates (Fig. 2b) of these time intervals in order to explore the modal spectrum of the signal.

Although Fourier analysis is well suited to the quantification of constant periodic components in a time series, it cannot recognize signals with time-variant frequency content. Whereas a Fourier transform analysis may determine all the spectral components embedded in a signal, it does not provide any information about the timing of occurrence. To overcome this problem, several solutions have been developed in the past decades that are able to represent a signal in the time and frequency domain at the same time.

The aim of these approaches is to expand a signal into different waveforms with local time–frequency properties well adapted to the signal structure (Cazelles et al., 2008). In order to get information on the amplitude of the periodic signals within the Fuego de Colima (MIROVA) time series, we performed a wavelet analysis by decomposing the weekly time series (Fig. 2a) into time–frequency space (Fig. 3).

Wavelet analysis is a powerful tool used in many scientific fields (i.e., ecology, biology, climatology, geophysics) and engineering. It is especially relevant to the analysis of nonstationary systems (i.e., systems with short-lived transient components; Cazelles et al., 2008). In particular, the wavelet analysis is well suited for investigations of the temporal evolution of aperiodic and transient signals (Lau and Weng, 1995; Mallat, 1998). For this study, practical the details of applying wavelet analysis were taken from Torrence and Compo (1998) and Odbert and Wadge (2009). It is worth noting that wavelet analysis considers a wave that decays over a finite time and whose integral over infinite time is zero. Many forms of wavelets (called wavelet “functions”, ψ(η), or “mother functions” that depend on a nondimensional time parameter η) have been designed for analytical use (Farge, 1992; Weng and Lau, 1994; Daubechies, 1992), each with its own characteristics that make it suitable for certain applications. The choice of the wavelet can influence the time and scale resolution of the signal decomposition. Wavelet analysis is popular in geosciences (Trauth, 2006) as it does not require any a priori understanding of the system generating the time series. Our time series (weekly average discharge rates acquired by the MIROVA system), called (x_n), has equal time spacing (δt=7 d) and an equal number of points $n=\mathrm{0}\mathrm{\dots }N-\mathrm{1}$. Using the approximately orthogonal Morlet function wavelet ψ(η) (it must have a zero mean and be localized in both time and frequency space; Farge, 1992), we define the wavelet transform W_n(s) as the convolution of x_n with a scale (s) and translated version of ψ₀(η) (mother function). The formula is as follows:

$$\begin{array}{}\text{(1)}& W_n\left(s\right)=\sum _{n^{\prime }=\mathrm{0}}^{N-\mathrm{1}}{\mathit{\chi }}_n^{\prime }{\mathrm{\Psi }}^{*}\left[{\displaystyle \frac{(n^{\prime }-n)\mathit{\delta }t}{s}}\right],\end{array}$$

where (*) indicates the complex conjugate. The scale s should be equal to approximately 2δt, according to the Nyquist theorem. Therefore, the smallest wavelet we could possibly resolve is 2δt, and thus we choose s₀=14 d.

Generally, ψ(η) is a complex function, and therefore the wavelet transform is also complex. It is possible to reconstruct the “local” wavelet power spectrum as the absolute value squared of the wavelet coefficients, $\left|\right(W_n\left(s\right)){|}^{\mathrm{2}}$. The way to compute the wavelet transform for a time series is to find the Fourier transform of both the wavelet function (Morlet in our case) and the time series.

Following Torrence and Compo (1998), we made the normalization by dividing by the square root of the total wavelet variance (σ²).

Usually, a periodic component in a time series may be identified in a power spectrum if it has distinctly greater power than a mean background level (that would correspond to Gaussian background noise) (Odbert and Wadge, 2009). However, the spectra generated from many geophysical systems indicate that the noise in time series data tends not to have a Gaussian distribution (Vila et al., 2006), but it can be better described by colored noise, specifically red noise (Fougere, 1985).

For this reason, we use a simple model for red noise given by the univariate lag-1 autoregressive or Markov process (Torrence and Compo, 1998) in order to determine the significance levels for our wavelet spectrum. These background spectra are used to establish a null hypothesis for the significance of a peak in the wavelet power spectrum. The null hypothesis is defined for the wavelet power spectrum considering that the time series has a mean power spectrum: if a peak in the wavelet power spectrum is significantly above this background spectrum, then it can be assumed to be a true feature with a certain percentage of confidence.

For definitions, significant at the 5 % level is equivalent to the 95 % confidence level (Torrence and Compo, 1998). The confidence interval is defined as the probability that the true wavelet power at a certain time and scale lies within a certain interval about the estimated wavelet power (Torrence and Compo, 1998). Because we deal with finite-length time series, errors occur at the beginning and end of the wavelet power spectrum. A solution is to pad the end of the time series with zeroes to bring the total length N up to the next-higher power of two, thus limiting the edge effects. However, padding with zeroes introduces discontinuities at the end points and, especially towards larger scales, decreases the amplitude near the edges as more zeroes enter the analysis (Torrence and Compo, 1998).

The cone of influence (COI) is the region of the wavelet spectrum beyond which edge effects become important. The criterion for applying wavelet analysis is very similar to those employed with classic spectral methods. In other words, the wavelet transform can be regarded as a generalization of the Fourier transform, and by analogy with spectral approaches, we compute the local wavelet power spectrum as described above. Successively, this can be compared with the “global” wavelet power spectrum, which is defined as the averaged variance contained in all wavelet coefficients of the same frequency (Torrence and Compo, 1998; Cazelles et al., 2008).

Numerical simulations have been carried out using the magma flow model of Melnik and Costa (2014), who generalized the model proposed by Melnik and Sparks (2005) and Costa et al. (2007a) for a magma chamber connected to a dyke that develops into a cylindrical conduit near the surface.

In particular, the model of Melnik and Costa (2014) accounts for the possibility of a dual magma chamber system. The model accounts for rheological changes due to volatile loss and temperature-driven crystallization. These processes are both effective during dome extrusion eruptions because of the typically low magma ascent velocities (from millimeters to a few centimeters per second), which can result in magma transit times from days to weeks. These ascent times are often comparable with those of crystal nucleation and growth (Melnik and Sparks, 1999, 2005; Costa et al., 2007a).

Figure 3 (a) Local wavelet spectrum normalized by 1∕(σ)² ((σ)²; (m³ s⁻¹)²). The left axis is the period (in years). The bottom axis is time (in years). The shaded contours are at normalized variances of 0.5, 1, 2, and 4 (m³ s⁻¹)². The thick black contour encloses regions of greater than 95 % confidence for a red-noise process with a lag-1 coefficient of 0.72. It shows three orders of periodicities: long-term (ca. 1.5–2.5 years), intermediate-term (ca. 5–10 months) during 2002–2006 and 2013–2016, and short-term (ca. 2.5–5 weeks) during 2001–2006 and 2011–2016. The blue line indicates the “cone of influence”, outside of which edge effects become important; (b) global wavelet power spectrum. The green dotted line represents the best-fitting red-noise spectrum at the 95 % confidence level.

Table 1 Input parameters used in numerical simulations.

4 Input and target data for numerical simulations

4.1 Geometrical configurations of the magma plumbing system

Within the physical framework used in Melnik and Costa (2014), the model (Fig. 1b) consists of two elastic magma chambers located at different depths, with chamber pressures P_chs and P_chd able to drive the magma ascent in an elliptical cross-sectional volcanic conduit (approximating a dyke). Near the surface the conduit develops into a cylinder at depth L_T (called the “transition level”).

Numerical simulations were carried out considering the shallower magma chamber (single magma chamber configuration) or the double magma chamber (see Tables 2–3). The single magma chamber model considers a conduit feeding system composed of a shallow dyke (d_s) that connects the magma chamber to a shallower cylinder, in agreement with geological and geophysical evidence from different volcanoes (Melnik and Sparks, 2005; Costa et al., 2007a; Melnik et al., 2009; Melnik and Costa, 2014). The double magma chamber model includes the addition of a deep reservoir connected to the shallow chamber through a deep elastic dyke (d_d) (Fig. 1b).

Table 2 Conduit geometry parameters using a single magma chamber model.

Table 3 Parameters used for simulations carried out with a dual magma chamber model.

In order to reproduce the observed fluctuations in discharge rates recorded in some periods of the 1998–2018 eruption activity, we considered a discharge rate regime in which the period of pulsations is controlled by the elasticity of the shallow dyke and a discharge rate regime in which the periodicity is controlled by the volume of the single or dual magma chamber(s) (Barmin et al., 2002; Melnik and Sparks, 2005; Costa et al., 2007a; Melnik and Costa, 2014). In Appendix Figs. A1 and A2 we report some test simulations in order to show the control of the most sensitive parameters (i.e., water content in magma, dyke dimensions, volume of magma chamber, magma influx rate into the magma chamber) affecting the model outputs in the case of the single magma chamber model. The volumes of the magma chamber (V_ch) range from 20 to 50 km³ and the width of the feeder dyke 2a varies from 200 to 400 m (Massaro et al., 2018b).

In Appendix Fig. A3 we show the sensitivity test aimed to explore a broad range of chamber volumes and aspect ratios in the case of the double magma chamber configuration. The deep chamber has its top at 15 km of depth; it is pressurized and fed from below by a constant influx Q_in,d.

The volumes of the shallow magma chamber (V_chs) range from 30 to 50 km³, and the volumes of the deep magma chamber (V_chd) from 550 to 750 km³, according to geophysical data (Spica et al., 2017). The aspect ratios for the shallow and deep magma chambers (AR_s−AR_d) varied from 1 to 2. For each run included in sections 1–3 of Fig. A4, we used a fixed influx $Q_{\mathrm{in},\mathrm{d}}=\mathrm{2.3}$ m³ s⁻¹ and variable widths of the deeper dyke (2a_0d) from 200 to 3000 m (representative of weak to strong coupling of the magma chambers; Melnik and Costa, 2014). The lower dyke thickness 2b_0d is not used as input data to the model as it changes as a function of local pressure conditions, and therefore it does not appear in the diagrams. In section 4 of Fig. A3 we show two sets of runs having Q_in,d equal to 1 and 3 m³ s⁻¹, respectively, and the following fixed parameters: AR_s and AR_d=1, V_chd=650 km³, V_chs=40 km³.

4.2 Petrological data

Erupted products at Fuego de Colima are chemically intermediate and primarily andesitic lavas with ca. 61 wt. % SiO₂ (Lavallee et al., 2012). The observed dome growth phases are usually fed by prolonged magma ascent times, which allow for efficient degassing and crystallization. This is in agreement with the low mean porosity (14,%–16 %, e.g., Lavallee et al., 2012; Farquharson et al., 2015; Mora et al., 2002; 0.1 wt. %–2.5 wt. % for 1998–2005 products, Reubi and J., 2008).

Dome lava currently erupted exhibits a range of crystallinities (phenocrysts, 20 % vol.–30 % vol.; microlites, 25 % vol.–50 % vol.), and the groundmass constitutes as much as 68 % vol. (Luhr, 2002). The andesites show a porphyritic texture with plagioclase (13 % vol.–25 % vol.), orthopyroxene (2 % vol.–4 % vol.), clinopyroxene (3 % vol.–4 % vol.), and minor hornblende (<0.5 %) and Fe–Ti oxides (ca. 2 % vol.). Olivine occurs rarely as xenocrysts (Lavallee et al., 2012). As reported in Melnik and Costa (2014), the magma viscosity μ is calculated according to Costa et al. (2007a) considering the melt viscosity, μ_m, times a correction for the effects of crystallinity, θ, and for the bubbles, η. The formula is

$$\begin{array}{}\text{(2)}& \mathit{\mu }={\mathit{\mu }}_{\mathrm{m}}(c,T)\mathit{\theta }\left(\mathit{\beta }\right)\left(\mathit{\eta }\right)(\mathit{\alpha },\mathrm{Ca}),\end{array}$$

which depends on the melt viscosity μ_m (that is a function of the water content c and temperature T), on the crystal content β, on bubble fraction α, and on bubble capillarity number (Ca). The rheological model is described in detail in Costa et al. (2007a). Table 1 summarizes the value ranges used for the input parameters of the model.

5 Results

In Fig. 2 we show the average discharge rates at the Fuego de Colima volcano from November 1999 to October 2018. Here we define “high” discharge rate values as >0.1 m³ s⁻¹ (highlighted as dark blue areas). All values below 0.1 m³ s⁻¹ are considered “low” discharge rates (light blue areas).

Volcanological observations are reported at the top and the bottom of the diagram. It is worth noting that the high and low explosive activity correspond to the high and low discharge rate, respectively. In addition, we distinguish between lava flows and lava domes according to the dominant emplacement style typical of each eruption, as well as between low-magnitude (i.e., ash plumes, gas emissions) and high-magnitude (i.e., strong explosions, Vulcanian eruptions) explosive activity.

The weekly average discharge rates represent the complete dataset used in this study and are reported in Fig. 2a. These data have been calculated by using the MIROVA data (black dots) for the 2000–2018 period and complemented with published data (blue crosses) for the 1998–1999 period (Navarro-Ochoa et al., 2002; Zobin et al., 2005). Even if the data detection of satellite thermal energy represents a continuous spectrum of information, it is worth noting that it suffers from some limitations connected to cloud covering, magma composition, rheology, and the emplacement of the investigated lava body due to topographic conditions (Harris and Rowland, 2009; Harris et al., 2010; Coppola et al., 2009).

Figure 2b shows the monthly discharge rate spectrum from 1998 to 2018 using the MIROVA dataset (black dots), integrated with available published data (blue crosses) (Navarro-Ochoa et al., 2002; Zobin et al., 2005; Capra et al., 2010; Varley et al., 2010a; Sulpizio et al., 2010; James and Varley, 2012; Hutchison et al., 2013; Reubi et al., 2013; Varley, 2015; Reyes-Dávila et al., 2016; Thiele et al., 2017; GVP, 2000, 2017). Figure 2c summarizes the yearly average discharge rates from the MIROVA dataset, highlighting the good agreement with the available average estimation of yearly discharge rates from the literature (Mueller et al., 2013; Reyes-Dávila et al., 2016; Aràmbula-Mendoza et al., 2018; GVP, 1998, 2017).

5.1 Fourier analysis

We applied Fourier analysis to the 1998–2018 dataset (Fig. 2a). In particular, we chose two time windows: (i) the 2002–2006 period, which showed two periodic components, T₀=24.70 and T₁=6.17, corresponding to ca. 24 and ca. 6 months, respectively (Appendix Fig. A4a), and (ii) the 2013–2016 period that provided similar results of T₀=24.94 and T₁=6.23, corresponding to ca. 25 and ca. 6 months, respectively (Appendix Fig. A4b).

5.2 Morlet wavelet analysis

The whole analyzed dataset is composed of 825 data points, representing the time evolution of the oscillating components of the 1998–2018 eruptive activity (Fig. 2a). Figure 3a shows the normalized local wavelet power spectrum of the signal. The color scales for power values vary from light orange (low values) to dark red (high values). The thick black contours represent the 95 % confidence level. The blue line indicates the cone of influence (COI) that delimits the region not influenced by edge effects.

From this analysis, it is easy to observe three main periodicities during the 2002–2006 and 2013–2016 periods: (i) long-term periodicity of ca. 1.5–2.5 years; (ii) intermediate-term periodicity of ca. 5–10 months; and (iii) short-term periodicity of ca. 2.5–5 weeks.

The volcanological observations (for high and low discharge rates) are also reported in order to provide a closer link between the observational datasets and the identification of a frequency change in the extrusion rate time series. The short-term periodicity is also present in 2011 (Fig. 3a). Figure 3b shows the global wavelet spectrum corresponding to the local wavelet power spectrum plotted in Fig. 3a. The green dashed line shows the position of the best-fitting red-noise model at the 95% confidence level.

5.3 Numerical simulations

Appendix Figs. A1–A3 provide some sensitivity tests in order to explore the effects of different parameters on discharge rate fluctuations for the single (Figs. A1–A2) and dual magma chamber models (Fig. A3). In particular, in Appendix Fig. A1 the general steady-state solution of the numerical model is reported, with both stable and unstable branches (e.g., Melnik et al., 2009; Nakanishi and Koyaguchi, 2008; the cyclic behavior can occur only between 2 and 4 m³ s⁻¹ for the fixed input data – panel a). Varying the width of the shallow dyke 2a (from 200 to 400 m) and water content in the melt phase, we observed how the unstable branch changes its shape. This implies different periods of possible oscillations in the discharge rate (panels b–c).

Appendix Fig. A2 provides a set of simulations carried out varying the width of the shallow dyke 2a (panel a). The resulting periodicities vary from ca. 1000 d (2a=200 m) ca. 500 d (2a=300 m) to ca. 250 d (2a=400 m). These results highlight the negative correlation between dyke widths and periods of oscillation (Costa et al., 2007a). In this case, the variable widths influence the intensity and periodicity of discharge rates: the wider the dyke, the lower the intensity and periodicity of discharge rates.

Differences in the amplitude of oscillations are observed in panel (b), highlighting a positive correlation between the volume of the magma chamber V_ch and periodicities. Periodicities of ca. 500 d correspond to 20–30 km³, while larger values of ca. 970 and ca. 1176 d are provided for 40 and 50 km³, respectively. In panel (c), we also report a set of simulations considering the modeled discharge rate controlled by the elasticity of the shallower dyke with fixed influx rates Q_in (in the range of 0.01–0.1 m³ s⁻¹).

Appendix Fig. A3 contains four sections dedicated to the sensitivity tests for the dual magma chamber model (Table 3). As reported in Melnik and Costa (2014), the dual chamber model shows cyclic behavior with a period that depends on the intensity of the influx rate and the chamber connectivity (described as the horizontal extent of the dyke connecting the two chambers).

For a weak connectivity, the overpressure in the deeper chamber remains nearly constant during the cycle and the influx of fresh magma into the shallow chamber is also nearly constant. For a strong connectivity between the two chambers, their overpressures increase or decrease during the cycle in a synchronous way. Influx into the shallow chamber stays close to the extrusion rate at the surface (Melnik and Costa, 2014). We explored different cases considering various fixed parameters as follows: (i) volumes of the shallow and deep magma chambers (V_chs=40 km³, V_chd=650 km³); (ii) aspect ratios (AR_s=1, AR_d=1) and the deep magma chamber volume (V_chd=650 km³); and (iii) aspect ratios (AR_s=1, AR_d=1) and the shallow magma chamber volume (V_chs=40 km³). For cases (i), (ii), and (iii), the deep influx rate Q_in,d has fixed values from 3 to 1 m³ s⁻¹.

In conclusion, these sensitivity tests showed the passage from weakly connected magma chambers (lack of simultaneous oscillation of Q_in,s and Q_out) when 2a_0d=200 m to strongly connected magma chambers (synchronous oscillations of Q_in,s and Q_out) when 2a_0d=3000 m.

Figure 4 reports the results of numerical simulations aimed to reproduce the Fuego de Colima fluctuations during 1998–2018. Figure 4a shows a representative example of a time-dependent solution for a discharge rate controlled by the elasticity of the shallower dyke.

Simulations were carried out using fixed values of pressure (blue line) and influx rate (green line) at the source region of the shallower dyke, which is ca. 6000 m long. The dyke has width 2a=400 m, thickness 2b=2 m, and a dyke–cylinder transition T_L at 1300 m of depth. The magma chamber volume is fixed to 30 km³. Solutions present periodicities from 16 to 40 d, in agreement with the weekly periodicities of ca. 38–18 d (ca. 2.5–5 weeks) derived from the wavelet analysis (Fig. 3a).

Figure 4b describes a representative example of the single magma chamber model simulations. We set up the magma feeding system to be composed of a 6500 m long dyke with a width 2a=600 m, thickness 2b=4 m, and a dyke–cylinder transition T_L fixed at 1000 m of depth. The chamber has a volume fixed to 30 km³ and receives a constant $Q_{\mathrm{in},\mathrm{s}}=\mathrm{2.3}$ m³ s⁻¹. The transient solution accounts for the discharge rate controlled by the magma chamber volume, showing an intermediate-term periodicity of ca. 220 d, in agreement with the intermediate-term periodicity of ca. 146–292 d (ca. 5–10 months) obtained from the wavelet analysis (Fig. 3a).

Figure 4c reports a representative example of the solution obtained with the dual magma chamber model in order to assess the effect of the deep chamber on the discharge rate. We fixed the volumes of the deep and shallow magma chamber at 40 and 650 km³, respectively. The shallow dyke is 6500 m long with a width 2a=260 m and thickness 2b=4 m. The deep dyke has a width 2a_0d=500 m and a deep influx rate $Q_{\mathrm{in},\mathrm{d}}=\mathrm{2.3}$ m³ s⁻¹. A cyclic behavior of ca. 825 d is observed, reaching a peak discharge rate of ca. 6 m³ s⁻¹. This result is in agreement with the long-term periodicity of ca. 547–913 d (ca. 1.5–2.5 years) derived from the wavelet analysis (Fig. 3a).

Considering uncertainties in both modeling results and parameters and the fact that the thickness and width of the dykes are functions of the local overpressure, the results are quite consistent, although with a single model configuration the current approach cannot reproduce the periodicity observed at different timescales.

Figure 4 Results of numerical simulations. The physical framework of the conduit feeding system has deep and shallow chambers connected to the surface via vertical elastic dykes evolving into a nonelastic cylinder. The length of the shallow dyke L_ds is in the range of 6000–6500 m. The passage to the cylinder conduit L_T occurs at ca. 1300–500 m below the cone. (a) Discharge rates vs. time considering the elasticity of the shallower dyke, with a width 2a=400 m and thickness 2b=2 m. The cylinder diameter is D=30 m. Two cases are shown: (i) constant pressure (blue line) and (ii) constant influx rate at the source region of the dyke, providing different periodicities of 16 and 40 d, in good agreement with the short-term (weekly) periodicities observed in Fig. 3a; (b) discharge rate vs. time using the single magma chamber model. The dyke width 2a=600 m and thickness 2b=4 m. The chamber has a volume V_ch=30 km³, receiving a constant influx $Q_{\mathrm{in},\mathrm{s}}=\mathrm{2.3}$ m³ s⁻¹; periodicity is ca. 220 d, in good agreement with the intermediate-term (monthly) periodicities observed in Fig. 3a; (c) discharge rate vs. time using the dual magma chamber model. The aspect ratios of the shallow and deep chambers (AR_s−AR_d) are equal to 1.3 and 1.4, respectively. The upper feeding system has a chamber (V_chs=30 km³) connected to a dyke (width 2a=260 m; 2b=4 m) evolving into a cylinder (D=30 m) at L_T=1000 m. The shallow chamber is connected to the deep one (V_chd=500 km³) through a feeder dyke (2a_0d=500 m). A constant $Q_{\mathrm{in},\mathrm{d}}=\mathrm{2.3}$ m³ s⁻¹ is injected from below. Periodicity is in the range of ca. 825 d, in good agreement with the long-term (yearly) periodicities observed in Fig. 3a.

6 Discussion

In recent years, many studies have focused on magma flow dynamics in volcanic conduits during lava-dome-building eruptions (Melnik and Sparks, 1999; Wylie et al., 1999; Barmin et al., 2002; Melnik and Sparks, 2005; Costa et al., 2007a, b; Nakanishi and Koyaguchi, 2008; Kozono and Koyaguchi, 2009, 2012), highlighting periodic variations in the discharge rate due to the transition from a low regime (allowing for efficient crystal growth and leading to an increase in magma viscosity) to a high regime (with negligible crystallization). This difference in discharge rates can be of orders of magnitude, with strongly nonlinear responses to the variation of governing parameters from the volcanic system.

This behavior allows for periodic oscillations of the discharge rate (Nakada et al., 1999; Watts et al., 2002), as observed in different dome extrusion eruptions (e.g., Mt. St. Helens, Santiaguito, Soufrière Hills; Melnik et al., 2009). Although each volcano usually shows its complex pattern of discharge fluctuations, the cause can be explained as the superimposition of long-, intermediate-, and short-term effects of the coupled magma chamber(s) and conduit dynamics.

The long-term oscillations in discharge rate are a function of magma chamber size, magma compressibility, the amount and frequency of magma recharge, and withdrawal (Barmin et al., 2002; Costa et al., 2007b, 2013; Melnik et al., 2009). The short-term and intermediate oscillation dynamics can also superimpose on the main long-term periodicity through small changes in magma temperature, water content, and the kinetics of crystallization during magma transit in the conduit (e.g., Melnik et al., 2009).

The aforementioned eruptive behavior also characterized the Fuego de Colima activity in the 1998–2018 period, as demonstrated by the wavelet analysis of satellite thermal data. It is important to stress that the oscillating behavior is not regular, having a period between 2007 and 2012 that does not show any significant periodicity (Fig. 3a), possibly indicating a damped oscillation (Appendix Fig. A2). During this period the volcano enters an almost quiescent state with very low discharge rates. This period of low discharge rates is punctuated by low explosive activity, triggered by dome collapse or pressurization of the upper conduit.

It is well known for Fuego de Colima that Vulcanian explosions can evacuate significant portions of the upper conduit and destroy the lava dome. The influence of these processes on the periodicity of at least short-term periodic regimes could be significant. However, it is expected that such processes should mainly affect sub-daily periodicities, as explained by Costa et al. (2012), who analyzed the periodicity variations due to the collapse of a 200 m high plug at Montserrat. These changes should also have significant effects on the multi-week periodicity analyzed here.

Certainly, we cannot exclude the possibility that an exceptionally large evacuation of the upper conduit would be able to influence longer periodicities such as those investigated here, causing a transition to a more explosive eruptive style (i.e., Plinian) (Massaro et al., 2018b).

In order to investigate the relationship between the periodic components observed in wavelet analysis and the dynamics of the Fuego de Colima feeding system, we run simulations using the numerical model of Melnik and Costa (2014) (Fig. 4). The model can reproduce the results of the wavelet analysis in terms of observed periodicities, allowing us to relate short-, intermediate-, and long-term oscillations in discharge rates to the dynamics of the upper conduit, shallow magma chamber, and coupled shallow and deep magma chambers, respectively. This implies that the pressurization of the deep magma chamber has cascade effects on the whole feeding system of the Fuego the Colima, similar to what has been observed in other recent lava dome eruptions (i.e., Montserrat, Melnik and Costa, 2014).

It is of particular interest that the best output with the dual magma chamber model indicates that chambers do not oscillate simultaneously (“decoupled oscillation”; Fig. 4c). Although the presented data provide, for the first time, a framework able to describe the periodic behavior of effusive activity at the Fuego de Colima volcano, both the numerical model and wavelet analysis suffer from some limitations that need to be taken into account in interpreting the results.

• i.

  The available data for discharge rates and dome volumes collected for the 1998–2018 period do not have the same quality. This led us to extract only averages of discharge rates for the entire period, with biasing effects to lower amplitudes.

• ii.

  A common weakness of the spectral and wavelet analysis techniques is their inability to distinguish the source of any given periodic component (i.e., whether it is a signal from a volcanic process, an external process, or noise in the dataset). Elucidating the exact mechanism requires competing robust models and multiple independent field observations (Odbert and Wadge, 2009).

• iii.

  Assumptions behind the numerical model imply several limitations, such as those due to the constant value of the dyke width and simplified Newtonian rheology.

The first assumption greatly oversimplifies the physics. In the case of large overpressures, stress at the dyke tips will exceed the fracture toughness of the rocks and the dyke will expand horizontally (Massaro et al., 2018a), reaching some equilibrium configuration. When the deep chamber deflates, overpressure in the deeper dyke will decrease and, as the flow rate decreases, magma at the dyke tips can solidify, leading to a decrease in 2a_0d (Kavanagh and Sparks, 2011; Melnik and Costa, 2014).

Thermal exchange with wall rock can also affect the nonlinear dynamics of the system (Costa and Macedonio, 2002; Melnik et al., 2009). In addition, a more realistic estimate of the magma viscosity during lava dome eruptions should account for the coupling with energy loss, viscous dissipation, and stick–slip effects (e.g., Costa and Macedonio, 2005; Costa et al., 2007c, 2013.

Although this study revealed that different periodic signals are controlled by different mechanisms occurring in the plumbing system, the current model approach is not able to describe the three periodicities (long, intermediate, and short term) using a unique model configuration. Nevertheless, we hope this work will motivate further numerical modeling approaches in order to develop more sophisticated models able to describe the three timescales together, by incorporating further physical aspects (e.g., full thermal effects) and considering 3-D geometries.

7 Conclusions

The coupling of wavelet analysis and numerical modeling allowed us to decipher the eruptive behavior of Fuego de Colima in the period 1998–2018, as revealed by satellite thermal data. Three periodicities emerged from the study: (i) long term, (ii) intermediate term, and (iii) short term. The long-term periodicity extracted from wavelet analysis is ca. 913–547 d (ca. 1.5–2.5 years). It was replicated by the dual magma chamber model, which provided a periodicity of ca. 1000–500 d. The intermediate-term periodicity obtained from wavelet analysis is ca. 146–292 d (ca. 5–10 months), fairly replicated by the single magma chamber model with a periodicity of ca. 220 d. The short-term periodicity of ca. 18–38 d (ca. 2.5–5 weeks) is matched by model outputs considering the dynamics of the upper conduit (ca. 16–40 d).

The depicted behavior of effusive activity at Fuego de Colima is presented here for the first time, showing how the volcano presents similarities to the eruptive dynamics of other recent lava dome eruptions (i.e., SHV, Montserrat; Costa et al., 2013).

Appendix A

Figure A1 Sensitivity tests for steady-state solutions of discharge rate vs. chamber pressure (top) and time evolution of discharge rates (bottom). These solutions refer to the following main input parameters: (i) dyke thickness 2b=40 m as the conduit diameter at the top (D=2b), the transition from the dyke to cylindrical conduit L_T=500 m below the surface, the length of the dyke L_d=6 km³, and the volume of the magma chamber V_ch=50 km³. (a) General solution showing the transient regime in which the periodicity can occur; (b) solutions influenced by the dyke width 2a (from 200 to 400 m); (c) solutions influenced by the proportion of the water content in the melt (H₂O from 4 % to 5 %).

Figure A2 Sensitivity tests for transient solutions using the single magma chamber model. As a reference these solutions have the same main input parameters used for Fig. A1. (a) Dependence of discharge rate on time considering the influence of the dyke width 2a from 200 to 400 m; (b) influence of the magma chamber volume V_ch (from 20 to 50 km³); (c) dependence of discharge rate on time considering the dyke elasticity. Each curve shows a solution with a constant influx rate Q_in (in the range of 0.01–0.1 m³ s⁻¹).

Figure A3 Sensitivity tests for transient solutions using the dual magma chamber model. The shallow feeding system has a dyke with a width 2a=200 m, 2b=40 m, and L_T=500 m. The cylindrical conduit diameter is D=2b. For each diagram, the outflow (Q_out; black red and green lines), the flux entering the shallower magma chamber (Q_ins; blue line), and periods in days (T) are indicated. Runs of section 1–2–3 have fixed $Q_{\mathrm{in},\mathrm{d}}=\mathrm{2.3}$ m³ s⁻¹. Section 1: the volumes of the shallow and deep magma chambers are fixed to 40 and 650 km³, respectively. A set of runs is carried out for three different aspect ratios (ARs) of the shallow and deep chambers (AR_s=1; AR_d=1, AR_s=2; AR_d=1, AR_s=2; AR_d=1.5) considering three widths of the deeper dyke (2a_0d=200 m – black line, 1000 m – red line, 3000 m – green line). Section 2: the volume of the deeper magma chamber and the aspect ratios of both shallow and deep chambers are fixed to 650 km³ and ${\mathrm{AR}}_{\mathrm{s}}={\mathrm{AR}}_{\mathrm{d}}=\mathrm{1}$. A set of runs is provided for three different shallow chamber volumes (V_chs=30, 40, 50 km³) considering three widths of the deeper dyke (2a_0d=200 m – black line, 1000 m – red line, 3000 m – green line). Section 3: the shallow chamber volume and the aspect ratios of both shallow and deep chambers are fixed to 40 km³ and ${\mathrm{AR}}_{\mathrm{s}}={\mathrm{AR}}_{\mathrm{d}}=\mathrm{1}$, respectively. A set of runs is carried out for three deep chamber volumes (V_chd=550, 650, 750 km³) considering three widths of the deeper dyke (2a_0d=200 m – black line, 1000 m – red line, 3000 m – green line). Section 4: the shallow and deep chamber volumes are fixed to 40 and 650 km³, respectively. Two set of runs are carried out for Q_in,d equal to 1 and 3 (m³ s⁻¹). The aspect ratios (ARs) of the shallow and deep chambers are both equal to 1, considering three widths of the deeper dyke (2a_0d=200 m – black line, 1000 m – red line, 3000 m – green line).

Figure A4 Results of the Fourier analysis. (a) The 2002–2006 period shows two main periodic components, T₀=24.70 and T₁=6.17 months, corresponding to ca. 2 years and ca. 6 months, respectively. (b) The 2013–2016 period shows similar results: T₀=24.94 and T₁=6.23 months, corresponding to ca. 2.1 years and ca. 6 months, respectively.