“ VOLCANIC AND VOLCANO-TECTONIC ACTIVITY FORECASTING: A REVIEW ON SEISMIC APPROACHES „

Forecasting volcanic activity is a difﬁcult problem that is being addressed worldwide from different perspectives. Signiﬁcant advances have been made after the introduction of non-linear dynamical systems theory and the use of power-law distributions of different geophysical parameters in the Earth Sciences. In particular, frequency-magnitude power-law statistics evidences the scale-invariance and self-organi-zation of seismicity, and brittle fracture models show that under certain conditions, a precursory causal evolution characterized by an accelerating strain rate culminates in a catastrophic failure of the system under stress. The precursory organization of the seismicity and the distinct characteristics of the seismic events have allowed the development of forecasting tools. In this work we present some examples of forecasting methods based on seismic observations at different volcanoes in the world, and how results and experiences has been used to improve both hardware and software tools developed for short-term forecasting of volcanic and seismo-volcanic activity.


INTRODUCTION
Forecasting volcanic eruptions is one of the most controversial issues in volcanology. Historically there are evidences that people abandoned hazardous areas before eruptions, even before the existence of modern instrumental monitoring methods. For example, the excavations conducted by Marinatos since 1967 in the Greek island Santorini found that little or no loss of lives occurred in the settlement of Akrotiri during the Theran eruption that destroyed the city around 1627 BC, thus suggesting that the town was evacuated be-fore the eruption [Pellegrino, 1991]. Similarly, the relatively small number of corpses found in Herculanum, destroyed by the Plinian eruption of Mt. Vesuvius in 79 AC, suggests a similar situation [Scandone and Giacomelli, 2014]. A move to prevent a tragedy such as abandoning endangered settlements depends on three basic factors: the capacity of observation to recognize significant changes in the patterns of volcanic activity, the perception that eruptive activity may cause fatal damage, and the capability of making a political decision resulting in a complex move as in the case of the evacuation of a community. Failure in any of these three factors may lead to a complete disaster, as it was the case in St. Pierre (Martinique) in 1902 [Tanguy, 1994], or in Armero (Colombia) in 1985 [Voight, 1990]. These three basic factors have been addressed for a long time in Japan. For example, the 1910 Usu volcano eruption was studied by scientists at that time, in particular by Omori and Satō [Omori, 1903[Omori, , 1911[Omori, , 1912Satō, 1913;Yokoyama and Matsushima, 2018], and the management of the emergency situation was properly handled. With this background Noguchi and Kamiya [1963], Mogi [1963] and Minakami [1964] purported the concept of volcano observatory as a mean for forecasting volcanic eruptions. Similarly, the analysis of the 1975 eruption of Tolbachik [Fedotov and Markhinin, 1983], the successive eruptions of Etna [e.g. Vinciguerra et al., 1999] and the 1980 eruption of Mt. St. Helens [Lipman and Mullineaux, 1981], provided evidence in setting some physical basis on the nature of eruptions and their precursors [De Natale et al., 1999]. An important boost of word-wide research works on forecasting and modeling of pre-eruptive processes from different perspectives followed afterwards.
The theory of nonlinear dynamic systems and the power-law distribution attributed to the energy released by different geophysical phenomena have greatly helped to deal with and model the complex nature of forecasting hazardous activity [e.g. Pshenichny et al., 2003Pshenichny et al., , 2009. It is very important to emphasize that the nature of the causal forecasting methods based on the concatenated increase of a geophysical observable implies dynamic time scales that decrease as the monitored system evolves (hours-days). In contrast, non-concatenated causal phenomena only provide "stationary" information on the state of a geophysical system which can only be used for the long-term (years-decades) assessment of the related hazard [Jones, 1996]. In a volcanic region, this may be exemplified as the difference between an accelerated release of seismic energy, probably leading to some type of volcanic or seismo-volcanic manifestation, and a volcanic hazards map.
In this work we aim to describe the evolution of the precursory unrest activity of active volcanoes and its relevance in forecasting eruptions through causal forecasting models and methods based on changes in the seismicity observed. Our experience is focused mainly on the measurement and data processing of volcano seismicity, and to a lesser extent, in real-time defor-mation measurement using GPS L1 and L2. In addition, we believe that these seismic models and methods are adequate for volcanoes where the assessment of new unrest detected with provisional monitoring equipment are critical for risk management, but due to the lack of resources, the deployment and maintenance of a comprehensive volcano monitoring network is not always possible.

METHODS BASED ON CHANGES IN SEISMI-CITY
As stated in the previous section, historical evidence indicates that the perception of a relationship between seismic and eruptive activity may be as old as humanity. However, in the XX century the advent of electronics, computers, more sensitive motion transducers, and the capabilities to program elaborate signal analysis tools, allowed the recognition of multiple types of signals produced by internal processes of active volcanoes, most of them with precursory value [Omori 1911[Omori , 1912Schick and Riuscetti, 1973;Schick 1981]. For example, during the 1982-1984 Campi Flegrei bradyseismic crisis [Barberi et al., 1984], the introduction of real-time digital instruments [Ortiz et al., 1984;1992] not only permitted a follow-up of the ongoing activity, but also initiated the development of hardware and software specifically designed for volcano monitoring that is still used [Barberi et al., 2008].
In general, seismic forecasting methods may be grouped in two types: those based on the time evolution of volcano-tectonic earthquakes, and those based on the analysis of waveform signals of volcanic origin.

PRECURSORY TIME EVOLUTION OF VOLCANO TECTONIC EARTHQUAKES
The rapid increase of the rate of volcano tectonic (VT) earthquakes before an eruption is one of the earliest manifestations recognized as an eruption precursor [Benoit and McNutt, 1996]. However, since many earthquake swarms in volcanic areas are not followed by eruptive activity, the causality between these signals and eruptions is weak [Gudmundsson, 2002;García et al., 2014a]. Migration of seismic foci has also been used since the 1960's as a forecasting technique but sometimes eruptions are not preceded by migrations and vice versa [Wiemer et al., 1998]. The time evolution of the seismic energy release by VT earthquakes and par-ticularly its cumulative value provides a more reliable tool for the eruption forecasting. Yokoyama [1988] found that VT cumulative seismic energy reaching a threshold of 10 11 J has significant forecasting value. It is also a relatively simple observable to evaluate because it requires only a basic seismic catalog with a magnitude of completeness of the catalogue Mc=2 or lower [García et al., 2014a]. Examples of such behavior are shown in Figure 1, where a clear increase of the seismic activity near the Teide volcano (Tenerife, Canary islands, Spain) in 2004 contrasts with the uniformity of the seismicity at regional level. Similarly, Figure 2 shows an example of how different seismic parameters evolved before and after an eruption during the El Hierro 2011 precursory unrest. The daily data processing allowed forecasting the eruption and the higher magnitude VT events [García et al., 2014a,b; using some of the methods described below. The daily data processing using different methods of analysis ((Mean Recurrence Time MRT, García et al. 2016, Real-time Seismic Amplitude Measurement RSAM, Endo andMurray, 1991, etc.) allowed to forecast the volcanic activity and pro-vide the authorities and relevant bodies with pertinent technical input [García et al., 2014b]. We implemented a variety of methods because the appropriateness (and inherent precision) of a given method depends on the type of changes in the volcanic behavior. The MRT based in the seismic catalogue is the most used method.

MAGNITUDE DISTRIBUTION
The distribution of earthquake magnitudes in a volcanic region may also be described by the frequencymagnitude distribution of tectonic earthquakes based on the Ishimoto and Iida [1939] and the Gutenberg-Richter [1944] law where N is the cumulative number of earthquakes with magnitude ≥M, and a and b are constants that describe the power law decay of occurrence frequency with increasing magnitude. The b-value is the slope of the line which best fits the earthquake-magnitude distribution data above the magnitude of completeness Mc, defined 3 FORECASTING VOLCANIC ACTIVITY FIGURE 1. Near and far field seismicity at Tenerife island during the 2004 Teide volcano seismic crisis: The cumulative seismic energy release plots reveal a significant difference between the local VT seismicity and the seismicity at distance greater than 25 km from the volcano.  as the lowest magnitude at which 100% of the events are detected [Rydelek and Sacks, 1989;Taylor et al., 1990;Wiemer and Wyss, 2000]. The b-value mostly depends on the properties of the seismic medium and on the nature of the stress source causing the earthquakes. Areas where seismicity is caused by extensive regional stress show b-values near 1.0 [Frohlich and Davis, 1993]. Stress concentration leading to clustering of seismicity and fracture over distances of a few kilometers seems to have a strong effect on the value of b [Ogata and Katsura, 1993;Wiemer and McNutt, 1997;Wiemer and Wyss, 1997]; particularly, stress concentration produced by increased thermal gradients, or rock fracturing caused by hydraulic stresses [Grasso and Sornette, 1998] cause a significant increase in b [Warren and Latham, 1970]. Both of the former processes are common in volcanic regions, and b-values larger than 1 are commonly associated to some type of internal volcanic activity [De la Cruz-Reyna et al., 2008]. The analysis of the b-values calculated over short time windows in a volcanic region permits calculating the changing "mean recurrence time" (MRT), a parameter measuring the probability of a VT earthquake exceeding a given (relatively high) magnitude. An application of this methodology is illustrated in Figure 3. The methodology was developed as a forecasting tool for the expected largest-magnitude earthquake in each magma injection process during the El Hierro seismo-volcanic crisis [García et al., 2016]. Combining the MRT with an associated alert level procedure (VAL) provides a robust method to objectively communicate the seismo-volcanic hazard [García et al., 2014b].

1 ANALYSIS OF CONTINUOUS SEISMIC SIGNALS
The seismic signals recorded by a seismometer result from the superposition of multiple signals from different sources, including those generated by the volcanic activity. A main goal is to discriminate the signals generated by the volcanic activity from the rest. This is done by means of decomposition techniques such as Singular Value Decomposition [Carniel et al., 2006a], Independent Component Analysis and Non-negative Matrix Factorization [see Carniel, 2014 for a review]. This task can be particularly difficult at the onset of volcanic activity that is, when the signals of volcanic origin are relatively weak. In such a situation, relevant information may be disregarded as "noise". The analysis of noise acquires then a particular importance [Peterson, 1993]. Recent developments in noise analysis based on continuous assessment of noise levels for quality control purposes [Sleeman and Vila, 2007], variations of signal-to-noise ratios [Bormann, 2002;McNamara et al., 2009], or specific applications to geology [Ladina, 2012], resulted in useful tools for the recognition of subtle changes in the internal state of a volcano. Changes systematically altering the amplitude and/or predominant frequencies of the background seismic noise can be detected using these methods of analysis in a given volcanic region. It is particularly important to mention the Base Level Noise Seismic Signal (BLNSS) [Vila et al., , 2008, which proved to be very useful in monitoring the 2003-2005 unrest episode of the Teide volcano. Through this method it was possible to report changes in the internal process of the Teide volcano. As an example, Figure 4 shows the base level noise spectra computed for different days. A significant increase of noise energy in the 1-10 Hz frequency band from 2004-12-01 through 2004-12-03 preceded by nine days a new fissure with fumarole emissions that appeared in La Orotava valley in the northeastern flank of Teide volcano, which returned to the normal levels a few hours after the fissure aperture [see García et al., 2006 for details].
To assess the capabilities of the method, we performed BLNSS analyses on other active volcanoes such as Soufriere Hills in Montserrat, West Indies, and in Llaima and Villarrica in Chile . Generally speaking, the BLNSS proved to be a powerful tool that also provides relevant information on the status of volcanic unrest.
It is important to emphasize that the methods based on the characterization of background noise detect changes in the evolution of a volcanic seismic signal but they do not provide any a-priori information about the sources, nor about the structures traveled through by transitory seismic waves. However, new developments adapted from neurosciences like the Self-Organizing Maps, SOM [Kohonen, 1982] have been recently applied in geophysics along with decomposition techniques such as the Singular Spectrum Analysis to identify specific evolutionary patterns [Carniel et al., 2006bRoden et al., 2015]. Similar techniques were applied in volcanology to identify states of the volcanic activity from continuous tremor data [Carniel et al., 2006a;Langer et al., 2009;Carniel et al., 2013]. The geostatistical approach at different time scales is also promising [Jaquet and Carniel, 2003;Jaquet et al., 2006]. of Kalman filtering [Kalman, 1960], particularly when applied to other strain-rate related parameters such as the terrain deformation, as it was the case in the recent volcanic process at the Canary island of El Hierro in each case on a preliminary assessment of the physics controlling the volcanic activity and the type of available data. This sometimes requires developing new methods to deal with particular situations [see e.g. Salvage and Neuberg, 2016]. A review of other methods of data reduction to characterize volcanic activity can be found in Carniel [2014], where particular emphasis is put on approaches based on the time evolution of parameters that consider the volcano as a non-linear dynamical system.
The identification and follow up of the evolution of different volcanic signals may have a high forecasting value. However, monitoring and recognition of the different signals is a difficult task, particularly during unrest periods, when signals of many coexistent processes are recorded. Recently, some methods of automatic recognition and selection of families of signals in real or quasi-real time have been developed. For example, the Hidden Markov Models, that have been applied to several volcanoes such as Merapi [Ohrnberger, 2001;Wassermann and Ohrnberger, 2001], Las Cañadas caldera [Beyreuther et al., 2008], Popocatepetl [Cortes et al., 2009], Stromboli, Etna  and Volcán de Colima, where the results of automatic recognition accuracy, reported by the group in charge of the monitoring activities amount to 82% [Cortes et al., 2009;Arámbula-Mendoza et al., 2018].

A CASE HISTORY: EL HIERRO 2011-2014 VOL-CANIC ACTIVITY AFTER A LONG REPOSE PE-RIOD
The recent volcanic and seismic unrest observed at El Hierro (Canary islands;2011 provided with an opportunity to observe its complex evolution since its onset in July 2011 using some of the methods described in previous sections, and collect relevant forecast information data to assist decision-making. Seismic data, including earthquake magnitudes and locations were obtained from the on-line official catalog of the Spanish Geographical Institute (http://www.ign.es/ign/lay-outIn/sismoFormularioCatalogo.do). To process this information, we first implemented a real-time calculation of the cumulative seismic energy and a stepwise estimate of the Gutenberg-Richter b-value [Gutenberg and Richter, 1944] on intervals showing stability in the rate of occurrence of different earthquake magnitudes [García et al., 2016]. From the beginning of the seismic analysis, the evolution of the seismic energy conformed to the Yokoyama criterion [Yokoyama, 1988] showing significant increments in the number and magnitude of earthquakes when the cumulative energy exceeded 10¹¹J (Figure 7).
Kalman-filtered data were used to track the SSEM evolution. In the initial stages of unrest, the b-value provided information about the process of magma migrations causing the unrest, but difficulties with the stability of the low-cut-off magnitude made it difficult to continuously track changes of the b-value. The Kalman-SSEM based forecasting of major earthquakes provided useful results until the onset of the submarine eruption (2011-10-10), when the strong tremor caused by the eruption masked most of the other signals. Nevertheless, the SSEM method still allowed forecasts using distant stations.
Volcano-tectonic seismicity at El Hierro stopped when the eruption began, resuming some days later, on 2011-10-29, when a new seismic swarm and significant deformations developed in a way similar to the preeruptive unrest. The new unrest episode culminated with the largest earthquake with magnitude 4.6. Figure  8 shows the improved inverse SSEM forecast of this new unrest episode using the Kalman-filtered data. This activity could be explained by a new deep magma injection process. Afterward, a sequence of similar episodes of seismicity and deformation were recorded, with temporal and spatial evolutions that are consistent with the model of magma injections along paths defined by stress controlled conical surfaces [De la Cruz-Reyna and Yokoyama, 2011;García et al., 2014a]. The model allows detection of shallow regions of increased seismic and related landslide hazard. The similarities among the different episodes of seismicity motivated the definition of a new family of magma-migration related seismic swarms [Farrell et al., 2009;García et al., 2016].

VOLCANIC ACTIVITY FORECASTING AND RISK MANAGEMENT
The assessment and management of developing volcanic activity evolve in two parallel lines; one is specifically scientific, while the other is focused on the management of the risk posed by the activity. The latter must provide decision factors in near real time and include criteria and communication procedures to advise authorities and warn the population sufficiently in advance of a critical situation. Timely needs are determined by the type of necessary preventive actions (e.g. hours in the case of closing of a road in a zone of probable landslides, or days in cases requiring to evacuate). Marrero et al. [2010] presented a new tool for simulating and optimizing large-scale evacuation processes: The Variable Scale Evacuation Model (VSEM), which allows simulating an evacuation considering different strategies for diverse impact scenarios. The decision to evacuate is one of the critical issues in managing a crisis. In order to make such an important decision, it is essential to estimate the potential loss of lives for each of the expected impact scenarios. To improve the emergency response planning and facilitate the decisionmaking process we also developed a methodology to estimate the number of potential fatalities, which has been applied to the Central Volcanic Complex of Tenerife and to El Chichón volcano (México) [Marrero et al., 2012.
Employing several of the above mentioned forecasting tools, and according to a time window defined by the expected impact scenarios, it was possible to design an automatic "watchdog" system capable to detect actual data trend changes, and reject spurious signals re-FIGURE 7. Evolution of the cumulative seismic energy from the onset of unrest at El Hierro island, to the onset of the submarine eruption. The rate of cumulative seismic energy release (blue thick line) clearly accelerates after reaching the Yokoyama threshold level of 10¹¹ J. The purple line represents the inverse of the smoothed derivative of the cumulative seismic energy (inverse energy release rate), which clearly shows the increased rate of energy release just above (below in the inverse curve) that energy threshold. lated to instrumental failure. Such "watchdog" system operates in steps according to the following routine: 1. Automatic and/or visual daily examination of the available seismic data 2. Computation of earthquake magnitudes. 3. VT hypocentral/epicentral locations or at least estimation of distance to station. 4. Cumulative seismic energy calculation, and assessment of the seismic energy release rate. 5. 24-hour spectral analysis of the continuous signal. 6. RSAM calculation and/or derived methods: SSAM, SSEM, BLNSS. 7. MRT calculation using VT located at less than 10-20 km from the volcano. 8. Identification of VT families. In this way unrelated signals are rejected, and only changes related to volcanic activity may be translated into a specific Volcanic Alert System (VAS) that is easy to understand and use by scientists, technicians, and decisionmakers [García et al., 2014b]. Although this procedure has been tested, its effectiveness is highly dependent on the cultural, political and social background of the area in which a crisis develops [Marrero et al., 2015].

DISCUSSION
There are indeed many more methods for volcano monitoring and eruption forecasting with increasing degrees of sophistication. However, a common problem, particularly in developing countries, or in regions with volcanoes with low rates of high-magnitude eruptions, is the difficulty and sometimes impossibility to obtain funding for extensive and permanent monitoring, particularly when other basic social and economic issues must be addressed. The methods discussed here may render useful results even with data obtained with low-cost instruments, like those designed for educational purposes. For example, the Arietta (www.acmesystems.it), Raspberry (www.raspberrypi.org) or Arduino (www.arduino.cc) embedded micro-computers allow assembling seismic monitoring stations at very low cost. Affordable 4.5 Hz geophones with a pre-amplifier can generate valuable data extending their signal response to 5 seconds [Ortiz et al., 2001;Havskov and Alguacil, 2006]. Data may be transmitted to a processing center by GSM, Ethernet or WIFI in quasi-real time or by specific timeperiods depending on the available links, and following 11 FORECASTING VOLCANIC ACTIVITY the procedure described in the previous section. These methods have been used at El Hierro and El Chichón volcanoes with good results

CONCLUSIONS
Awareness of manifestations precursory to eruptions has existed for a long time. Disasters happen when precursors are not recognized, or misunderstood, or ignored, and a wrong decision follows. It is thus necessary to increase the reliability of the precursor's recognition and interpretation, and improve the communication of the hazards and the confidence of the decision-makers. This task is not simple, as the number of parameters involved in a volcanic process exceeds the number of measurable quantities that can be accurately measured at the surface. It is now recognized that even simple dynamical systems can show a complex behavior when they become non-linear [Lorenz, 1963;Mandelbrot, 1977;Turcotte, 1987;Newman, 2011] and this triggered the introduction of novel methods to characterize [Carniel, 2014] and model [Pshenichny et al., 2003[Pshenichny et al., , 2009] the time series representing the dynamics of volcanic systems. For instance, the determination of powerlaw distributions for different parameters, in particular frequency-magnitude power-law statistics, may reveal some general scale-invariance and self-organization characteristics of the magmatic systems. Such information can help to recognize a process that will culminate as a catastrophic failure. Considering the technical developments of sensors, transducers, computers and communication devices, and the persistently increasing capabilities of software to process that type of information, it is presently feasible to monitor volcanoes and maintain a sound forecasting capability at an affordable cost. What is needed is to promote the political will to invest in such systems.