Forecasting Volcanic Eruptions: Beyond the Failure Forecast Method

Volcano-tectonic seismicity and ground movement are the most reliable precursors to eruptions after extended intervals of repose, as well as to flank eruptions from frequently active volcanoes. Their behavior is consistent with elastic-brittle failure of the crust before a new pathway is opened to allow magma ascent. A modified physical model shows that precursory time series are governed by a parent relation between faulting and elastic deformation in extension, subject to independent constraints on the rate of crustal loading with time. The results yield deterministic criteria that can be incorporated into existing operational procedures for evaluating the probability of crustal failure and, hence, levels of alert during an emergency. They also suggest that the popular failure forecast method for using precursory time series to forecast eruptions is a particular form of the parent elastic-brittle model when rates of stress supply are constant, and that magma transport and crustal fracturing during unrest tend toward conditions for minimizing rates of energy loss.


INTRODUCTION
Most of the world's active volcanoes are not regularly monitored (Tilling, 1995;Sparks et al., 2012), especially those that have been in repose for centuries or more. When such volcanoes reawaken, short-term forecasts of eruption are normally based on data from rudimentary monitoring networks installed after the start of unrest. Local microseismicity, or volcano-tectonic (VT) events, and ground deformation are the pre-eruptive phenomena most frequently detected (Zobin, 2003;McNutt, 2005;Dzurisin, 2007;White and McCausland, 2016) and forecasting strategies use accelerations in both signals to estimate when their rate will reach a critical pre-eruptive value. A popular approach follows the failure forecast method (FFM) developed by Voight (1988) and Cornelius and Voight (1994,1996. However, the best-fit trends for individual time series are not unique and statistical fits can yield ambiguous results (Bell et al., , 2013Boué et al., 2015;Vasseur et al., 2015). Here we argue that the precursory time series are governed by a deterministic parent relation between seismicity and deformation, subject to independent constraints on how the crust is loaded with time. The FFM then emerges as a particular form of the parent relation when rates of stress supply are constant. Because the parent relation is generic, it can be applied in the absence of information about previous unrest and offers the prospect of enhancing the reliability of forecasts by integrating deterministic estimates of eruption time with existing probabilistic evaluations. The physical basis for the parent relation is reviewed and updated, before it is used to propose new operational procedures for emergency forecasts of eruptions. FIGURE 1 | Resistance to magma ascent is governed (A) by magma rheology at open volcanoes that maintain a connection between their feeding body and the surface, but by the strength of the crust at closed volcanoes, when a new pathway (arrows) must be formed before (B) a flank eruption (even when an open conduit is also available) and (C) eruptions after an extended interval of repose.

SEISMIC AND DEFORMATION PRECURSORS TO ERUPTION
Precursors to eruption describe how magma is able to overcome resistance to its ascent through the crust. When magma can utilize an existing open conduit, resistance is provided by the magma's rheology and friction along the conduit's walls (Figure 1). When an existing conduit is not used, resistance is dominated by the need to break open a new pathway for magma ascent; this so-called closed condition applies to long-dormant volcanoes, which have sealed previous conduits, and to volcanoes creating new pathways in addition to an existing conduit, such as flank eruptions through volcanic edifices that already have an open summit vent (Figure 1).
A crustal control on opening new magmatic pathways favors similar patterns of precursory behavior, independent of a volcano's tectonic setting, magma composition and style of eruption; similar patterns support the view that short-term forecasts of eruptions are feasible and can be quantified by rock mechanics; and the control of rock mechanics suggests that repeatable patterns are likely to be scale-independent. Voight (1988Voight ( , 1989 was among the first to promote the importance of scale-independence by recognizing the similarities between precursors to eruptions and rock deformation in the laboratory. He argued that restricted ranges of precursory behavior are driven by a positive feedback between the rate (d /dt) and acceleration (d 2 /dt 2 ) of a precursory signal with time: where A is a constant and α, which usually lies between 1 and 2, describes the strength of the feedback. As the rate increases, the acceleration becomes larger and increases the rate more quickly than before. The deterministic forecasting potential of the relation depends on estimating the time at which the rate tends to infinity. This singularity is interpreted to represent when a major change occurs in the structure of the deforming rock that ultimately favors an eruption, such as bulk failure and the formation of a pathway along which magma can reach the surface (Voight, 1988;Cornelius and Scott, 1993;Cornelius and Voight, 1994, 1996Kilburn and Voight, 1998;De la Cruz-Reyna, 2001;Chastin and Main, 2003;Collombet et al., 2003;Kilburn, 2003;Bell and Kilburn, 2011). Constrained by the limited test data available at the time, Voight (1988) proposed that individual precursory sequences are characterized by a single value of α, that α can take values from 1 to 2, and that the associated range of accelerations can be expected for any precursory signal related to ground deformation and fracturing. Subsequent studies, however, suggest that α can evolve from 1 to 2 during an individual precursory sequence (Figure 2; McGuire and Kilburn, 1997;Kilburn and Voight, 1998;Kilburn, 2003Kilburn, , 2012, and that different signals need not accelerate throughout a full sequence. For example, closed volcanoes and rock-physics experiments show that rates of fracturing can accelerate while the deformation rate remains constant (Figure 3; Kilburn, 2012;Robertson and Kilburn, 2016). Under such conditions, Eq. (1) may accommodate the VT event rate, but not the contemporaneous rate of deformation, and so describes only part of the precursory processes that lead to eruption.
(1) has been fundamental in advancing the quality of short-term forecasts: it has revealed the control of rates of change in a precursory signal on determining the approach to eruption (rather than a threshold value of the precursor); it has showed that precursory mechanisms are expected to be scale-independent; and it has demonstrated that numerous, independent empirical relations for rock failure share common dynamic constraints (Main, 1999;Turcotte et al., 2003;Ojala et al., 2004;Davidsen et al., 2007;Schmid and Grasso, 2012). In addition, the preferred values of 1 and 2 for α correspond to exponential and hyperbolic increases with time. Both types of increase are common among self-accelerating processes in biological (Monod, 1949) and physico-chemical (Frank-Kamenetskii, 1939;Bowden and Yoffe, 1952;Gruntfest, 1963;Shaw, 1969) systems and so, in retrospect, ought not to be a surprising feature of accelerating rock failure.

PRECURSORS CONTROLLED BY AN ELASTIC-BRITTLE CRUST
Precursory unrest is usually detected ∼1-10 months before eruption at closed strato-volcanoes (Zobin, 2003;McNutt, 2005;Dzurisin, 2007;Bell and Kilburn, 2011;White and McCausland, 2016), but may continue for decades at closed large calderas (Robertson and Kilburn, 2016). At strato-volcanoes, it commonly involves as many as ∼10 4 VT events, which occur without systematic changes in hypocentral locations, and ground displacements of ∼0.01-1 m over distances of ∼1-10 km (Kilburn, 2003(Kilburn, , 2012Bell and Kilburn, 2011); at large calderas, the total number of detected VT events can be at least an order of magnitude greater (e.g., ∼10 5 events at Rabaul before its eruption in 1994; Robertson and Kilburn, 2016) and ground displacements reach ∼1-10 m over distances of kilometers (Bellucci et al., 2006;Acocella et al., 2015;Di Vito et al., 2016;Robertson and Kilburn, 2016). Most of the detected VT events have magnitudes of 0-2 and are triggered by the movement of faults ∼10 −2 -10 −1 km across, or ∼0.1-10% the size of deforming crust. VT events FIGURE 2 | The andesitic-dacitic volcano Soufriere Hills, Montserrat, erupted on 15 November 1995 after 350 years of repose. VT data from three sources are consistent with initial and final episodes of accelerating VT event rate separated by an interval of constant event rate. (A) 01 January-31 July (time t = days 0-212 since 01 January). The mean daily VT event rate (dashed curve) increased exponentially as 16 exp (t/55.6) [data from Shepherd et al. (2002)]. (B) 13 August-30 September (t = days 225-273). Constant mean VT rate at c. 17 events per day [line; data from Gardner and White (2002)]. (C) 01-14 November (t = days 305-318). The mean daily rate increased hyperbolically (curve) as (0.037-0.0024t) -1 until the eruption [curve; data from Kilburn and Voight (1998) and Kilburn (2003)]. Additional stations were added to the monitoring network between 18 and 31 July [days 200-212;Gardner and White (2002)] Most VT events had magnitudes between 1.3 and 2.7 (Power et al., 1998;Gardner and White, 2002). (D). The three datasets suggest that the mean VT trends evolved with time from quasi-elastic (yellow), through steady inelastic (orange) to accelerating inelastic (magenta); dashed portions of the curve are qualitative interpolations. and ground deformation can thus be viewed as proxies for the inelastic and total deformation of a crust that contains a dispersed population of small faults.
The trends described by Eq. (1) can be incorporated into a general model by recognizing that deformation and VT events are mutually dependent and that their time series are the result of a parent relation between inelastic and total deformation, constrained by specified changes with time in loading the crust. The parent relation describes the potential for rock to fail when supplied with a mean differential stress S sup . For each increment of supplied stress, S sup , a proportion S is used elastically to deform atomic bonds, while the remainder, S loss , is lost inelastically by breaking bonds: The total elastic component determines the mean, or bulk, differential stress S that is established across the deforming volume V of crust and, hence, also the amount of elastic energy being stored. Under idealized conditions, the dominant mode of deformation is initially elastic (without seismicity; S sup = S) and, once faulting has begun, evolves with increasing differential stress from quasi-elastic to inelastic (Figure 4), as the additional energy being supplied is at first mostly stored elastically (quasielastic, S sup ≈ S) and, later, mostly consumed by faulting (inelastic, S sup ≈ S loss ). In the quasi-elastic regime, the mean deformation approximately increases in proportion to S until it reaches its maximum value S F , when the amount of stored energy has reached its maximum capacity. Any additional stress (e.g., from a pressurizing magma body) is consumed in fracturing and fault movement and deformation continues under a constant maintained stress, which defines the inelastic regime. Bulk failure begins within one or more portions of the crust, where the stresses at fracture tips remain large enough to persistently overcome rock resistance and so allow fractures to grow and link together (Griffith, 1921).
Deformation, rather than stress, is measured in the field, so that the quasi-elastic and inelastic regimes have to be identified indirectly from relations between VT seismicity and total ground movement. In the quasi-elastic regime, the number of VT events (a measure of inelastic deformation) is expected to accelerate with the total amount of ground movement (a measure of total deformation) as the proportion of inelastic deformation becomes larger. In the inelastic regime, net deformation continues by faulting, so that the ground movement is expected to increase in proportion to the number of VT events. A complete failure sequence is thus expected to show an accelerating and then  Bell and Kilburn, 2011). In this case, the VT number, N, increases exponentially with time t (in days; broken curve) as 1.17 exp (t/23.3); the mean daily VT event rate increases 0.05 exp (t/23.3). The combination of exponential VT rate and constant deformation rate is consistent with quasi-elastic brittle behavior under a constant rate of stress supply ( Table 1). The eruption occurred almost immediately after the end of the quasi-elastic regime when t/23.3 ≈ 4.
linear increase in the number of fracturing events with ground movement (Figure 4) (Nakata, 2006;Bell and Kilburn, 2011) and El Hierro (Istituto Geográfica Nacional [IGN], 2011;Sagiya, 2011;Kilburn et al., 2017) and for 23 years at Rabaul (Johnson et al., 2010;Robertson and Kilburn, 2016); they include both strato-volcanoes (Kilauea and El Hierro) and large calderas (Rabaul); and the resulting eruptions ranged from basaltic effusive at Kilauea (Swanson et al., 1971) to andesitic phreato-Plinian at Rabaul (Rabaul Petrology Group, 1995). All three VT-deformation trends show the expected evolution from quasielastic to inelastic behavior, independent of the type of volcano, the duration of precursor, magma composition, style of eruption and, as discussed below, how VT event rate and deformation rate vary with time. Such similarity is compelling evidence that elasticbrittle failure of the crust determines the pattern of seismic and deformation precursors to eruptions at closed volcanoes.
Elastic-brittle failure additionally favors the scale-independent behavior implicit by Voight's analysis (Voight, 1988). VT events are produced by fault slip; slip is accommodated by extending the margin of a fault; and fault extension is controlled by cracking within the zone of stress concentration, or process zone, that develops around a fault's tip (Atkinson, 1984;Main and Meredith, 1991;Lockner, 1993;Cowie and Shipton, 1998). Smaller discontinuities (both cracks and slip surfaces) grow in the process zone until it is broken, so extending the parent fault. Similarly, the discontinuities grow by cracking at a still smaller scale around their own process zones. Whatever the scale of observation (e.g., in the field or laboratory), slip can be viewed as an hierarchical process, in which the larger movements are the result of cracking at successively smaller scales. It is therefore expected that large-scale movements can ultimately be described by processes operating contemporaneously at the smallest scale (Kilburn, 2003(Kilburn, , 2012.

QUANTIFYING REGIMES OF DEFORMATION Parent Relations for Bulk Failure
Bulk failure begins when fracture growth becomes selfperpetuating within part of the stressed crust. It may occur when the failure stress, S F , is first achieved, or after an extended interval of inelastic deformation under a stress maintained at S F , during which the externally supplied stress replaces the stress lost by fracturing and fault movement, so resembling rock creep (Figure 4).
Applying well-known methods from classical statistical mechanics (Reif, 1985;Ruhla, 1992;Guénault, 1995), the magnitudes of local stresses are expected to follow a Boltzmann distribution about the bulk value, for which inelastic deformation ε in increases with supplied stress S sup as (Kilburn, 2003(Kilburn, , 2012: where the activation stress S act is the additional stress required to initiate bulk failure, the characteristic stress S ch is the atomic free energy per volume available for deformation, (dε in /dS sup ) a is the rate at which natural fluctuations in atomic configuration attempt to initiate failure, and exp (−S act /S ch ) measures the probability that an attempt is successful (Reif, 1985;Ruhla, 1992;Guénault, 1995). The activation stress is S F −S, the difference between the failure and applied differential stresses. In the quasi-elastic regime, it decreases to zero as S F is approached and all local attempts to fracture are successful [dε in /dS sup = (dε in /dS sup ) a ]. Fracturing continues to accelerate in the inelastic regime, owing to the local redistribution of stress around fractures (Lawn, 1993;Valkó and Economides, 1995). The redistribution concentrates stress at fracture tips and, as long as the bulk stress is maintained at S F , increases the stress remaining after each increment of growth. A decreasing proportion of S sup is required to continue fracture growth, so that the interval between growth steps persistently decreases, leading to runaway growth when the interval becomes infinitesimally small.
If S R is the mean amount of stress redistributed around fractures, then, from the viewpoint of fracture tips, the activation energy becomes (S F − S R )−S = (S F − S R )−S F = − S R . The negative sign indicates the local surplus of stress at fracture tips (which is balanced by a local deficit where stress has been transferred from around the sides of a fracture). Although the failure strength is not changed, it appears to be decreasing by an amount equal to S R and so mimics fracturing promoted by progressive rock weakening. Interpretations that have qualitatively related accelerated fracturing to stress corrosion at fracture tips (Anderson and Grew, 1977;Atkinson, 1984;Main and Meredith, 1991;Lockner, 1993;McGuire and Kilburn, 1997;Kilburn and Voight, 1998;Kilburn, 2003) may thus alternatively be viewed in terms of increasing stress concentration around fractures in rock with unchanged strength. Setting ε in proportional to the number N of VT events, Eq. (3) yields for the quasi-elastic regime: and for the inelastic regime: Equations (4 and 5) are the basic expressions for deriving relations between VT event rate and both deformation and time. A first comparison with Eq. (3) suggests that the probability of bulk failure is proportional to exp [(S−S F )/S ch ] and to exp ( S R /S ch ) in the quasi-elastic and inelastic regimes, respectively. The exponential terms, however, do not represent the same type of probability. The quasi-elastic regime describes conditions before the start of bulk failure, so that exp [(S−S F )/S ch ] measures the probability that bulk failure will begin under an applied differential stress S. As a result, the probability becomes one when S = S F . The inelastic regime, in contrast, describes when the failure process will be completed, so that exp ( S R /S ch ) measures the probability that failure will occur within a given time interval.

The Characteristic Stress
The value of the characteristic stress in Eqs (4 and 5) changes with the conditions of loading. Independent of the applied stress, atoms deform naturally as they oscillate around their equilibrium positions. The oscillations may randomly yield configurations that assist stress-induced deformation and so promote failure at a stress smaller than S F . The total free energy per volume for oscillations is S * = (3φT+P c −P p )/3, where T is absolute temperature (K), P c and P p are the confining and porefluid pressures, and the molecular energy per unit volume per temperature, φ, has a notional value of (7 ± 1) × 10 4 J m −3 K −1 for common silicate rocks (Kilburn, 2012). How much of the total free energy is utilized determines the characteristic stress. Failure in compression is limited by shearing between atoms, but in tension by the tensile failure of bonds. Shearing requires the integrated deformation of all bonds and so utilizes the full amount of S * . Tension, however, requires deformation only among bonds in a particular direction and so utilizes a fraction of S * ; the amount utilized in tension defines the rock's tensile strength, σ T .
The condition that S ch = S * in compression has been verified against laboratory experiments (Kilburn, 2012). It was also applied to VT precursors of the 1991 eruption of Pinatubo, in the FIGURE 4 | (A) The amount of deformation due to fault movement is given by the difference between the equilibrium (unbroken curve) and ideal elastic (broken line) deformation trends. As the differential stress increases to its failure value S F , behavior evolves from elastic (white), through quasi elastic (yellow) to inelastic (magenta). This yields (B) accelerating quasi-elastic and linear inelastic increases in VT event with deformation. (C) For a constant rate of supplied stress, the deformation rate remains constant until a hyperbolic increase during inelastic behavior. The hyperbolic increase may be preceded by an interval of steady inelastic VT rate (orange). (D) The contemporaneous VT event rate increases from exponentially to hyperbolically with time, again possibly separated by an interval at constant rate (see Figure 2). (E) The model VT-deformation trends are consistent with normalized pre-eruptive sequences observed at El Hierro, Canary Islands, 18 July-12 October 2011 [squares; data from Istituto Geográfica Nacional [IGN] (2011) and Sagiya (2011)], Rabaul, Papua New Guinea, 1971-1994data from Robertson and Kilburn (2016)], and Mauna Ulu on Kilauea, Hawaii, 14 November 1971-04 February 1972data from Bell and Kilburn (2011)]. See also Figure 3. Measures of ground movement (λ; displacement and tilt) and their characteristic values (λ ch ) are used as proxies for deformation. In all three cases, inelastic behavior (magenta) began when λ/λ ch ≈ 4 and eruptions occurred (violet symbols) when λ/λ ch ≥ 4. The trends were normalized using the total number of VT events at the end of the quasi-elastic regime, N QE,m , and the characteristic ground movement λ ch . See Figure 5 for the non-normalized trends. Typical magnitude ranges for VT events were 0.5-2.0 for Rabaul, 1.0-2.5 for El Hierro and 1.5-4.0 for Mauna Ulu. ]. The inelastic regime shows linear trends. The data are consistent with transition between regimes (star) when the ratio of total deformation to characteristic value (λ/λ ch ) is about 4. See Figure 4 for the normalized trends.
Philippines, and 1995 eruption of Soufriere Hills, on Montserrat (Kilburn, 2003), implying the unlikely scenario of eruptions through crust in compression. Extensional stresses are instead anticipated in crust being deformed by a pressurizing magma body and, as demonstrated below, the precursory sequences at Pinatubo and Soufriere Hills can be explained more simply by equating S ch instead with the tensile strength.

Changes in Seismicity With Deformation
Equations (4 and 5) describe the change from an exponential to linear increase in VT number with deformation, corresponding to the evolution from quasi-elastic to inelastic behavior. In the quasi-elastic regime, S supp ≈ S [Eq. (2)] and both S sup and S can be set to Eε, the product of Young's modulus and bulk deformation. If ground movement, λ, is proportional to the bulk deformation, Eq. (4) yields exponential increases with ground movement of both the change in VT seismicity and the total number of VT events [ Table 1, Eqs (T1) and (T2)]. In the inelastic regime, the additional elastic strain supplied by S sup is consumed in faulting (Section 3). Changes in N and λ now both measure VT events, so dN/dλ is constant and total VT number increases linearly with ground movement [Figures 4, 5 and Table 1, and Eq. (T3)].
The transition from quasi-elastic to inelastic behavior occurs when the applied differential stress reaches its failure value, S F . Since S ch = σ T for tensile deformation, the transitional value for S/S ch in Eq. (4) is expected to be S F /σ T for crust being stretched. Applying the Mohr-Coulomb-Griffith criterion for bulk failure, S F /σ T is four or less in tension, between 4 and 5.6 in extension (combined tension and shear) and 5.6 or more in compression (Figure 6; Secor, 1965;Shaw, 1980;Sibson, 1998). In the quasielastic regime S is proportional to ground movement, so that S F /σ T can be obtained from λ/λ ch at the end of the regime. The examples in Figures 4, 5 are consistent with values for λ/λ ch of about four and, hence, with bulk failure in tension.

Quasi-Elastic Changes in Seismicity and Deformation With Time
The VT-deformation trends describe how the proportion of inelastic movement varies with total deformation. To transform the trends into time series, independent constraints must be introduced to specify how loading changes with time. The simplest condition is a constant rate of stress supply. Under this condition, Eq. (4) for quasi-elastic behavior leads to exponential increases with time in VT rate and the total number of VT events [ Table 1, Eqs (T4)-(T5)]. The contemporaneous rate of ground movement is constant because S sup is proportional to deformation [ Table 1, Eq. (T8)]. Quasi-elastic behavior at a constant deformation rate thus naturally yields the exponential VT-time trends frequently observed at closed volcanoes (Figures 2, 3; Shepherd et al., 2002;Bell and Kilburn, 2011;Kilburn, 2012;Wall, 2014).
Under a constant deformation rate, the ratio λ/λ ch is measured by t/t ch , where t ch is the characteristic timescale of the exponential VT time series. If T is the duration of the full exponential sequence, then T/t ch = S F /σ T . Observations at closed volcanoes suggest preferred values for T/t ch of 2-4 ( Figure 6). The range is consistent with the limiting values for λ/λ ch obtained from the VT-deformation trends (Figures 4, 5), so reinforcing the interpretation of deformation in tension (Figure 6).
The particular value of S F /σ T in tension changes with the geometry of failure. The generic Mohr-Coulomb-Griffith criterion applies to failure along a plane (e.g., creating a new failure plane or pulling opening a sealed fault) and values for S F /σ T are ≤4. Pressurized bodies, in contrast, rupture their margins at values of S F /σ T that change according to their shape; for example, S F /σ T is three for a sphere, but two for a long vertical cylinder (Jaeger, 1969;Saada, 2009). The preferred field values of 2-4 are thus consistent with the onset of bulk failure at the margins of magma bodies, as well as the opening of healed faults in crust being stretched by those bodies. In all cases, tensile bulk failure requires the effective principal stresses (normal stress−pore-fluid pressure) to be less than three times the tensile strength (Figure 6). For a notional strength of 10 MPa, therefore, the effective principal stress cannot exceed about 30 MPa, which corresponds to maximum lithostatic depths of about 1.2 and 2 km in dry and water-saturated crust. Tensile failure at greater depths thus implies deformation of super-saturated crust with pore-fluid pressures greater than hydrostatic (Shaw, 1980).

Inelastic Changes in Seismicity and Deformation With Time
In the inelastic regime, a constant rate of stress supply generates hyperbolic increases in VT and deformation rates with time [ Table 1; Eq (T9)]. After a time τ, the rates tend to an infinite value, which is taken to be the mathematical equivalent of continuous fracture coalescence and successful bulk failure (Voight, 1988;Kilburn, 2003); τ therefore defines the duration of accelerating inelastic deformation, which is analogous to tertiary brittle creep under a constant external stress.
The hyperbolic trends are equivalent to linear decreases with time in the inverse VT and deformation rates. Applied to VT seismicity, principal episodes of fracture coalescence coincide with local peak rates (dN p /dt) -or local minima in inverse-rates [(dN p /dt) −1 ] -during which the mean increase in stress concentration per event is B(S F 2 /E), where E is Young's modulus and the geometric constant B = π for straight-edged fractures (Kilburn, 2003). Hence S R = B(S F 2 /E)N p and so, after manipulation, the parent relation for inelastic deformation [Eq. (5)] yields dN p /dt = (dN p /dt) a exp [B(S F 2 /E)N p /σ T ] for a constant rate of stress supply in extension. This form of relation can be re-expressed to show changes in VT event rate with time (Voight, 1988), leading to: or, alternatively, to: (7) where the subscript 0,IN denotes values at the start of the inelastic regime.
The gradient for the inverse-rate VT minima in Eq. (6) has a magnitude γ * = B(S F /σ T ) 2 (σ T /E) and units of "per VT event" ( Table 1). For common crustal rocks, σ T /E∼10 −4 (Mogi, 2007;Heap et al., 2009) and so, with B = π and S F /σ T between 2 and 4 (Figure 6), the gradient is expected to be ∼10 −3 −10 −2 per event, which is consistent with the values of 6.4 × 10 −3 and 2.4 × 10 −3 B, Geometric factor for cracks growing in the inelastic regime; E, Young's modulus; L, Reference length for determining mean deformation; N, Number of VT events; N d , Number of detected VT events; N p , Value of local peaks in the number of VT events; S, Mean maintained differential stress; S ch , Characteristic stress or free energy per unit volume (in extension, S ch is the tensile strength, σ T ); S F , Mean failure stress; S sup , Mean supplied differential stress; S R , Mean stress redistributed around fractures; t, time; t ch , characteristic time for exponential trends, S ch /(dS sup /dt); β, Mean ground movement per VT event; ε, Mean strain; γ * , B(S F /σ T ) 2 (σ T /E); λ, Measure of ground movement as proxy for mean deformation; per event found for inverse-rate minima at Pinatubo in 1991 and Soufriere Hills, Montserrat, in 1995 (Figure 7; Kilburn, 2003). In Eq. (7), the duration of inelastic acceleration is given by τ = [γ * (dN p /dt)] −1 [so when (t−t 0 , IN ) = τ, the inverse-VT rate tends to zero and the VT rate tends to infinite values]. At andesitic strato-volcanoes, observed values of dN p /dt are ∼10 events per day, which, with typical values of γ * , give the observed durations on the order of 10 days (Kilburn, 2003). Hence, if about 30-50% of the total time is needed to confirm a hyperbolic increase in rate, realistic warning times at such volcanoes are on the order of days.

ELASTIC-BRITTLE FAILURE AND VOIGHT'S RELATION
Equations (4 and 5) quantitatively capture the essential features of VT-deformation trends before eruptions. Their good agreement with field data confirms the mutual dependence of VT and deformation precursors and suggests that (a) precursory deformation evolves from quasi-elastic to inelastic, (b) when the rate of stress supply is constant, the increase in VT event rate with time changes from exponential (quasi-elastic behavior) to hyperbolic (inelastic behavior), and (c) precursory behavior can be interpreted in terms of elastic-brittle deformation, without the need to invoke additional rheological responses in the crust (such as plastic flow); this does not preclude the operation of additional processes, only that they need not be invoked unless field data indicate otherwise.
For a constant rate of stress supply, Equations (4 and 5) yield time series (Eqs. T4 and T6) that can be re-expressed as which has the form of Voight's original FFM Relation [Eq.
(1)]. The exponential and hyperbolic VT rates coincide with Voight's relation for α = 1 and 2, but it has yet to be confirmed whether intermediate values of α have any physical meaning. Indeed, it is possible that intermediate values are artifacts from seeking best-fit single trends to partial data sets that extend across the quasi-elastic and inelastic regimes. Comparison with Eq. (1) also shows that the term A in Voight's FFM Relation is not a constant, but changes from 1/t ch for α = 1 to γ * for α = 2; in other words, 1/A is a characteristic timescale when α = 1, but a characteristic number of VT events when α = 2.
Applied to contemporaneous changes in rates of ground movement with time, an expression similar to Eq. (8) is available only for inelastic behavior, when the rate increases hyperbolically (α = 2) to yield d 2 λ/dt 2 = (γ * /β)dλ/dt (where β is the mean ground movement per VT event). Hence, the deformation and VT rates both favor α = 2 immediately before bulk failure, a feature noted in Voight's original analysis (Voight, 1988). During quasi-elastic unrest, however, the deformation rate remains constant while the VT event rate increases and so Voight's relation is not equally applicable to the two precursory signals.
When the rate of stress supply varies with time, the VT and deformation time series derived from Eqs (4 and 5) no longer follow Voight's relation. Hence, the fact that the FFM has  Secor (1965); Shaw (1980), and Sibson (1998)]. The Mohr circles (continuous and dashed semi-circles) describe states of stress between the maximum and minimum effective normal stresses; positive values are in compression, negative values in tension. The diameter of a circle gives the applied differential stress, S. Failure occurs when the minimum normal stress meets the failure envelope. Tensile failure corresponds to where the circles meet the envelope at zero shear stress (hence intersect the horizontal axis) and the minimum effective normal stress equals the tensile strength (with a normalized value of -1). To satisfy this condition, the Mohr circle must have a diameter less than or equal to 4σ T (or normalized diameter of four; arrows beneath graph) The maximum value (larger continuous semi-circle) may represent the opening of a failure plane. For comparison, failure at the margins of a long vertical cylinder is shown with a normalized diameter of two (smaller continuous semi-circle). The dashed semi-circle (normalized diameter of 5.6) shows the Mohr circle at the transition between failure in extensional (tension and shear) and compression (shearing alone). (B) Values of T/t ch for quasi-elastic pre-eruptive behavior at Etna, Sicily (blue squares; Wall, 2014), Kilauea, Hawaii (green triangles; Bell and Kilburn, 2011), Mauna Loa, Hawaii (magenta diamond;Lengliné et al., 2008), Piton de La Fournaise, Réunion (small dark red triangle; Lengliné et al., 2008) and Soufriere Hills, Montserrat (black diamond; Figure 2). Most of the values for T/t ch lie between 2 and 4 (lower and upper dashed lines).
such sequences may often develop under a constant rate of stress supply.
Pre-eruptive sequences at strato-volcanoes commonly have durations of years or less. Longer sequences provide a greater opportunity for significant variations in the rates of stress supply, in particular the decadal time intervals for precursory sequences at large calderas, such as Rabaul (Robertson and Kilburn, 2016) and Campi Flegrei (Kilburn et al., 2017). For example, the full 23-year sequence of unrest before Rabaul's 1994 eruption did not show changes with time in rates of VT event and of uplift consistent with Voight's relation -even though it followed the trend between VT events and deformation expected for an evolution from quasi-elastic to inelastic behavior (Figures 4, 8). Thus the VT event rate and uplift rate both fluctuated by an order of magnitude or more (∼10 2 −10 4 events per month and c. 0.015−0.4 m y −1 ), and the peak rates occurred half-way through the sequence (Figure 8). During its final 2 years, however, the sequence did show changes in rates that were consistent with hyperbolic increases with time (Figure 8).
The peak rates in the Rabaul sequence have been associated with the arrival of magma at a depth of about 2 km, possibly intruding into an existing magma chamber (McKee et al., 1984). The contemporaneous inverse VT rates showed a less clearly defined trend (Robertson and Kilburn, 2016).
The final rates have instead been associated with segments of the ring fault being torn open to allow magma to erupt (Robertson and Kilburn, 2016). The full sequence thus represents the caldera being stretched to breaking point under a variable rate of increasing magmatic pressure. A constant rate of stress supply was established during the final 2 years of deformation, when the crust had already entered the inelastic regime, so promoting the hyperbolic increase in rates with time (described by the Voight's relation with α = 2). However, rates of stress supply were not constant for the first 20 years of unrest, during which time Voight's relation could not be used to describe the precursory time series.

DEVIATIONS FROM MODEL ELASTIC-BRITTLE CONDITIONS
By describing common trends among VT and deformation precursors, Eqs (4 and 5) provide a starting point also for identifying when the model cannot be applied without adjusting its underlying assumptions. For example, the model assumes that the crust's mean behavior follows that of an elastic medium containing a dispersed population of small faults. It does not explicitly accommodate conditions when bulk deformation is controlled by movements of faults with lengths similar to that of the crust being deformed. Slip along such faults would favor the occurrence of a small number of large-magnitude VT earthquakes instead of a large number of small VT events; as a result, VT trends may be dominated by only a few earthquakes, whose total number is too small to yield repeatable mean behavior.
Even in the presence of large faults, VT seismicity may remain controlled by small earthquakes. For example, before Rabaul caldera's 1994 eruption, VT seismicity was concentrated in an annulus related to ring faults extending to depths of 2-4 km and produced by caldera collapse about 1,400 years beforehand (Nairn et al., 1995;Jones and Stewart, 1997). Most of the VT events had magnitudes of 0.5-2, triggered by the slip of small faults in the crust surrounding the ring faults, rather than movements along the ring faults themselves (Robertson and Kilburn, 2016). Similar behavior has been argued for VT seismicity in the vicinity of ring faults at the Campi Flegrei caldera in southern Italy (Troise et al., 2003). In both cases, the large ring faults appear to have constrained the size and geometry of the crust being deformed, rather than to have contributed directly to rates of VT seismicity.
The initial application of Eqs (4 and 5) implicitly assumes that self-accelerating crack growth begins at the start of the inelastic regime. In this case, precursory VT event rates are anticipated to accelerate with time across the transition from quasi-elastic to elastic regimes. Such behavior is indeed observed in the field (Bell and Kilburn, 2011) and laboratory (Kilburn, 2012). However, the transition may also be marked by an interval of constant VT event rate between the quasielastic and inelastic accelerations (Figure 2). Such steady rates may develop in the inelastic regime when the coalescence of fractures, which favors accelerating rates, is temporarily retarded by stress barriers in rock between non-interacting fractures (Main, 2000;Heap et al., 2009;Girard et al., 2010). An accelerating transition may thus represent the limiting case when fracture coalescence dominates from the start of the inelastic regime. Hence, should a steady rate appear after quasi-elastic acceleration, it cannot be taken to indicate an approach to stability, but must be viewed as a temporary lengthening of unrest until the onset of hyperbolic accelerations to bulk failure (Figure 2).
Although bulk failure in the crust is necessary before eruption at a closed volcano, it does not ensure that an eruption will occur (Kilburn, 2003;Grasso and Zaliapin, 2004). For example, similar precursory sequences have been identified before eruptions and intrusions at Kilauea (Bell and Kilburn, 2011) and at Krafla in Iceland (Blake and Cortés, 2018). Failure may be initiated at the margin of a magma body, after which an eruption requires a magmatic overpressure large enough to drive magma to the surface (Jellinek and DePaolo, 2003). It may also begin at a more distant location where stresses are concentrated (e.g., the tips of a major fault) or where the crust is locally weak (e.g., altered rocks in a hydrothermal system), in which case the propagating fracture must additionally extend to the magma body itself. The probability of eruption can thus be expressed more generally as the product of the probability of bulk failure, the probability that failure breaches a magma body, and the probability that magma can erupt through the new breach. The elastic-brittle model addresses only the first of these and so, by itself, is a necessary but insufficient condition for guaranteeing an eruption.
In common with most analyses of VT and deformation precursors, the elastic-brittle model does not specify conditions within the magmatic, or other, systems driving volcanic unrest. One exception is the model by Lengliné et al. (2008), who investigate how specific conditions of magma accumulation may control rates of precursory signal at basaltic volcanoes. The acceleration in VT number with time may also increase with the size of the following eruption (Schmid and Grasso, 2012). Future models that couple crustal stresses with magmatic processes promise to yield new insights for refining and enhancing forecasting procedures.

SCALE-INDEPENDENT FRACTURING
The elastic-brittle model is inherently scale-independent, because it expresses rates of VT seismicity and deformation in terms of the ratio of mean applied differential stress to characteristic stress. Applied stress measures changes in the free energy per volume available for atoms to do work (in this case to deform and break bonds); the characteristic stress measures the total free energy per volume already available for a specific type of deformation (e.g., tension or compression) before a differential stress is applied. Thus, if e is the change in free energy and e the free energy available, then S/S ch = ( e/V)/(e /V), where V is the volume of deformation. At the atomic scale, V is the volume of an atom itself. At the macroscopic scale, it refers to the size of rock being deformed, including the volume of discontinuities. At any given scale, the same V is used to calculate both stresses from their respective energy terms, so that S/S ch always measures the ratio of change in free energy to the reference free energy available. Hence, even though the values of individual stress terms may change with scale and amount of fracturing, their ratios remain the same. It assumes a transition from exponential to hyperbolic VT event rates after an amount of ground movement 4h ch and, for a constant rate of stress supply, after a time 4t ch .
The influence of fracturing is shown by the change from quasi-elastic to inelastic deformation regime, which appears to correspond to a change from weak to strong interactions among fractures (Vasseur et al., 2017) and favors the onset of self-accelerating crack growth. Exploring the interactions has generated an independent class of fracture models, most of which have focused on strong interactions immediately before bulk failure (Main, 1996(Main, , 2000Molchan et al., 1997;Guarino et al., 1998;Amitrano et al., 2005;Girard et al., 2010). Applying statistical methodologies, they have shown that increases in seismic event rate with time in the inelastic regime are expected to follow a power-law trend of the form: where t 0 is the starting time, τ f is the duration of the sequence before bulk failure and p is an undetermined empirical factor. The result is scale-independent because it depends on the interaction among fractures, regardless of size, and values of p between 0.5 and 2 have been obtained from field data (Amitrano et al., 2005;Schmid and Grasso, 2012). When p = 1, the result is identical to that obtained for the acceleration to failure in the inelastic regime [Eq. (7)] and corresponds to the limit of α = 2 in Voight's analysis (Voight, 1988). The similarity of Eqs (7) and (9) suggests agreement between the different modeling strategies. Their future integration may thus yield enhanced procedures for general forecasts of bulk failure.

PRECURSORS AS INDICATORS OF ENERGY LOSS BEFORE ERUPTIONS
In addition to demonstrating the elastic-brittle features of unrest, precursory time series indicate preferred energy states in magmatic systems before eruption. Common quasi-elastic behavior suggests a preference for steady rates of deformation, which are a natural consequence of minimizing rates of energy dissipation. Primary controls on energy dissipation are VT seismicity in the crust and frictional resistance during the transport of magma. Under a constant deformation rate, the VT event rate increases exponentially with time and so is not an obvious minimizing factor. A constant deformation rate, however, is favored by a constant rate of pressurization in the magmatic source and this, in turn, is promoted by a constant flux of magma from depth. Approximately steady deformation rates may thus reflect a preference to minimize rates of energy loss during magma transport.
FIGURE 10 | The elastic-brittle trends can be transformed into probabilities of bulk failure and levels of alert. The example here assumes a steady rate of stress supply and an immediate transition from exponential to hyperbolic VT event rates after a time 4t ch . Inverse trends are plotted because changes from an exponential to a linear decrease is easier to visualize that a change from exponential to hyperbolic increase. The transition from quasi-elastic to inelastic regime marks a step increase in alert level. Each level may itself consist of gradational changes, as implied by the colored gradient. The probability being quantified changes from the probability that bulk failure begins in the quasi-elastic regime to the probability that failure will be complete before the time interval τ in the inelastic regime.
A second remarkable feature is the agreement between theory and observed inverse-VT rate gradients in the inelastic regime. The model VT expressions implicitly assume that the number of detected events is similar to the number of essential events for bulk failure. In principle, however, fault movement may be triggered in parts of the crust that do not directly affect the volume that will ultimately fail, so that the number of detected events is greater than the essential number. The difference will not affect expressions for VT events in the quasi-elastic regime, because it changes only pre-exponential terms and so cancels from both sides of the relevant equations. In the inelastic regime, however, N appears within the exponent itself [Eq. (T6)], so that the inverse-rate field gradient is given by γ * N es /N -that is, the model gradient multiplied by the ratio of the number of essential (N es ) to detected (N) VT events -which is ≤γ * . At Pinatubo and Soufriere Hills, Montserrat, the observed gradients are within the expected theoretical range, suggesting that the proportion of non-essential VT events is small. Although such a condition may not be universal, the fact that it appears for these two examples implies that precursory deformation may tend to involve the minimum amount of crust possible for initiating bulk failure, which again indicates a preference for minimizing energy loss.

POTENTIAL APPLICATION TO OPERATIONAL FORECASTS
Event trees have become a popular aid for deciding levels of alert during a volcanic emergency (Newhall and Hoblitt, 2002;Aspinall, 2006;Sandri et al., 2009;Marzocchi and Bebbington, 2012;Selva et al., 2012;Sobradelo et al., 2013Sobradelo et al., , 2014. They pose a sequence of interconnected questions about interpreting precursory signals, the combined answers to which are used to estimate the probability of eruption. Critical pre-eruptive values for individual signals (e.g., number of VT events and amount of ground deformation) are identified empirically from observations of previous eruptions at the same volcano or at apparently analogous volcanoes elsewhere (Sobradelo and Martí, 2015).
The progression from quasi-elastic to inelastic behavior yields objective pre-eruptive criteria that complement the existing empirical evaluations. This is especially important at volcanoes for which no data are available from previous unrest. An example of how additional procedures may be applied is shown in Figure 9 and summarized below. It is based on the simplest conditions for a constant rate of stress supply and no steady VT rate between quasi-elastic and inelastic regimes.
1. Identify the regime of deformation. Quasi-elastic behavior will show an exponential increase in VT number with ground movement and yield a characteristic movement, λ ch . VT events and deformation will follow distinctly different time series. For a constant rate of stress supply, the deformation rate is constant and the VT rate will increase exponentially with time, yielding the characteristic timescale t ch ; following the VT rate here is crucial, because a constant deformation rate on its own may erroneously be interpreted as a sign of dynamic stability. Inelastic behavior will show a constant increase in VT number with ground movement and VT event rate maxima and deformation rate will both increase hyperbolically with time (equivalent to linear decreases with time in inverse event rate for VT minima and inverse deformation rate). 2. If behavior is quasi-elastic, the transition to inelastic behavior is expected when the total amount of ground movement is between 2λ ch and 4λ ch and, if deformation rate is constant, after corresponding times of 2t ch and 4t ch . 3. If behavior is inelastic, the magnitudes of the inverse-rate gradients will be in the range 10 −3 -10 −2 for VT minima. Larger gradients may indicate that the linear trends are spurious. Extrapolating the trends to zero inverse-rate (i.e., infinite rate) provides the preferred time τ for bulk failure to be completed.
In all cases, observations must cover a minimum range of values to confirm that a trend is not spurious. For exponential trends, the minimum range is one characteristic interval (λ ch or t ch ), because shorter observations may yield trends statistically indistinguishable from linear; at least one quarter of a trend's total duration must therefore be used to identify its existence. For the inverse-rate trends, the smallest range is a decrease in inverserate minima by a factor of two (hence, at least half of a trend's total duration) before a linear trend can be proposed.
Under the ideal conditions when unrest data follow a complete sequence from lithostatic equilibrium (without a differential stress) to bulk failure, the elastic-brittle model provides estimates of the maximum ground movement (2λ ch -4λ ch ) before the emergence of inelastic deformation and the onset of bulk failure (and, for a constant rate of stress supply, also the maximum time to the transition, 2t ch -4t ch ). For a continuous acceleration in VT rate across regimes (without at intervening interval at steady rate), bulk failure is expected after an additional time τ at most.
The probability that quasi-elastic behavior will lead to bulk failure increases exponentially with the amount of ground movement (and with time for a constant deformation rate); once it has started, the probability that failure will be completed before an additional interval τ increases hyperbolically with time ( Figure 10). The ends of the quasi-elastic and inelastic regimes thus coincide with critical changes in the potential for eruption. They are natural stages for step-like increases in alert (Figure 10) and can readily be incorporated into existing forecasting event trees (Selva et al., 2012;Sobradelo et al., 2013Sobradelo et al., , 2014 or measures of unrest (Potter et al., 2015).
At volcanoes reawakening after long repose, emergency measurements are often gathered systematically a significant time after the start of unrest, so that the transition from quasielastic to elastic behavior will occur sooner than the ideal maximum time. In the least-favorable case, measurements may become available only when the crust is already in the inelastic regime.

CONCLUSION
The good agreement between field observation and parent model for changes in seismicity with deformation at closed volcanoes suggests that their behavior before eruptions is normally governed by the stretching of elastic-brittle crust. Immediate implications are that there is no starting requirement to invoke deformation mechanisms other than elastic-brittle, and that VT signals reflect how a population of small faults responds to changes in stress.
The trend between VT events and deformation persists independently from how seismicity and ground movement vary with time, as shown most clearly by the case for Rabaul (Figure 8). A corollary is that VT and deformation time series are mutually dependent and controlled by the rate at which the crust is being stressed. In particular, they reveal how a constant rate of stress (e.g., from magma pressurization) will yield the observed combinations of exponential VT event rate with constant deformation rate, and of hyperbolic VT event rate with hyperbolic deformation rate. These combinations show that the popular FFM model of Voight (1988) is a particular form of the parent VTdeformation model under a constant rate of stress supply. They also provide new criteria to be tested for integrating deterministic physical constraints into probabilistic forecasts of eruption.
The parent model, however, is still not a complete description of pre-eruptive conditions at closed volcanoes. It identifies conditions for bulk failure in the crust, which, although necessary to open a new pathway for magma ascent, do not guarantee that magma will reach the surface. It also focuses on conditions that favor continued acceleration in VT seismicity with time, without significant intervals of steady VT rate. The parent model thus provides a starting point for identifying additional precursory trends and preeruptive criteria. Outstanding goals include identifying whether precursory signals can distinguish between pre-eruptive and non-eruptive outcomes; whether VT rates will accelerate to bulk failure without an interval of steady behavior; and, indeed, whether final VT accelerations are inevitable a significant time before eruption. Achieving these goals will provide new constraints for coupling changes in crustal stresses with specific magmatic processes and, hence, yield greater insights for refining and enhancing current forecasting procedures.

AUTHOR CONTRIBUTIONS
CK developed the elastic-brittle model for precursors to eruption and prepared the manuscript.

ACKNOWLEDGMENTS
Careful reviews by Jérémie Vasseur, Jean-Robert Grasso, Raffaello Cioni and Valerio Acocella clarified and improved the original manuscript. The research was privately funded. Thanks are also due to Barry Voight, whose support and advice during the emergency at Soufriere Hills, Montserrat, in 1996 laid the foundations for the present work.