Detecting multiscale periodicity from the secular effusive activity at Santiaguito lava dome complex (Guatemala)

Santiaguito, Guatemala, represents one of the best cases of active lava dome complex in the world, producing lava flow effusion, weak explosive activity, and cycles of lava dome extrusion over varying timescales. Since the inception in 1922, it has shown a remarkable constant eruptive activity, characterized by effusion of blocky domes and lava flows punctuated by moderate explosions of gas-and-ash and pyroclastic flows. In this study, we reconstruct the time evolution of discharge rates of Santiaguito across one entire century, from 1922 to 2021, combining, for the more recent activity, new satellite thermal data. By using discrete Fourier transform (DFT) and Morlet wavelet analyses, we identify three fundamental periodicities in subsets of the 1922–2021 time-series: (i) long term (ca. 10 years), (ii) intermediate term (ca. 3.5 years), and (iii) short term (from ca. 1 year to ca. 3 months), which are comparable with those observed at other lava dome eruptions at calc-alkaline volcanoes. Such inferred periodicities provide a powerful tool for the interpretation of the non-linear eruptive behaviour and represent a pivotal benchmark for numerical modelling aimed to reconstruct the dynamics of the magma feeding system based on a time-averaged discharge rate dataset.

Changes in eruption dynamics can be driven by complex variations in magma viscosity, pressurization of the feeding system, or changes of local stress loading. Changes in viscosity are governed by many physical factors as temperature (e.g., Hess et al. 1996), dissolved volatile content (Richet et al. 1996), distribution and shape of bubbles (e.g., Llewellin et al. 2002;Pal 2003) and kinetics of crystallization (e.g., Melnik and Sparks 1999;Barmin et al. 2002;Kozono and Koyaguchi 2012). Pressurization of the whole magmatic feeding system (shallow and deeper magma reservoirs plus conduit) due to injection of fresh magma has been claimed at different volcanoes (e.g., Sparks 1999, 2005;Costa et al. 2007a,b;Koyaguchi 2010, 2012), especially for long-term variation of eruptive style (e.g., Sulpizio et al. 2016). Changes in local stress due to tectonics, gravitational spreading of the volcanic edifice, or gravitational unloading may in some cases trigger sudden transitions from effusive to explosive eruptions (e.g., Korringa 1973;Costa et al. 2011;Sulpizio and Massaro 2017;Massaro et al. 2019).
The extrusion rate observed at lava domes may be complex, with sudden changes in magma discharge rate that cannot be assigned to any well-defined regular pattern (e.g., Lascar volcano, Chile, Matthews et al. 1997). On the other hand, discharge rates may fluctuate regularly, reaching maximum values of the order of 20-40 m 3 s −1 , but usually about 10 m 3 s −1 (Sparks 1997; Melnik et al. 2008). Periods of slow extrusion or even pause are common (e.g., 1980-1982on Mount St. Helens, Swanson and Holcomb 19901998-2019activity at Fuego de Colima, Mueller et al. 2013Cassidy et al. 2015;Massaro et al. 2019).
In particular, the temporal changes in discharge rate have been studied on different timescales through analytical (e.g., Barmin et al. 2002;Gonnermann and Manga 2003) and numerical models (e.g., Sparks 1999, 2005;Costa et al. 2007a,b;Melnik and Costa 2014), showing superimposed short-, intermediate-, and longterm periodicities (e.g., Melnik et al. 2008;Costa et al. 2012;Massaro et al. 2019) which are generally explained by the upper conduit pressurization related to the nonconstant ascent of magma flow (Denlinger and Hoblitt 1999;Melnik and Sparks 1999;Voight et al. 1999;Wylie et al. 1999;Ozerov et al. 2003;Costa et al. 2007aCosta et al. ,b, 2012Kozono and Koyaguchi 2012) or controlled by pressure variations in magma reservoirs (Costa et al. 2007a;Melnik and Costa 2014). This means that the magmatic feeding system not only represents the transient path for magmatic fluids towards the surface, but it is also the theatre of fundamental processes influencing the style of activity (Spina et al. 2017).
The investigation on the cyclic activity at lava dome volcanoes has been improved in the last decades, with the availability of satellite thermal data (e.g., Harris et al. 2003;Coppola et al. 2009Coppola et al. , 2013Ramsey and Harris 2013) that can be studied using spectral analyses (e.g., Odbert and Wadge 2009;Odbert et al. 2014). Using data acquired by the space-based moderate resolution imaging spectroradiometer (MODIS) sensor, Wright et al. (2015) estimated the total amount of energy radiated into the atmosphere by 95 volcanoes that erupted between 2000-2014 suggesting how these volcanoes differ with regard to the amount of energy they radiate and how the total amount of energy released varies with eruption style, specifically whether the volcano is characterized by the emplacement of mafic lava flows or felsic lava domes, the presence of persistently active mafic lava lakes, or ventconfined explosive activity.
In this study, thermal remote sensing data along with decades of direct measurements, are used to detect the fluctuations in the magma discharge rate at the Santiaguito lava dome complex (Santa Maria volcano, Guatemala) during the last century . Santiaguito represents a unique case-study for deciphering periodicity in lava dome extrusion rate because of the available long period of historical (Harris et al. 2003) and satellite thermal data (Coppola et al. 2009(Coppola et al. , 2013. Magmatological data are abundant (Sahetapy-Engel et al. 2008;Scott et al. 2012;Singer et al. 2011;Cisneros de Lèon et al. 2019;Wallace et al. 2020), but the knowledge about the volcano feeding system is still poor due to the lack of accurate geophysical investigations (i.e. tomographic images). The aim of this work is to decipher the multi-term periodicity from the secular dome extrusion time-series at Santiaguito by using spectral techniques (i.e. discrete Fourier transform and Morlet wavelet analysis), with particular focus on the recent period  for which satellite thermal data are available and can be integrated with published data (Harris et al. 2003;Melnik et al. 2008;Scott et al. 2012;Drust, 2018;Rhodes et al. 2018;Zorn et al. 2020).

The historical activity of Santiaguito
Santiaguito is a ca. 1.1 km 3 dome complex belonging to Santa Maria volcano, located in Guatemala, in the Central American Volcanic Arc, which extends from Mexico to El Salvador (Fig. 1a;Rose 1973;Duffield et al. 1993;Bennati et al. 2011). The Santa Maria eruptive activity began around 103 ka, with four main stages before a 25-thousand-year period of quiescence prior to 1902 (e.g., Rhodes et al. 2018) when a devastating Plinian eruption occurred, erupting 8.5 km 3 of dacitic magma as fall deposits, and leaving an explosion/collapse crater on the south-west side of the edifice (Williams and Self 1983). In 1922, after a period of quiescence of 20 years, volcanic activity resumed with the extrusion of El Caliente dome in the crater excavated in 1902 that, without interruption, is still continuing today. Thus far, four domes have been constructed: El Caliente (1922-19391972-present day), La Mitad (1939-1949), El Monje (1949-1958), and El Brujo (1958-1986, extruding a combined total of ca. 1 km 3 of magma ( Fig. 1b; Harris et al. 2003). Following dome growth at El Caliente, La Mitad, El Monje, and El Brujo grew sequentially, moving progressively to the west until activity renewed at El Caliente coincident with activity at El Brujo in 1972 (Harris et al. 2003). Activity prior to 1929 was endogenous, characterized by subsurface build-up of magma to inflate and uplift the carapace of the El Caliente dome. Between 1929 and 1958 a combination of endogenous and exogenous behaviour occurred. Post-1958 growth has been solely exogenous, where lava extrudes at the surface (Rose, 1987;Harris et al. 2003). Since 1977, activity has been focused at El Caliente vent and consists of semi-continuous extrusion of blocky lava flows interspersed by frequent gas-andash explosions (Rose 1987;Harris et al. 2003). For the past two decades, explosions have generally been of small to moderate size with volatile-rich, ash-poor plumes Recently, the culmination of the inflation phases at Santiaguito results in either an outgassing event or a gasand-ash explosion, where the latter outcome is accompanied by inflation occurring at higher rates, as well as by the occurrence of a very long period seismic event. In December 2012, activity consisted of regular explosions (≤ 2 h −1 ) producing gas-and ash-plumes reaching heights of < 1 km De Angelis et al. 2016) and simultaneous lava flow extrusion accompanied by frequent incandescent rock falls.
A new lava dome appeared within the summit crater of Caliente in October 2016 and it has continued to grow, producing frequent block avalanches down the flanks.
During January 2018-September 2020, the activity remained centred at Caliente lava dome with the extrusion of blocky lava along with the emissions of gas and ash several times a day (Gottschämmer, et al. 2021). Daily explosions with ash plumes and block avalanches still continued during February-July 2021 (Global Volcanism Program 2021).
Since 1922, the volcanic activity showed variations in effusion rates (Rhodes et al. 2018). In particular, Harris et al. (2003) showed that lava extrusion is cyclic, with 3to 6-year-long episodes of high extrusion rate (0.5-2.1 m 3 s −1 ) followed by 3-to 11-year periods of low extrusion rate (≤ 0.2 m 3 s −1 ), estimating a time-averaged growth rate of 0.46 m 3 s −1 . More recently Zorn et al. (2020) measured extrusion rates of 0.04-0.06 m 3 s −1 , in accordance with previous estimations during low activity phases (Harris et al. 2003), while during high active phases extrusion rate values have been observed to exceed 2 m 3 s −1 . Other studies have found higher extrusion rates ( (Ebmeier et al. 2012). The extrusion rates at Santiaguito have then slowed down significantly, although activity was still intermittently higher during 2012 (Ebmeier et al. 2012).
According to Harris et al. (2003), nine different cycles have occurred with the latest emplacement of four lava flows > 2 km during 2011-2015, signalling the start of a new cycle. Rhodes et al. (2018) argued that a typical eruption cycle at Santiaguito begins with an increase in extrusion rate and the early extrusion of spines, followed by lava flows. Toward the end of a cycle the extrusion rate declines again, and the trend is reversed in spine growth. It has been observed that the longest lava flows were emplaced during constant, high-extrusion rates that can last up to 2 years (Harris et al. 2003). During periods of low extrusion rates, lava flows of decreasing length developed (Rhodes et al. 2018). Prior to 1990 the dacitic composition favoured conditions suitable for spine growth (Rhodes et al. 2018).
Visual observations between 2012 and early-2015 showed that the shallow dome structure in the El Caliente vent remained mostly intact during explosions, and lava flows descended the southeast flanks. Effusion rate declined through 2014, ceasing in December . From late 2015-2016, the long-term eruptive behaviour at Santiaguito underwent a transition to more powerful and less frequent explosions, producing plumes up to 8 km high that destroyed the summit dome and formed a deep crater in the El Caliente vent . From late-2016 onward, the activity has returned to the long-term behaviour with a new lava dome that slowly increased in size in 2017, generating constant steam and gas emissions along with daily explosions (Global Volcanism Program 2017). Similar activity was observed in 2018-2021 and is still in progress (Global Volcanism Program 2018, 2019. In recent years, geophysical monitoring campaigns have shown that the long-term, predominantly weak explosive activity is accompanied by repetitive "inflation-deflation" cycles as monitored by long-period seismic and tilt signals (e.g., Sanderson et al. 2010;Johnson et al. 2014;Scharff et al. 2014;Lamb et al. 2019). The culmination of the inflation phases at Santiaguito results in either an outgassing event or a gas-and-ash explosion (e.g., Ebmeier et al. 2012;Lavallée et al. 2015), where the latter outcome is accompanied by inflation occurring at higher rates, as well as by the occurrence of a very long period seismic event (Johnson et al. 2014). The shallow magmatic system thus resides in a delicate yet repeatable balance between pressurization, fracturing events, and outgassing regimes, where the rate and timescale of deformation ultimately characterize the eruptive style .

Data
In this work, we reconstructed the discharge rates from 1922 to 2021 using published data (Harris et al. 2003;Melnik et al. 2008;Scott et al. 2012;Drust et al. 2018;Rhodes et al. 2018;Zorn et al. 2020), integrated with satellite thermal data (middle infrared observation of volcanic activity-MIROVA system, 2022) available for the 2000-2021 period (Fig. 2). 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 (2022) data system. It focuses on the middle infrared (MIR) radiance to better detect and analyse thermal radiation emitted from volcanic sources, in terms of volcanic radiative power (VRP, Watt).
As in Massaro et al. (2019), the heat flux radiated from the Santiaguito dome complex was converted into averaged magma discharge rates to reconstruct the yearly, monthly, and weekly time-series from February 2000 to August 2021 (Fig. 2 b-d). The conversion was carried out considering the relationship between the volcanic radiated energy (VRE; time-integration of VRP through Fig. 2 a Averaged yearly discharge rates from 1922 to 2021 obtained combining literature data (Harris et al. 2003;Melnik et al. 2008;Scott et al. 2012;Drust, 2018;Rhodes et al. 2018;Zorn et al. 2020) and MODIS data (Coppola et al. 2016); b-c-d averaged yearly, monthly, and weekly discharge rates from 2000 to 2021 considering only MODIS data the eruption) and the erupted volume (Vol) by using an empirical parameter called radiant density (c rad = VRE/ Vol, J m −3 ; Coppola et al. 2013). Here we calibrated the c rad value using the estimates of the lava volume erupted between 2000 and 2009 calculated by Ebmeier et al. 2012 (120 × 10 6 m 3 ) and then compared to a VRE of 3.32 × 10 15 J related to the same period (derived from MIROVA data). The resulting c rad value, equal to 2.77 × 10 7 J m −3 , is compatible with lavas of acidic composition (Coppola et al. 2013) and it was therefore used to convert the VRP values into time-averaged lava discharge rate (TADR) values, through the relationship: TADR = VRP/c rad . The c rad value obtained for this conversion is the result of a complex interaction between the rheology of the lava body, the rate of cooling, topography, and conditions of emplacement. However, on the basis of numerous eruptions analysed, Coppola et al. (2013) estimated that for each eruptive case (such as the Santiaguito's case), the natural variation of this coefficient is limited, and produces an uncertainty on the single measure of effusion rate equal to ± 50%.
Overall, our data indicate that between 2000 and 2021 Santiaguito erupted 202 (± 101) × 10 6 m 3 of lava, with a mean output rate equal to 0.33 (± 0.16) m 3 s −1 (see Appendix 1). It is important to note that the instrumental limit of the MIROVA system is not able to detect thermal anomalies below 0.5-1 MW. Since we used a crad = 2.77 × 10 7 J m −3 , the minimum reliable value of the discharge rate is ca. 0.02 m 3 s −1 . Moreover, the thermal data obtained from MIROVA are not corrected for the presence and/or attenuation of clouds, therefore the effusion rates and volumes are to be considered minimum estimates (Coppola et al. 2016).

Spectral analyses
Time-series analyses are generally used to extract information and to identify patterns in fluctuations. In this study, we used both the discrete Fourier transform (DFT) and wavelet analysis to detect the multi-term cyclic behaviour during the last century of eruptive activity at Santiaguito lava dome complex.
Generally, the time-and frequency-domain approaches are two different ways of presenting the interaction between signals and systems. An arbitrary signal can be considered as a linear combination of frequency components, and the time-domain representation is the superposition sum of the frequency components. In this perspective, the DFT can calculate a signal's frequency spectrum. This is a direct examination of information encoded in the frequency, phase, and amplitude of the component sinusoids. Second, the DFT can find a system's frequency response from the system's impulse response, and vice versa. This allows systems to be analysed in the frequency domain, just as convolution allows systems to be analysed in the time domain. Third, the DFT can be used as an intermediate step in more elaborate signal processing techniques (e.g., Steven 2013). Before calculating DFT spectrum, we detrended the considered time-series using a smoothing procedure consisting of a running average of a time series, using a window size of 5.
Although Fourier analysis is well suited to the quantification of constant periodic components in a timeseries, it cannot recognize signals with time-variant frequency content. On the contrary, the wavelet analysis has the advantage to resolve the signal on a local time-scale domain by providing time and frequency, overcoming the problems of non-stationarity in timeseries (e.g., Lau and Weng, 1995;Torrence and Campo 1998). Using this approach, we can track how the different scales related to the periodic components of the signal change over time (Cazellas et al. 2008).
Following Torrence and Compo (1998) and Massaro et al. (2019), we applied the wavelet analysis on our yearly, monthly, and weekly averaged discharge rate time-series (Fig. 2). Each time-series is called xn, showing an equal time spacing δt and number of points n = 0… N − 1, respectively. Using the approximately orthogonal Morlet function wavelet ψ(η) (characterized by a zero mean and localized in both time and frequency space; Farge 1992; see Eq. 1), we have defined the wavelet transform Wn (s) as the convolution of x n with a scale s and the translated version of ψ 0 (η) called "mother function" (Torrence and Compo 1998: see Eq. 2). In formula, we have: where the (*) indicates the complex conjugate. The scale s should be equal to approximately according to the Nyquist theorem. Therefore, the smallest wavelet we could possibly resolve is 2δt, thus we choose as s 0 a period of 14 days.
It is possible to reconstruct the "local" wavelet power spectrum as the absolute value squared of the wavelet coefficients, |Wn (s)| 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.
Since the spectra generated from many geophysical systems are well described by red noise (e.g., Fougere 1985), we used the model for red noise given by (1) ψ 0 (η) = π − 1 4 e iω 0 η e −η 2 /2 , the univariate lag-1 autoregressive or Markov process (Torrence and Compo 1998) in order to determine the significance levels for each 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, which is defined as the probability that the true wavelet power at a certain time and scale lies within a certain interval of the estimated wavelet power (i.e. significance at the 5% level is equivalent to the 95% confidence level; Torrence and Compo 1998).
The criterion for applying the local wavelet analysis is very similar to those employed with classic spectral methods since it is a generalization of the Fourier transform. 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 (i.e. Torrence and Compo 1998;Cazelles et al. 2008;Massaro et al. 2019).

Results
In this section, the results obtained from the DFT and wavelet analyses are shown on the averaged discharge rates from 1922 to 2021 (Fig. 2). In Fig. 3, we reported the periodograms (or spectra) of the signal, considering different temporal ranges. Peaks indicate the main contribution in frequency to the variance of the series, from which periodicities, if present, can be detected. Fig 3a refers to the results of DFT analysis of the yearly averages of the whole discharge rate time-series from 1922 to 2021. It shows that the prominent frequencies vary from 0.023 to 0.102 year −1 , which are equivalent to periods in the range of 9.80-43.48 years. Fig 3b concerns the monthly averaged discharge rate time-series obtained with thermal data, from February 2000 to December 2021. In this case, we detected three main frequencies of ca. 0.025, 0.055, 0.070 and 0.084 month −1 corresponding to 3.33, 1.52, 1.19 and 0.99 years, respectively . Fig 3c is about the weekly averaged discharge rate time-series highlighting frequencies of ca. 0.019, 0.039 and 0.058 week −1 , corresponding to periods of ca. 1 year, 6 and 4 months. Figure 4 shows the normalized local wavelet power spectra of the signal from the 1922 to 2021. The colour scale for power values varies from light orange (low values) to dark red (high values). The thick black contours represent the 95% confidence level while the blue dotted line indicates the cone of influence that delimits the region not influenced by edge effects (Torrence and Compo, 1998). Different periodicities can be observed: (i) long-term periodicity of ca. 10 years, well defined from 1942 to 2002 (Fig. 4a); (ii) intermediate-term periodicity of ca. 3.5 year during 2007-2017 (Fig. 4b); and (iii) short-term periodicities of ca. 1-1.5 year during 2001-2003, 2005-2007 and 2011-2015, ca. 6 months during 2011-2014, ca. 3-4 months detected during 2001-2003, and more weakly in 2011-2016 (Fig. 4c). The 1-year periodicity is also detected in Fig. 4b between 2001-2003 and 2011-2014. It is worth noting that the red dashed lines drawn on the DFT spectra (Fig. 3) indicate the frequencies in agreement with the principal periodicities found in the normalized local wavelet power spectra of the signal, within the COI (Fig. 4).

Discussion and conclusions
Spectral analyses were used to detect multiscale periodicities in the secular Santiaguito's discharge rates, using satellite thermal data and published extrusion rates. The time-averaged discharge rates were analysed in the time and frequency domains using the wavelet transform (Morlet function; Torrence and Compo 1998). The wavelet power spectra were computed with a significance level of 95%, determined against a red noise background spectrum. Results point out three orders of periodicity in the 1922-2021 period: (i) a nearly constant periodicity of ca. 10 years over the whole investigated period (Fig. 4a); (ii) an intermediate periodicity of ca. 3.5 years, visible from 2007 to 2017 (Fig. 4b, c), and (iii) short-term periodicities of ca. 1 year, 6 months and 3-4 months detected in 2001-2003, 2005-2007, 2011-2015 (ca. 1 year), in 2011-2014 (ca. 6 months) and in 2001-2003 (only 3-4 months).
Multiscale periodicities have been described at other volcanoes (e.g., Soufrière Hills, Montserrat, Odbert and Wadge 2009;Odbert et al. 2014;Fuego de Colima, Massaro et al. 2019), and attributed to superimposition of the dynamics of deep and shallow components of the feeding system (e.g., Melnik and Costa 2014). Massaro et al. (2019) explained the multiscale periodicity of discharge (See figure on next page.) Fig. 3 The discrete Fourier transform analysis showing the a yearly (main frequencies at 0.023, 0.055, 0.078 and 0.102); b monthly (0.025, 0.055, 0.070, 0.084) and (c) weekly (0.019, 0.039, 0.058) power density spectrum of the discharge rate from 1922 to 2021. The equivalent periods for the detected frequencies are given in each plot. For the monthly and weekly time-series, the data cover the last 21 years . The red dashed lines indicate the frequencies showing the best agreement with the periodicities observed in the Morlet analysis (Fig. 4) Massaro et al. Earth, Planets and Space (2022)  Geophysical and petrological data of sufficient quality are not available for Santiaguito, and this prevents the constraint of any detailed numerical model for explaining the observed periodicities. Despite this, our spectral analysis may be used for gaining some hints on pressurization and depressurization cycles of the magmatic feeding system at Santiaguito in the 1922-2021 period. Similar to Colima and other volcanic systems, such as Soufriere Hills, Montserrat, the long-term (10 years) periodicity may be controlled by the dynamics of the deep feeding system, the intermediate (3-4 years) periodicity may be attributed to the pressurization/depressurization of a potential shallow magma chamber, and the short-term periodicities (3-4 months up to 1 year) to the dynamics of upper conduit system. We have to note that the 10-year periodicity is present over the entire investigated period, while the intermediate and short periodicities appear scattered during 2007-2017 (Fig. 4b), 2001-2003, 2005-2007 and 2011-2015 periods (Fig. 4c). This suggests that the deep feeding system was active almost constantly during the 1922-2021 period, whereas the contribution to the discharge rate of the intermediate and shallow parts of the feeding system occurs only sporadically. Unfortunately, this suggestion cannot be tested over the entire period, due to the limited time span covered by thermal satellite data.
These inferences are in good agreement with the results of Wallace et al. (2020), who explored the shift in eruption style at Santiaguito from early-2015 to late-2016, showing an escalation in explosivity which could have been influenced by the injection of a higher temperature, relatively volatile-rich magma into the feeding system. With respect to other studied volcanic systems (Melnik and Costa 2014;Massaro et al. 2019), the behaviour of Santiaguito lava dome complex can be even more complicated as the upper feeding system is connected to multiple vents. This is because clustered vents may reflect different conduits characterizing the shallow feeding system, which may contribute to dissipation of pressurization differently from a single conduit system.
In conclusion, discrete Fourier transform (DFT) and Morlet wavelet analyses of discharge rates at Santiaguito volcano for the 1922-2021 period were used to explore the non-linear eruptive behaviour of its lava dome complex at different timescales, providing a framework of the secular discharge rate fluctuations. Three orders of periodicities were detected: a nearly constant periodicity of ca. 10 years over the whole investigated period (Fig. 4a)