Near-Real time analysis of seismic data of active volcanoes : Software implementations of time sequence data analysis

This paper presents the development and applications of a software-based quality control system that monitors volcano activity in near-real time. On the premise that external seismic manifestations provide information directly related to the internal status of a volcano, here we analyzed variations in background seismic noise. By continuous analysis of variations in seismic waveforms, we detected clear indications of changes in the internal status. The application of this method to data recorded in Villarrica (Chile) and Tungurahua (Ecuador) volcanoes demonstrates that it is suitable to be used as a forecasting tool. A recent application of this developed software-based quality control to the real-time monitoring of Teide – Pico Viejo volcanic complex (Spain) anticipated external episodes of volcanic activity, thus corroborating the advantages and capacity of the methodology when implemented as an automatic real-time procedure.


Introduction
One of the main goals of volcanologists is to forecast episodes of major activity.This study addresses the capacity of real-time (RT) monitoring of seismic data to provide information on which to base forecasts of volcanic activity.
Nowadays, seismic acquisition systems record data in RT and in continuous mode at high sampling rates, thus providing data instantaneously.High sampling rate measurements contain information that can be extracted at a com-Correspondence to: J. Vila (jvila@am.ub.es) parable rate to other parameters.Some of the most relevant examples of these measurements are the Real-time Seismic Amplitude Measurement (RSAM) (Endo and Murray, 1991) and the Real-time Seismic Spectral Amplitude Measurement (SSAM) (Rogers and Stephens, 1995).These techniques offer the advantage that the time series obtained have lower sampling rates that can be readily compared with nonseismic data, thus simplifying and accelerating interpretation.
RT computation has proved to be an excellent tool to obtain data to forecast volcanic activity.In this regard, RSAM techniques are widely used in many volcanoes and have shown excellent performance (Sparks, 2003).At present, seismic data acquisition systems are based mainly on software programmes for computers and their CPU clock allows time to perform RT analysis and data acquisition simultaneously.
In the present paper we describe the adaptation of a software-based quality control system to monitor the internal status of a volcano.This tool was developed initially in 2002 to test permanent stations at the "Observatori Fabra" in Barcelona (Llobet et al., 2003).The extension of this system to large networks was accomplished in 2004 and 2005 and it is currently running as a quality control device in several seismic observatories (Sleeman and Vila, 2007).In the following sections, we describe a number of examples of the kind of signal variation that can be detected by this monitoring system and the results of its application to data obtained in Villarrica (Chile), Tungurahua (Ecuador) and Teide -Pico Viejo complex (Spain) volcanoes.

Methodology
Our methodology is based on a seismic data quality control (QC) procedure that obtains measurements of several parameters for which alterations can be attributed to variation of incoming signals.Problems related to full system intrinsic anomalies (e.g.malfunctioning, changes in the installation), influences of new incoming signals or simply corrupt data are signal characteristics than can be monitored and identified.This methodology can also be an excellent complement of any implemented state-of-health channels.Here we present a proposal of a preliminary analysis that selects fixed length segments from the already built basic database format, assigns time on the basis of each segment and applies various data analysis subroutines.
The QC implementation subroutines compute parameters that involve both time and frequency domain values.After the extraction (or collection) of raw data from the pool of continuous data into segments of a fixed length, average amplitude and average square amplitude in the time domain are obtained for each segment.The frequency domain procedure of the analysis uses a Power Spectral Density (PSD) estimation of the ground acceleration.For each segment, the PSD is estimated using periodogram averaging (Welch, 1967).Only positive frequencies are taken into account (so-called onesided PSD) to compare the PSD with the USGS's High Noise Model (NLNM) and Low Noise Model (NLNM) (Peterson, 1993).PSD values are smoothed slightly by taking the average PSD values in a constant relative band width of 1/10 of a decade.The PSD is deconvolved with the instrument response to convert the PSD from digital counts [counts 2 /Hz] into ground acceleration units [m 2 /s 4 /Hz].
The PSD of each segment is integrated in diverse frequency bands to give information on how the energy is distributed in the spectrum and it is stored as function of time.For selected frequencies the PSD is also stored as function of time.Both the selected frequencies and frequency bands can be adapted to the specific requirements of the study.Finally, the bottom envelope of PSD is tracked and updated after the processing of a new segment, thus giving an estimation of all non-transient signals present during the period of operation.This envelope is stored after one full day of operation, shows the lowest noise measured at each frequency and is a reference of the minimum levels of daily activity (Vila et al., 2006;Sleeman and Vila, 2007).
The newly created time series (selected frequencies and PSD integrations as function of time) have a sampling rate established by the length of the time interval considered as input data.These time series are suitable to be used as new input for any supplementary near real-time (NRT) analysis and to be compared with state-of-health, weather/environmental data or data from other surveillance systems, such as deformation, temperature and gas emission, for multi-disciplinary analysis.The lower sampling rates of the new time series also allow fast and easier representation, thereby facilitating visual inspection and interpretations, and thus approaching the third level of observatory automatic procedures (ESF-EVOP Working Group, 1994).Examples and applications of RT QC in seismic networks of diverse nature can be found in the web pages -http://www.orfeus-eu.org/data-info/dataquality.htm,-http://bbnet.gein.noa.gr/QC/htdocs/index.html, -http://sismic.iec.cat/stations/qcm.html.

Application to volcano monitoring
Although volcanology involves many disciplines, the study of seismic activity is still the most used surveillance methodology.Rock fracturing, magma transfer, bubble interactions, etc., cause seismic signals that can be monitored and analyzed in RT.The classical methods to study seismicity associated with volcanic phenomena address discrete information (events) and its classification in categories (Minakami, 1974).However, not only the transitory seismic signals provide insights on the internal status of a volcano, but also the continuously generated background noise Seismic phenomena usually generate signals with well characterized frequencies.The behavior of these frequencies may provide early warning of internal changes, for instance, occurrence and/or shifts of predominant frequencies.Many times signs of activity reported in frequency domain can be inappreciable in a time domain analysis.The most adequate frequency components and frequency bands to be monitored may depend on particular situations that are not evident until the data analysis is in progress.
In the following section we present two "a posteriori" examples (Villarrica, Chile and Tungurahua, Ecuador) in which the potential capacity of QC implementation is shown and its application to Teide -Pico Viejo volcanic complex (Spain) in RT is described.

Villarrica volcano
Figure 1 displays the application of a QC method to data recorded in Villarrica volcano (Chile) in 2000.Plot (a) shows the evolution of the amplitudes of several selected frequencies, and plot (b) the evolution of the frequencies in which the maximum occurs in bands 1-5 Hz and 5-10 Hz.Plot (c) displays the inverse of the integrated square amplitude in time domain and plot (d) the evolution of the maximum spectral amplitude in the frequency band 1-10 Hz.All graphs display changes in behavior before an explosive episode as a consequence of plug removal (arrow 1) and as a culmination of increased internal activity (arrow 2), as reported by Ortiz et al. (2003) and Tárraga et al. (2006a).
The outputs related to amplitudes (Fig. 1a, c and d) show continuous variations, with the changes on the slopes being potential indicators of variations in internal activity that may produce an explosion episode (arrows).As presented in Ortiz et al. (2003), the eruption was forecasted by means of a material failure forecasting method (FFM) (Voight, 1988; De la Cruz-Reyna and Reyes-Dávila, 2001).
In the case of outputs related to predominant frequencies (Fig. 1b), the analysis also shows that changes are displayed as sudden shifts on the plots.Before the reported eruptive episodes, clear changes are observed.In a RT analysis these changes may be useful to provide a first-order alarm.
Inspection of the behavior of predominant frequencies in diverse frequency bands can determine the frequency range that is most sensitive to variations in volcanic activity in each stage.Once determined, a focused analysis of these sensitive frequencies can be performed.Apart from the general trends of precursory activity, Fig. 1b reveals further significant information.Large dispersions in the frequencies at which the maximum occurs imply low differences between values of amplitude in the selected frequency interval.Low dispersions indicate that an evident source dominates the frequency band, independently of its absolute amplitude.Before the moment in time marked by arrows (1) and (2), Fig. 1b displays incoming new and well defined input signals.This is not only seen in the frequency band 1-5 Hz, but also in the 5-10 Hz band.In the case of the second episode, this new incoming signal appears more than two days before its detection in the 1-5 Hz band.
In Fig. 1b the shift of the predominant frequency between 1.0 Hz and 5.0 Hz (time axis 299.6) coincides with the inflexion of the maximum amplitude in the same frequency band, whereas the shift that appears in the 5-10 Hz band coincides with the change in the behavior of Fig. 1c, which is related to the energy released.

Tungurahua volcano
The 2003 annual report of Tungurahua volcano (IGEPN, 2004) describes that the seismic activity before 20 August 2003 was low and the unrest activity of the volcano was classified as of low energy.Figure 2 shows the seismic record of 20 August 2003, in which a tectonic earthquake (A) of magnitude 4.5 located close to the volcanic edifice is registered.Seven hours later (C), an increase in volcanic activity started and it lasted for several days.
Figure 3 displays the evolution of the unrest activity of Tungurahua from 20 to 25 August 2003, on the basis of the methodology presented in this paper.Figure 3a and b represent the evolution of the maximum amplitude and the frequency at which the maximum occurs in the frequency bands 1-5 Hz and 5-10 Hz respectively.Figure 3c shows the integrated square amplitude and the integrated spectral square amplitude in the frequency band 1-10 Hz.
The behavior described in Fig. 3 (Tungurahua volcano) differs significantly from that observed for Villarrica volcano.In the present case, the perturbation produced by the tectonic earthquake modifies the "stationary" unrest.This is detected mainly in the seismic amplitudes both in time and in frequency domain but the frequency of the maximum peak is not sensitive to the increase in volcanic activity, Fig. 3b.However, Fig. 3a and c indicate a change in behavior, which can not be identified by visual inspection of the daily seismogram in time domain presented in Fig. 2. Figure 3 indicates that 5 h after the tectonic earthquake (arrow A), the parameters related to seismic energy start increasing continuously (arrow B) and the change in the behavior can be appreciated in the seismic record 2 h later (arrow C).This is a clear example of a volcanic episode triggered by an external mechanism (Tárraga, 2007) and can be considered similar to those that characterized many other eruptions or reactivations related to seismic interactions (see Hill et al., 2002;Manga and Brodsky, 2006).Although this volcanic episode was the result of an external influence, the results shown in Fig. 3 demonstrate that the reactivation can be monitored with the methodology presented in this paper.
Upon observation of the amplitudes of the seismic signals, both in spectral and in time domain, Fig. 3a and c, subsequent oscillations can be appreciated.These could indicate that the system is returning to levels of internal activity close to those present before 20 August.However, Fig. 3b shows a new predominant frequency around 8.5 Hz clearly visible on 23 August.This frequency persists noticeably while amplitudes appear to return to their initial stage.The presence of new frequencies may indicate that the system is being affected by another phenomenon.According to IGEPN (2004), the biggest explosions of the cycle occurred on 25 August, just after the time interval shown in Fig. 3.

Teide -Pico Viejo volcanic complex
After more than 30 years of seismic and volcanic quiescence, at the end of 2003 the region around the Canary Islands (Spain) started to show signs of seismic and volcanic activity (García et al., 2006), in particular in the northern part of the island of Tenerife.Moreover, the seismic noise recorded in Tenerife since 2004 contains signals of possible volcanic origin (Carniel et al., 2008a(Carniel et al., , 2008b;;Tárraga et al., 2006b;Tárraga, 2007) and a clear interaction between seismic events and the continuous background noise seismic signals has been reported (Tárraga et al., 2006b).These observations provided the starting point of a project to detect any sign that could provide an indication of an early stage volcanic episode.Given the increase in seismic activity observed in Tenerife, the first part of the study addressed the development and implementation of a seismic monitoring system.
The experience and promising results obtained from the "a posteriori" test performed in Villarrica and Tungurahua volcanoes led to the installation of a RT QC system in Teide -Pico Viejo volcanic complex.This system is designed to monitor changes in the behavior of the seismic signals in order to determine whether changes in the background noise status reveal variations in local response caused by alterations in volcanic activity.
The particular adaptation of the methodology for the Canary Islands (named Teide Information Seismic Server -TISS) involves the use of incoming data in 60-min long segments.The analysis and generation of output time series is only delayed by the completeness of the segment to be analyzed.The bottom envelope of the PSD is tracked, updated daily after processing 24-h intervals, and compared with previous days.The system has been monitoring the activity of Teide -Pico Viejo volcanic complex since June 2004 and has proven very useful to track changes in the behavior of processes (García et al., 2006).
Figure 4a displays the evolution of the integrated PSD amplitude and Fig. 4b the evolution of the frequency at which the maximum occurs in diverse frequency bands, together with the occurrence of seismic events, as reported by the Spanish National Agency IGN (http://www.geo.ign.es).A relevant example of the utility of RT monitoring can be appreciated from the output series of the TISS on 29 November 2004, which shows a systematic increase of energy in every frequency band.On 5 December 2004, a fumarole emission (largest green arrow 1) appeared in La Orotava valley (northeastern flank of Teide -Pico Viejo volcanic complex) and several local seismic events with magnitudes from 1.4 to 2.0 were registered during the following days.The changes in the internal status of the volcanic system (smaller blue arrows) are identified in Fig. 4b.Thanks to the RT monitoring, warnings (red medium size arrows) were issued in advance of the occurrence of the two phenomena.

Discussion and conclusions
A volcanic eruption occurs when magma or gas, which have migrated through the lithosphere and the volcanic edifice, are released from the earth's crust.Before the external surface manifestation, this matter and heat transfer produce physical and chemical phenomena that can be detected by instruments.These are precisely the signals that contribute to forecasting eruptions.
QC procedures to identify changes in seismic behavior of active volcanoes provide an excellent tool to complement other surveillance techniques based on diverse monitoring methodologies.The off-line application of a QC system to data recorded in Villarrica and in Tungurahua volcanoes demonstrated that changes in seismic behavior show good correspondence to eruptive episodes.After the reactivation of Teide -Pico Viejo volcanic complex, a RT monitoring system was set up with the aim to collect results almost simultaneously with the data acquisition processes.
Preliminary results of its application in RT demonstrate the usefulness of the QC methodology to provide indications of changes in the internal status of a volcano.It is worth mentioning that the current activity status of Teide -Pico Viejo volcanic complex can not be considered as high as that of Tungurahua or Villarrica.Thus the detection capacity of the QC procedures based on the analysis of background noise at its lowest levels is accentuated.
On the basis of the results obtained in the present study, we propose that frequency be considered one of the most significant indicators of changes in dynamical systems, as shown in Figs.1b, 3b and 4b.Analysis of the shifts of predominant frequencies provides unambiguous information and may indicate the start of internal phenomena, as previously described in seismic quality control systems (Llobet et al., 2003) and also reported in Rogers and Stephens (1995).
Here we have reported a close association between changes in the behavior of the spectral amplitude and volcanic unrest, see e.g.Figs. 1d,3a,and 4a.The conjunction of several of the outputs displayed in Figs. 1, 3 and 4, and a selection of the most appropriate channels may result in a considerable improvement in the confidence levels of RT analysis for volcano surveillance purposes.
Our study does not aim to determine any algorithm to be applied to volcanoes.Here we present an applied methodology that, from a generic point of view, is a potentially useful tool for monitoring the status of a seismic station.From a scientific perspective, this tool can provide accurate data and characterization of the behavior of a seismic site.From a technical point of view, the RT analysis is an excellent QC procedure to detect any type of variation in seismic activity, including internal malfunctioning or site modifications that affect the waveform.A single artifact may be enough to indicate emergency situations when anomalous behavior is detected or when a dynamical system is stable or not.QC subroutines, together with modern data transmission procedures allow the collection of RT information without the complex logistics required by permanent observatory infrastructures.

Fig. 1 .
Fig. 1.Villarrica volcano: Evolution of several output parameters provided by the QC system.The two arrows indicate the end of reported eruptive episodes.Time axis is common to all four plots (see text for details).

Fig. 3 .
Fig. 3. Tungurahua volcano: Evolution of several output parameters provided by the QC system.Arrows A, B and C coincide with A, B and C in Fig. 2. Time axis is common to all three plots (see text for details).

Fig. 4 .
Fig. 4. Time evolution of the integrated spectral square amplitude (a) and predominant frequency (b) in several frequency bands, as provided by the RT QC system specifically developed to monitor the status of Teide -Pico Viejo volcanic complex.Vertical lines indicate the occurrence of seismic events, represented by circled crosses (see text for details).