Short-term volcano-tectonic earthquake forecasts based on a MRT algorithm: the El Hierro seismo-volcanic crisis experience

Under certain conditions volcano-tectonic (VT) earthquakes may pose significant hazards to people living in or near active volcanic regions, especially on volcanic islands; however, hazard arising from VT activity caused by localized volcanic sources is rarely addressed in the literature. The evolution of VT earthquakes resulting from a magmatic intrusion shows some orderly behavior that may allow forecasting the occurrence and magnitude of major events. Thus govern5 mental decision-makers can be supplied of warnings of the increased probability of larger-magnitude earthquakes in the short term time-scale. We present here a methodology for forecasting the occurrence of large-magnitude VT events during volcanic crises; it is based on a Mean Recurrence Time (MRT) algorithm that translates the Gutenberg-Richter distribution parameter fluctuations into timewindows of increased probability of a major VT earthquake. The MRT forecasting algorithm was 10 developed after observing a repetitive pattern in the seismic swarm episodes occurring between July and November 2011 at El Hierro (Canary Islands). From then on, this methodology has been applied to the consecutive seismic crises registered at El Hierro, achieving a high success rate in the real-time forecasting, within 10 day time-windows, of volcano-tectonic earthquakes


El Hierro volcanic process
The morphology of El Hierro island (27.7 o N; 18.0 o W) has been attributed to a triple volcanic rift with the NE, NW and S axes of these rifts diverging by about 120 o from each other (Carracedo et al, 1999).The island landscape displays abundant dikes and fissures (Gee et al, 2001), as well as major 60 landslides (Masson et al, 2002).The subaerial island displays a high concentration of volcanoes with a spatial distribution that may be explained by the model proposed by Stroncik et al (2009), in which the volcanic activity at El Hierro is controlled by a complex array of magma pockets.Magma injections from the pockets through dikes or sills may generate VT seismicity.
In July 2011, a sudden increase of seismicity was detected at El Hierro Island, prompting several 65 Spanish institutions to deploy recording instrumentation.In addition to the official seismic network operated by the Instituto Geográfico Nacional (IGN), in September 2011, the IGEO-CSIC Institute (see authors' affiliations) installed, for research purposes, a seismic network consisting of 3 permanent, real-time, vertical, short-period stations, and 5 similar, temporary, manual-access stations.
Seismicity, deformations and diffuse gas emissions increased persistently until October 2011, when 70 a submarine eruption at 27 o 37.18' N; 17 o 59.58' W was first detected (López et al, 2012;Ibáñez et al, 2012;Prates et al, 2013).The eruptive activity persisted for months, until March 2012.However, neither the seismic activity nor the deformations stopped, and have in fact continued to the time of this submission (García et al, 2014b, a).
The evolution of the seismic activity at El Hierro island for the whole period of unrest (2011-75 2015) is summarized in Figure 1, in terms of the time distribution of earthquake magnitudes, the number of earthquakes exceeding M1.5 per day, the cumulative number of events, and the cumulative seismic energy released, estimated from the IGN public earthquake catalog (http://www.ign.es/ign/layoutIn/sismoFormularioCatalogo.do as downloaded on 2015-08-29), and from records of the above-mentioned IGEO-CSIC seismic research network.The uncertainties in calculating the dura-80 tion and magnitude of very low-amplitude signals recorded by the IGEO-CSIC have been treated according to Del Del Pezzo et al (2003).

Figure 2 near here
Inspection of Figure 1 reveals that the seismicity, deformations and eruption evolution observed during the 2011 unrest episode not only is consistent with the model of Stroncik et al (2009) (García et al, 2014b), but it also reveals some other interesting features of the volcanic unrest process dynamics.First, as expected, a marked rise of the VT event occurrence rate preceded the onset of the 10 October 2011 eruption.However, during the continuing seismic unrest, the event rate and the cumulative number of events have also markedly increased before the occurrence of major VT earthquakes.But it is the cumulative release of seismic energy, shown in the uppermost plot of 90 Figure 1, the parameter that best illustrates the precursory character before major VT events.Additionally, the M-T diagram in the lower plot of Figure 1 shows a peculiar behavior.The minimum magnitude rises as the number of higher-magnitude events increases.This may be caused by com- Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2015-273, 2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.pleteness fluctuations of the catalog, derived from the difficulties in counting very small events when larger earthquakes are concurrent, or it may be an actual condition related to the physical process 95 of seismic energy release.Several other cases in which the completeness is maintained at the lowest measurable value during the duration of the VT swarm suggest that the latter condition may be real.
For example, Benoit and McNutt (1996)  To date, five seismo-volcanic crises or VT earthquake swarms have been identified at El Hierro (García et al, 2014b, a), showing in all cases an increasing trend of the maximum detected magnitudes, as shown in Figure 2. It is significant that the acceleration of the seismic energy release rate begins in all cases when the cumulative seismic energy exceeds the 10 11 J energy threshold criterion 105 of Yokoyama (1988).Hours to days after the onset of this acceleration stage, the swarm sequences culminated with earthquakes exceeding magnitude 4 (Figure 3).Currently, VT earthquake swarms are classified according to three family types proposed by Benoit and McNutt (1996), and later updated by Farrell et al (2009) and Brumbaugh et al ( 2014): Type A-mainshock and aftershocks; type B-foreshocks, mainshock and aftershocks; and type C-swarms without a mainshock.Further exam-110 ination of the lower part of Figure 2 reveals that the evolution of the daily rate of seismicity does not correspond to any of these established earthquake families.We thus introduce a fourth family, type D, characterized by an increasing number of growing magnitude foreshock events culminating with a mainshock and followed by a rapid decay of seismicity lasting just a few hours.The differences between these models are illustrated in Figure 4. Therefore, a type D sequence does not follows 115 an ETAS model (Ogata, 1992), nor follows the Omori law (Utsu and Ogata, 1995), but it reveals a causal process that allows forecasting of major VT events.

Forecasting major VT earthquakes at El Hierro
In a time-evolving seismicity, the Gutenberg-Richter Law (GRL, Gutenberg and Richter, 1944) scales the activity, with respect to its magnitude, as: Equation 1 here where N (M )/∆T is the number of earthquakes with magnitude M or larger detected within a 125 given region, in a time interval ∆T .In a process that remains stationary over ∆T , the Gutenberg-Richter Parameters (GRP) a and b remain constant.Here, we use the time-dependent GRP to estimate the Mean Recurrence Time (MRT) of volcanic earthquakes having or exceeding given magnitudes, and thus the likelihood of major events occurring in the shorter time-scale (hours to days).The MRT Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2015-273, 2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.
(t T ) between events with magnitude equal or greater than M may thus be estimated from Equation 1130 as: Equation 2 here Increasing likelihood of major events may thus be estimated from the increasing rates of occurrence of smaller earthquakes.In practice this may prove to be a difficult task since the available seismic catalogs frequently present completeness uncertainties derived from technical factors such 135 as the stability of the seismic networks and links, weather conditions, background noise, event counting and magnitude assignment routines, among others.Problems of catalog completeness in general mean an under-counting of the lowest magnitudes (Bell and Kilburn, 2012) that define a low-end Magnitude Cutoff (MC) above which the GRL holds.A second major problem is that the evolution of VT seismicity in the El Hierro crisis is highly variable, and the estimates for the GRL parameters 140 may vary in time and space, causing local distributions of shocks quite different from those of the full catalog.This problem has been amply discussed in the literature about seismicity from different origins (e.g.Wiemer et al, 1998;Lombardi, 2003;Marzocchi and Sandri, 2003;Helmstetter et al, 2007;Bengoubou-Valerius and Gibert, 2013;Alamilla et al, 2015;Márquez-Ramírez et al, 2015).
We have addressed those problems considering a moving and variable time-window ∆T in equa-145 tions 1 and 2. In most cases, the time-window starts running when the cumulative energy released by a VT swarm reaches the 10 11 J energy threshold of Yokoyama (1988).Experience has shown that at El Hierro such swarms precede a stage of accelerated seismicity with increasing magnitudes that culminate in major (M>4) VT events.
The MRT algorithm starts a forecast when a swarm of at least 200 earthquakes with M≥1.5 are 150 detected in a time span of 5 days.Then, a completeness magnitude MC is calculated using the Maximum Curvature Method (MCM, Wiemer and Wyss, 2000).The GRL parameters are estimated for each time-window using a Least-squares Regression Estimator (LRE) and/or a Maximum Likelihood Estimator (MLE) (Peishan et al, 2003;Bengoubou-Valerius and Gibert, 2013).In general, the LRE performs better for seismic swarms, as it allows the contribution of the dominant magnitudes of the 155 swarm to be emphasized underscoring only the linear part of the distribution.In successive swarms, some patterns appear as the time to a mainshock approaches: MC may increase from the initial 1.5 to as much as 3 or more, and the time window shortens.Algorithm 1 shows a schematic description (pseudocode) of the algorithm.The algorithm has been implemented in a batch processing mode that may be activated every 24, 12 or 6 hours, depending on the level of activity.

160
Algorithm 1 near here Table 1 shows the main features of the swarms detected at El Hierro, and Figure 5 shows the evolution of the MRT for the initial four cycles of El Hierro activity.A warning window for the likely occurrence of an earthquake exceeding M4 is issued when the calculated MRT for that magnitude drops below 10 days, and a statement saying that there is an increased likelihood of an M4 (or higher)   An 8-month lull followed those four cycles, and as frequently happens, a belief that the seismo-170 volcanic crisis had ended spread throughout the island's population.Additional factors, such as a particularly strong winter storm and the Christmas holidays, caused failures of the IGN official seismic network that translated into a rather incomplete seismic catalog.When a new seismic swarm began on 2013-12-20, the MRT could only be calculated using predominantly data from the IGEO-CSIC seismic research network.Figure 6 shows the MRT evolution using magnitudes from both 175 catalogs.On 2013-12-24 a new warning was issued communicating to decision-makers the increased likelihood of a M5 or larger earthquake occurring in the following days.On 2013-12-27, at 17:46:02 UTC, the largest earthquake of the entire episode, with magnitude 5.4, caused numerous landslides across the island, fortunately without casualties.It is interesting to note that no M4 earthquakes occurred in this stage.Earthquake forecasting is a problem that has been addressed in the past from various different perspectives.Significant advances were obtained after the introduction of complexity theory (Lorenz, 1963;Mandelbrot, 1977) and the development of predictive algorithms (Gabrielov et al, 1986; 185 Molchan, 1990;Keilis-Borok, 1996;Turcotte, 1997).Currently, the need to implement operational short-term earthquake forecasting based on predictive and forecasting algorithms has been stated by Jordan et al (2011).The MRT is one of such algorithms, based on a property that is typical of complex systems: seismic events cluster according to a simple scaling law, the Gutenberg-Richter relationship (Equation 1).This empirical relationship contains information about the occurrence of 190 large magnitudes embedded in the number of small-magnitude events (Helmstetter et al, 2007), and the time variations of the distribution of the small-magnitude events provides information about the time-scale of large-magnitude occurrences.In the case of VT earthquakes, two additional factors help to simplify the forecasting process: one is the limited geographic extent of the volcanic source of seismicity, confining the spatial extension of the forecast.The other factor is the short time-scale, 195 measured in days, at which orderly fluctuations of the GRPs evolve reflecting stress accumulations induced by the magma injections, and allowing rapid recognition of significant changes of the MRT as defined in Equation 2, and therefore viable short-term forecasts (García et al, 2014a).In contrast, a similar method for tectonic earthquakes would require recognizing orderly evolution patterns of the 6 Nat.Hazards Earth Syst.Sci. Discuss., doi:10.5194/nhess-2015-273, 2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.
GRPs over times scales of years to decades, making short-term forecasting non-viable (e.g.Jordan We envision the El Hierro process as the result of magma injections causing stress distributions capable of triggering seismic swarms.Such swarms culminate with mainshocks that release stress concentrations, and are thus not followed by aftershock sequences.Further magma migrations cause new swarms and mainshocks, in repeating, but non-periodic cycles of activity.For example, the 205 largest earthquakes of the March and December 2013 swarms (see Table 1), were located at less than 20 km West from the westernmost shore of the island, suggesting a nearby localized stress source possibly resulting from subsurface magma injections under that region (Figure 7).From the initial stage of unrest in July 2011, five earthquake swarm cycles have been identified, but only the first of these culminated in an eruption.The subsequent cycles show two main features: first, the 210 swarms are separated by increasingly longer lull intervals, and second, they tend to culminate in mainshocks of increasing magnitudes (see Figure 2).

Figure 7 near here
The forecasting algorithm we are proposing here is based on the above-mentioned envisioning of the process that is causing the seismic swarms, and begins when the background seismic energy 215 exceeds the 10 11 J energy threshold criterion of Yokoyama (1988).The acceleration of the earthquake rate (Figure 4-D) is interpreted in terms of piecewise estimations of the Gutenberg Richter parameters, The variations of those parameters are, in turn, translated as time-windows of increased likelihood of a large-magnitude earthquake.The MRT algorithm may be, and has been, used in real time.In the absence of swarms, the GRPs are calculated with a long database of the seismic cata-220 log, and updated on a daily basis (Algorithm 1), usually yielding fluctuating, but always very long values of the MRT (that may exceed 106 days).At the onset of a seismic swarm, the MRT drops rapidly to values of the order of 103 days.When this happens the GRPs are calculated from shorter time-windows of variable duration, provided that they contain a large enough number of events (at least 200) to yield stable values of the GRPs.When the MRT for events of magnitude 4 or higher is less than 10 days, a warning is issued to the decision makers, expressed as a 10-day window of increased likelihood for larger, potentially damaging earthquakes.Up to the time of this submission, all the recent earthquakes of M≥4 at El Hierro have been forecasted several days in advance.Application of the MRT to other volcanoes may be relatively simple to implement, provided that the observed seismicity corresponds to swarms of type D (Figure 4): increasing number of VT events 230 with increasing magnitudes followed by a mainshock, and by a rapid decay of the seismicity afterwards.The algorithm may be readily adapted to any real-time monitoring system.In addition, the concept of time-windows for the increased likelihood of a large-magnitude earthquake should be easily understood and accepted by governmental decision-makers, making the management of risk more effective.
report an increase of the minimum magnitude when the VT activity approaches a mean magnitude 3. Likewise, De la Cruz-Reyna et al (2008) and Hurst et al (2014) report a similar behavior for the VT seismicity of Popocatépetl and Tongariro volcanoes 100 respectively.

Figure
Figure 2 near here Figure 3 near here Figure 4 near here 120 Earth Syst.Sci.Discuss., doi:10.5194/nhess-2015-273,2016   Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.earthquake in the following 10 days is delivered to the decision-makers.The warning remains active until the MRT returns to a value exceeding 10 days.

Figure
Figure 5 near here 180

Figure
Figure 6 near here

Figure 2 .
Figure 2. Time evolution of the daily count of VT events (red dots) and maximun magnitudes (gray bars) during the earthquake swarms of the El Hierro seismo-volcanic process.

Figure 3 .
Figure 3. Seismo-volcanic crises identified at El Hierro Island (2011-2014).The upper plots show the evolution of seismic energy for each of the cycles of unrest.The lower plots show the daily event count.The plots illustrate how the seismic activity drops rapidly after the largest magnitude earthquakes.Particularly significant is the absence of aftershocks after the largest (M5.4) earthquake of 2013-12-27.Data are from IGN public earthquake catalog and the IGEO-CSIC seismic research network.

Figure 4 .
Figure 4. Seismic swarms are usually classified in three family types (A, B, C, Farrell et al, 2009; Brumbaugh et al, 2014).From our results at El Hierro in 2011-2015, we propose a fourth family type D, characterized by an increasing number of VT earthquakes and their magnitudes with time, and by a quick decay of the seismicity, in a matter of hours after a mainshock usually exceeding M4.

Figure 6 .
Figure 6.MRT evolution in December 2013 calculated from the IGN official catalog (blue) and from the IGEO-CSIC catalog (red).The magnitude of the earthquake on 2013-12-27 17:46:02 was estimated by IGN as 5.1, and as 5.4 by the IGEO-CSIC network and other international catalogs.

Figure 7 .
Figure 7. Swarms and mainshock epicentres of the March and December 2013 seismicity cycles.The M5.4 mainshock of 2013-12-27, at 17:46:02 UTC is located in the same region of the March 2013 cycle.