Volcanic eruption time forecasting using a stochastic enhancement of the Failure Forecast Method

In this study, we use a doubly stochastic model to develop a short-term eruption forecasting method based on precursory signals. The method enhances the Failure Forecast Method (FFM) equation, which represents the potential cascading of signals leading to failure. The reliability of such forecasts is affected by uncertainty in data and volcanic system behavior and, sometimes, a classical approach poorly predicts the time of failure. To address this, we introduce stochastic noise into the original ordinary differential equation, converting it into a stochastic differential equation, and systematically characterize the uncertainty. Embedding noise in the model can enable us to have greater forecasting skill by focusing on averages and moments. In our model, the prediction is thus perturbed inside a range that can be tuned, producing probabilistic forecasts. Furthermore, our doubly stochastic formulation is particularly powerful in that it provides a complete posterior probability distribution, allowing users to determine a worst-case scenario with a specified level of confidence. We verify the new method on simple historical datasets of precursory signals already studied with the classical FFM. The results show the increased forecasting skill of our doubly stochastic formulation. We then present a preliminary application of the method to more recent and complex


Theoretical introduction
We have introduced a new method for performing short-term eruption forecasting, when eruption onset is related to a rupture of materials.
• The method enhances the well known FFM equation. We allow random excursions from the classical solutions. This provides probabilistic forecasts instead of deterministic predictions, giving the user critical insight into a range of failure or eruption dates, and allowing retrospective evaluation and improvement of forecasting methodology.
• Our doubly stochastic formulation can consider the "worst case scenario" with a probability of occurrence of at least 5%. This was not possible in the classical formulation.
• We tested the method on historical datasets of precursory signals. The data show the increased forecasting skill of the doubly stochastic formulation, expressed as the likelihood in the day of the actual eruption.
• We also described an assessment of failure time on present-day unrest signals.
The new formulation enables the estimation on longer time windows of data, locally including the effects of variable dynamics.
This approach is the subject of ongoing and future work, with the purpose to further testing forecasting robustness, for example, exploring the sensitivity on a linear evolution of α with time, or a more general structure of the noise.  Cornelius & Voight, 1995) Paper Number V23G -0155 Abstract ID: 413317

Preliminary application to present-day unrest 5. Conclusions
The Failure Forecasts method (FFM) is a well-established tool in the interpretation of monitoring data as possible precursors, providing quantitative predictions of the eruption onset t e , commonly represented by inverse rate plots (Voight, 1988;1989 Fig. 1).
The model (Fig. 2) represents the potential cascade of precursory signals leading to a large-scale rupture of materials culminating at the failure time t f . The equation was originally developed in landslide forecasting (Fukuzuno, 1985).
Laboratory experiments and theoretical models demonstrated the use of the FFM under constant stress and temperature (Hyman et al., 2018). Without this assumption, more fundamental relations between rock fracture and deformation imply time dependent changes in the power law properties (Robertson&Kilburn, 2016;Kilburn, 2018).
Our aim is to produce probability forecasts with the FFM, instead of deterministic predictions. We embed noise in the model and enable sudden changes in the power law properties (Fig. 3).
We assume that the FFM equation is not exactly satisfied, but there is a transient difference produced by additive noise, which however decreases exponentially in time. This is called a Hull-White model in finance (Hull&White, 1990).
FFM is known to be affected by sources of uncertainty, like: • the occurrence of multiple phases of acceleration in the signals (Bouè et al., 2015) • varying stress on the system, or temperature (Kilburn, 2012) • the superposition of signals originating from different causes (Salvage&Neuberg, 2016) • heterogeneity in the breaking material, producing changes in the signals (Vasseur et al., 2015).
In addition, the statistical fitting of model parameters can be poorly constrained (Bell et al., 2011).
In this study, we enhance the classical FFM by systematically characterizing the uncertainty. We incorporate stochastic noise in the equations, and a mean-reversion property to constrain the noise.
We describe three different methods for estimating t f , the ODE-based Method 1, the new SDE-based Method 2, and their combined doubly stochastic formulation Method 3.
Doubly stochastic models describe the effect of epistemic uncertainty in the formulation of aleatory processes (Cox&Isham, 1980;Bevilacqua2016).

CLASSICAL FFM FORMULATION 1.2 ENHANCED FFM FORMULATION
Log-rate vs Log-acceleration Technique (LLT), and Hindsight Technique (HT) are classical estimators of α (Cornelius & Voight, 1995). LLT is less accurate than HT and needs a second order derivative of data. HT needs to know t e and can only be used retrospectively.
We apply three different forecast methods on the datasets in Voight, (1988), and we test t f as an estimator of t e .
• Method 1 (Fig. 4) characterizes the epistemic uncertainty related to the parameter fitting in the classical FFM.
Blue lines assume α as from LLT, red as from HT. The bold line is g tf , the probability/day scale bar is related to it. In Fig.6 bold dashed lines are its 5 th and 95 th percentile values.
Thin dashed lines bound the 90% confidence interval of the ODE paths of 1/X, and a thin continuous line is the mean path.
Black points are inverse rate data.
LLT is not well constrained in this case.
• Method 2 (Fig. 5) allows excursions from the classical FFM solutions, modeled as aleatoric uncertainty sources. The numerical solution of the SDE is performed as 5,000 discretized random paths by the Euler-Maruyama method.

Fig 7.
Column plots of the likelihood g tf (t e ), i.e. the probability of failure time t f in the correct day t e . In plot (a) the black bars assume Method 1, the colored bars Method 2. Plot (b) assumes Method 3, and the full bars are the mean values, and the shaded bars are the 95 th percentile values of the likelihood. If α is based on LLT Method 1 provides low likelihoods, below 1% in some case (Fig. 7). The 95 th percentile values of Method 3 clearly outperform other methods. We remark that the HT method cannot be used in forward forecasting, but only retrospectively.
Estimators based on the whole sequence of signals are not forecasts. We display forecasts of t e obtained from parameters inferred from the data collected in a limited time window T=[t 2 ,t 1 ]. For simplicity, α is still based on the entire sequence of processed data in Voight, (1988).
Method 1 (Fig. 8) and mean pdf in Method 3 simulate 20,000 samples, Method 2 (Fig. 9) is based on 5,000 random walks. The percentile values in Method 3 (Fig. 10) are based on 150,000 samples We compare two time windows with extremes reported in figure. They include different subsequences of data. Forecasts can be significantly uncertain, because based on fewer data. In Methods 1 and 3, there is a chance that the solution path never hits the real axis.
When the forecast is poorly constrained, Method 2 typically reduces the uncertainty affecting t f , compared to Method 1. It tends to give a correct forecast only when the eruption is close. The doubly stochastic formulation of Method 3 appears to have an impact.
Method 3 mean forecasts provide consistent likelihoods with Method 1 (Fig. 11). The 95 th percentile values are significantly higher than other forecasts, from 5% to 10% in the first and second time windows, and above 15% in the third.  Fig 8,9,10. Forecasts of t f based on Method 1 (Fig. 8), Method 2 ( Fig. 9), and Method 3 (Fig. 10).

Probability forecasts on historical datasets
In (a,b) and (d,e) two examples are tested on different time windows T.
The bold line is g tf , the probability/day scale bar is related to it. In Fig.6 bold dashed lines are its 5 th and 95 th percentile values.
Thin dashed lines bound the 90% confidence interval of the ODE paths of 1/X, and a thin line is the mean path.
Points are inverse rate data, those in T are colored.  The green line is mean value of g tf , the probability/day scale bar is related to it. Dashed lines mark its 5 th and 95 th percentiles.
Thin blue dashed lines bound the 90% confidence interval of the ODE paths of 1/X, and a thin line is the mean path. Grey dotted liner display 50 sample paths of the SDE solution after 10 Sep 18.

PROBABILISTIC FAILURE TIME
• Method 3 (Fig. 6) combines the other methods. The mean solution is affected first by the parameter uncertainty and then perturbed by the stochastic noise. g tf is reported as 5 th percentile, mean, and 95 th percentile curves. Mean is based on 10,000 samples. Percentile values on 60,000 samples.
When α ≈ 2, all methods provide good estimators of t e . When α < 1.6, Method 1 can overestimate t e . Noise reduces this issue and can push 1/X to zero earlier. Method 3 performs well, and its formulation allows users to determine a "worst case scenario" with a specified level of confidence .
We assume γ=1/15 days. This is a choice based on the empirical observation that the total length of temporal sequence is at the scale of 45 days, and the duration of well-aligned observations is at the scale of 15 days. Doubled or halved γ produced minor effects on the results.
Campi Flegrei caldera is characterized by prolonged unrest. Since 1950, it has undergone four episodes of caldera-wide uplift and seismicity, which raised the central region by 4.5 m (Troise et al., 2019).
After about 20 years of subsidence, following the uplift peak reached in 1984, the caldera started a new, low rate uplift episode, accompanied by low magnitude increasing seismicity (Fig. 12) and marked geochemical changes in fumaroles. Although the interpretation of the current unrest is a matter of debate, Campi Flegrei may be evolving towards conditions more favorable to eruption (Kilburn et al., 2017) We focus on the earthquake count after 2007, available on www.ov.ingv.it. Seismicity is characterized by repeated swarms, superimposed on a background rate (Chiodini et al., 2017). Swarms are believed to be related to fluid-transfer episodes (D'Auria et al., 2011). Further analysis may test the effect of separating the swarms from background rate.
In Figure 13 we test our Method 3 on this seismic dataset. We remark that the time window is much longer than in the classical applications, and spans over 50 years.
The interpretation of t f as the onset of a volcanic eruption is purely speculative. However, it is the time when accelerating signals as observed in the last 10 years will diverge to infinity. A few pauses in the seismic rate are filtered out, and are reported in Fig.13a,c. Locally decreasing rates prevent the LLT to calculate α, and we uniformly sampled α in [1.4, 2.0]. Here γ = 1/120 days.
Estimates of t f are in [2019,2032] with 90% chance, at the scale of 0.02% probability per day. These are robust against the choice of the bin on which the rate is calculated (Fig. 13a,b). However, if we consider only the last 5 years, they become sensitive to that choice.
LLT is not well constrained in this case. LLT is not well constrained in this case.