Pre‐Eruptive Damage, Weakening and Magma‐Edifice Coupling at Piton De La Fournaise Volcano

Eruptions in basaltic volcanoes are often preceded by increasing seismicity and surface deformation, which progressively damage and weaken the volcanic edifice. We show how damage and crack interaction produce the inverse Omori‐Utsu law for earthquakes during pre‐eruptive periods. Rock mass continuity, representing damage, is shown to decrease exponentially with the earthquake number; we interpret it as a general form of the Omori‐Utsu law. Pre‐eruptive earthquake time series are shown to be controlled by heterogeneity distribution, finite‐size effect and crack interaction, and by the feeding system characteristic time. Magma‐edifice coupling is described by state variables that depend on the continuity and the feeding system characteristic time. Pre‐eruptive seismicity of the 2004–2017 24 summit/proximal eruptions of Piton de la Fournaise volcano was well modeled by an inverse Omori‐Utsu law. It allowed identifying two cases: (a) strong crack interaction and earthquake number acceleration, when failure in strong intact rock and finite‐size effects dominate the brittle fracture process; in that case the magma‐edifice interaction power exhibits a maximum before the eruption; (b) weak crack interaction, generating an almost constant earthquake rate and corresponding to a brittle fracture process at constant strain in a weak, fractured rock mass. In this latter case eruptions occurred when the continuity reached a critical value, close to 0.25. Specific times are identified, from the time variations of the state variables; they define estimators that provide values of the eruption time within 10% of the true value in 60%–75% of the cases studied, from the complete time series.


Introduction
During the pre-eruptive phase, volcanoes often show precursory increase in deformation and seismicity, sometimes at mid/long term.This particularity could be used to anticipate or forecast the eruption.This has led to an intense research in this domain, often devoted to forecast the eruption time (Bell, Greenhough, et al., 2011;Bell, Naylor, et al., 2011;Bell et al., 2021;Boué et al., 2016;Chastin & Main, 2003;Collombet et al., 2003;Kilburn, 2018;Kilburn & Voight, 1998;Main, 1999;McGuire & Kilburn, 1997;Schmid et al., 2012;Voight, 1988; see also Vasseur et al., 2017) from long (10-100 days) pre-eruptive time series.In practice, attention is often paid to the occurrence of earthquake swarms (Roman & Cashman, 2006), which often immediately precede eruptions.Duputel et al. (2019), for example, show that volcano-tectonic earthquake swarms at 1.5-2 km depth have most often preceded the aseismic dyke propagation and the eruptions by 1-5 hr, at Piton de la Fournaise volcano, from 2014 to 2018.Such earthquake swarms reveal the final step of the brittle fracture process, when the volcanic edifice becomes locally out-of-equilibrium.However, similar swarms also occur without being immediately followed by eruptions.This stresses the fact that the exact time of the final rupture and eruption is controlled by medium heterogeneities and details of the stress field, this latter being controlled by the magma pressure and rheology, and by the edifice weight, topography and rheology.Medium heterogeneities and details of the stress field can be represented only by their (unknown) probability density, and rupture and rupture time can be known only in terms of probability.Tang (1997), Tang and Kaiser (1998) and Tang and Hudson (2011) have shown the relation between medium heterogeneity (described by a Weibull function), stress-strain curve and acoustic emission.Long-term deformation and seismicity therefore contain information that may reveal when the conditions for final rupture are met; such an information can be used to characterize and model the state in which the volcano is, relative to these conditions.This mid-term characterization may be thought to be less sensitive to small-scale heterogeneities, the effect of which will be found in the model residual.Such a mid-term characterization is the aim of this work.
In this frame, a volcano may be understood as a system described by state variables that evolve with time, due to the coupling between magma pressure and edifice response.It is thus necessary to study, at the scale of the volcano, the constitutive laws of the magma and the edifice, and to link them to geophysical observables.Preeruptive increase in volcano deformation and seismicity may be, for example, linked to the pressurization of the magma feeding system, and to the progressive fracturing and weakening of the edifice.Our work will be first devoted to establish firmly the relation between seismicity, deformation and weakening of the elastic properties of the edifice, through a damage approach.In a second part, the effect of this weakening on the coupling with the pressurized magma, especially on the state variables (pressure, magma flow) will be investigated.It is, however, not a complete modeling of the pre-eruptive process, since the evolution of magma properties with time or deformation is not modeled, nor constrained by using geophysical observables in this study.
Joint inversion of seismicity and deformation data was used to infer the propagation of the 2007 Father's Day intrusion in Kilauea volcano (Segall et al., 2013), and of the 2014 Bardarbunga dike (Heimisson & Segall, 2020), by using Dieterich (1994) seismicity-rate law.More recently, Bell et al. (2021) carefully studied the uplift and seismicity driven by magmatic inflation at Sierra Negra volcano, which provided perfect pre-eruptive time series from 2006 ro 2018.Damage approach was used to model the propagation of a dyke pressurized by gas-saturated magma (Maimon et al., 2012) and allowed to distinguish two (fracture-and magma-) controlled regimes.This approach was also used to model hydraulic stimulation (see, e.g., Lyakhovsky & Shalev, 2021;Shalev & Lyakhovsky, 2013).Carrier et al. (2015) and Got et al. (2017) modeled the February-March 2007 Piton de la Fournaise and the 2004-2011 Grimsvötn pre-eruptive surface deformations, by using a seismicity-based damage model.The 2004-2011 Grimsvötn pre-eruptive earthquake series has been shown to follow an inverse Omori law, what has allowed finding an analytical solution for the magma overpressure, surface displacement, magma flow and magma-edifice interaction power.However, pre-eruptive earthquake time series do not generally follow a simple inverse Omori law and this simple analytical model has to be adapted to more frequent cases.Piton de la Fournaise volcano is characterized by a shallow pressurized magma reservoir (Peltier et al., 2007) so that eruptions are often preceded by intense deformation and seismicity (see, e.g., Duputel et al., 2019).In this paper, using the 2004-2017 seismicity and deformation data from Piton de la Fournaise, we have studied the earthquake response of the edifice, especially the relation between seismicity, damage, earthquake interactions and stress transfer that produce inverse Omori-Utsu laws.Understanding these relations and laws is necessary to interpret earthquake time series, and their relations with the pre-eruptive dynamics.We also studied the changes induced in the magma-edifice coupling by the occurrence of earthquakes, and we evidenced three main stages in the volcano deformation.We have identified the key parameters controlling this coupling, and expressed the state variables (magma overpressure and flow, interaction power) as a function of these parameters; we have searched for estimators of the eruption time.This work is not focused on eruption forecasting, where eruption time is inferred from the knowledge of a first part only of the pre-eruptive time series.Instead, we have used complete pre-eruptive deformation and earthquake time series, up to the eruption time.Results of the estimations have been compared to the real eruption times and used to infer probability density functions for these estimators, to quantify the modeling error.

Data
The data are constituted by both the volcano-tectonic earthquakes and surface displacements time series recorded by the Observatoire Volcanologique du Piton de la Fournaise (OVPF) seismic and GPS networks (Figure 1 and Figure S1 in Supporting Information S1).The earthquake time series begins in 1998 (Figures S2-S4 in Supporting Information S1), and continuous GPS surface displacement records begin in 2003, so that 34 pre-eruptive  (Peltier et al., 2007).earthquake times series and 27 GPS pre-eruptive time series were available.The first complete (GPS and seismicity) pre-eruptive time series was that of 2 May 2004.In order to keep the same geometrical conditions during this study, we limited the eruption set to summit and proximal eruptions, occurring in the summit cone or at its base after a direct vertical path generating volcano-tectonic earthquakes in the summit volume located between the shallow-level reservoir and the surface, excluding distal eruptions occurring in the rift zones after a horizontal path.Finally data from 24 summit/proximal eruptions were used, between 2004 and 2017 (Figure 2).These summit/proximal eruptions break the same portion of the edifice; their rupture area, located between the reservoir and the summit cone (Figure 1 and Figure S1 in Supporting Information S1), and the corresponding damage and continuity parameters defined in the model section are thus comparable.Most (>95%) of the pre-eruptive earthquakes related to these eruptions and recorded by the seismic network are located (Duputel et al., 2021) within 1 km above the sea-level magma reservoir inferred from geodetic displacement analysis (Peltier et al., 2007).Earthquake magnitudes are carefully computed since 2012, and completeness magnitude is 0.8 since this time (Figure S5 in Supporting Information S1).
Earthquake counts and Global Navigation Satellite System (GNSS) time series issued from the OVPF analysis were represented on a daily basis.GNSS solutions are processed by using the GAMIT/GLOBK (Herring et al., 2018) postprocessing software package.

Model Setting and Physical Justifications
The model used in this study has been inspired by the seismicity and surface displacement data recorded at Piton de la Fournaise volcano and some other basaltic volcanoes (Figures 3 and 4; see, e.g., Bell et al., 2021;Got et al., 2017).These data show a progressive acceleration of the recorded earthquake number, and a simultaneous remarkable pattern for the surface displacement, comprising a first rapid inflation stage, with a decreasing rate with time, followed by a constant inflation rate stage, and then by a terminal stage showing an increasing inflation rate.The first rapid inflation stage is often interpreted as due to the pressurization of an out-of-equilibrium shallow-level reservoir; during this stage the local earthquake rate is close to the background seismicity rate, whereas it increases during the two following stages.
Displacement data recorded close to the summit of a volcano do not allow to discriminate the detailed geometry of the magma feeding system.Increasing rates of earthquakes occurring during pre-eruptive periods shows that the medium is stressed above its linear elastic limit; it is fractured and therefore cannot be considered merely as linear elastic.As a consequence, we will represent our model with the simplest possible geometry, acceptable for a volcano, and a damage approach.This choice does not exclude that more complex geometries (see, e.g., Kozono, 2021;Reverso et al., 2014) may exist; these geometries are not constrained by the seismicity and displacement data used, but they are not incompatible with damage models.
In this study we therefore use a model represented by a pressurized spherical magma reservoir fed by a cylindrical vertical magma conduit in a damageable elastic half-space (Figure 5; see also Carrier et al., 2015;Got et al., 2017;Lengliné et al., 2008).Got et al. (2017) have shown that this geometry was a member of a broad class of models, described by the same equations regardless of the precise description of their geometry, which is represented by the characteristic time.These models are damageable Kelvin visco-elastic models: Kelvin visco-elastic models with elastic modules decreasing with fracturing.We use an effective medium approach: fracturing evolution is not described geometrically, but using the earthquake number and the corresponding displacement data; fracturing increases the volume of the pressurized magma reservoir by decreasing the shear and Young's modules of the volcanic edifice, according to the earthquake and displacement data.This model is axisymmetric, homogeneous and isothermal.Shear modulus, however, decreases with the earthquake number and surface displacement following a damage approach described in the next sections.This decrease accounts for the elastic potential energy consumed in fracturing.The recorded earthquake number is an input data of the model; during the pre-eruptive deformation process, it is controlled by both the stress field and the medium heterogeneity.The time evolution of the model is therefore constrained by the medium heterogeneity, and not only by the visco-elastic coupling of the edifice with the feeding system.The model is constrained by earthquake and displacement data; we will see that it depends on three parameters only: the linear elastic limit displacement u el , the characteristic time of the magma feeding system τ 0 (controlled by the magma viscosity and feeding system geometry), and the proportion δ of the total area to be ruptured, ruptured by each earthquake, whatever this rupture is pre-, co-, or post-seismic, and whatever the orientation of the rupture plane.
The magma pressure is considered as constant at the base of the conduit.This condition is fulfilled when a large magma chamber collects mantle magma through a sufficiently large area, in such a way the magma inflow in the chamber through this area is equal to the magma outflow through the conduit.The magma pressure in the chamber  Bell et al., 2021).The GPS station is located at the center of the caldera.cannot drop strongly, to ensure chamber stability.The shallow-level reservoir is initially under-pressurized, and the pressure gradient generates a magma flow from the magma chamber.The magma feeding system is initially closed by the surrounding rock mass that is pressurized and progressively fractured in tension, by the magma reservoir overpressure.Rock mass is damaged by fracturing; its strength exhibits a peak, and post-peak weakening.Eruption corresponds to the rupture of the rock mass; this rupture is a post-peak process (see, e.g., Li, 1987, Figure 9.13).
Magma is considered incompressible.At Piton de la Fournaise, viscosity of the basaltic magma is low and the edifice exhibits heavily fractured and gaspermeable rift zones, so that gas content and magma compressibility are considered low during the pre-eruptive process.Bulk modulus of degassed basaltic magmas is around 10 GPa (Rivalta & Segall, 2008;Spera, 2000) (compressibility 10 10 Pa 1 ) and may be compared to the bulk and shear moduli of fractured rock masses in volcanic areas, which are in the range 0.1-1 GPa (see, e.g., Dzurisin, 2007, p. 281;Heap et al., 2020;Rivalta & Segall, 2008).The incompressible-magma approximation holds as far as the bulk modulus of the magma is greater than the post-peak weakening rate of the rock mass, at the edifice scale; this is the case when the post-peak rupture process remains quasi-static (Figure 6).At Piton de la Fournaise, the beginning of the dynamic process (fast aseismic magma propagation) has been shown to take place less than 12 hr before the eruption, between 2014 and 2018 (Duputel et al., 2019), whereas the time scale of our study is 10-100 days, during which the deformation process is quasi-static.
Our model is a semi-coupled model, in which the solid edifice properties evolve with deformation and damage, whereas the magma properties remain constant.The fully coupled problem, in which the magma viscosity and compressibility simultaneously evolve due to the coupling with the deforming and damaging edifice, remains to be solved.
We consider that there is no leak nor addition of magma outside the conduit so that at any time the magma flow (given by Poiseuille's law, left member of Equation 1) is equal to the rate of linear elastic volume variation of the magma reservoir (given by, e.g., Delaney & McTigue, 1994, right member of Equation 1; see Figure 5 for the meaning of the symbols): ΔP(t) is the magma overpressure in the reservoir as a function of the time t and P is the overpressure when the equilibrium is reached, that is when the volume of the reservoir no longer varies, in the linear elastic case.Note that this equation means that the edifice loading is controlled by the strain rate (magma flow), not by the stress rate.
In that case magma reservoir overpressure is controlled by a simple equation (see the complete derivation of the equations in Section A in Supporting Information S1): where Figure 5. Physical model used in this study (from Carrier et al., 2015;Got et al., 2017;Lengliné et al., 2008).A magma reservoir (radius a r ) embedded in a homogeneous isotropic elastic half-space (shear modulus G) is fed by magma through a cylindrical conduit (radius a c , length H c ). Shear modulus G is assumed to decrease homogeneously with damage (see text for details) and therefore with time t.Pressure P s at the base of the conduit is assumed to be constant.The magma is characterized by its viscosity μ.In the text, the feeding system will be represented by its characteristic time where Journal of Geophysical Research: Solid Earth 10.1029/2023JB027595 GOT ET AL.
The characteristic time τ(t) results from both the feeding (characteristic time τ 0 ) and the damage process, represented by the Kachanov's continuity , where S(t), E(t), G(t), and D(t) are respectively the intact area, the Young's modulus, the shear modulus and the damage parameter at time t and S 0 = S (t = 0), and ψ(t = 0) = 1.Magma reservoir overpressure, surface displacement u(t) and strain ϵ(t) are related by (Got et al., 2017): (where u el and ϵ el are respectively the limit displacement and limit strain reached at the equilibrium when the edifice is considered linear elastic).Using Equations 2-4 we can find that which is the strain form of the equation governing a damageable Kelvin model (Got et al., 2017), with τ 0 = η E 0 , where η is the viscosity in the Kelvin model.

A Solution for the Time Evolution of Normalized Overpressure, Surface Displacement, Strain, Magma Flow and Power Only Depending on τ 0 and ψ(t)
The solution of the homogeneous equation associated to Equation 2is Determining the constant C by the Euler-Lagrange method, we find If we write the magma reservoir overpressure writes The displacement u(t) and strain ϵ(t) are Notice that Equation 11 may be rewritten, in a compact form: (a1, green line): case where the magma pressure drops more sharply with strain than the edifice strength; the edifice is unloaded during the post-peak deformation, which remains a quasi-static incremental process.The magenta lines correspond to the visco-elastic re-loading of the edifice by the magma reservoir pressure.(a2, red line): case where the magma pressure drops less with strain than the edifice strength.An instability initiates (beginning of the red line), which leads to the dynamic macro-rupture of the edifice and dyke propagation (see, e.g., Duputel et al., 2019).(b) Black line: stress-strain curve of a weak, fractured, volcanic edifice.In that case the probability that magma pressure drops more sharply with strain than the edifice strength (case b1, orange line) is high, and the post-peak deformation of the edifice remains quasi-static.
As quoted by Got et al. (2017), magma flow may be derived from Poiseuille's law under the form: so that we can define (Got et al., 2017) a normalized magma-edifice interaction power: where Π 0 = PQ 0 .
Equations 9-12 constitute a general analytical solution of Equation 2. This solution allows expressing the system state variables ( ΔP P , Q Q 0 , Π Π 0 ) as a function of the characteristic time τ 0 of the feeding system and the continuity ψ(t) only.The linear elastic displacement u el is a scaling parameter.Notice that the characteristic time τ 0 controls the viscous component of the model and ψ the brittle or plastic component.In the following section, we will see that ψ (t) can be computed directly from the event number N(t), so that this solution (Equations 9-12) doesn't need any seismicity model; the model parameters τ 0 , δ and u el are computed by fitting u(t) to the displacement data.The solutions (Equations 9-14) are identical to the Runge-Kutta solution of Equation 2 (Figure S7 in Supporting Information S1).

Expression of the Continuity ψ(t) as a Function of the Earthquake Number N(t)
Changes in the material elastic properties are contained in the damage or continuity evolution law (see, e.g.Costin, 1987;Kachanov, 1958), and determine the reaction of the volcanic edifice.The Kachanov's definition of the continuity in continuum damage mechanics is based on the tensile rupture of metal rods of finite, known, section S. In this work we will use summit eruptions that break the same, finite and known, volume of the edifice, located between the shallow-level magma reservoir and the surface; the state of stress in this volume is tensile.The use of these eruptions warrants the continuity, as the area to be broken, to be well-defined.
We reproduce the progressive failure of rocks by using Amitrano and Helmstetter (2006)'s incremental elastic damage approach, in which the intact area S(t) and the Young's modulus E(t) decrease for each earthquake: S (t + dt) = S(t) dS(t) = (1 δ)S(t) where δ = dS S is the relative incremental decrease of the intact area S, for each earthquake; δ is taken constant as the earthquake magnitude distribution is dominated by magnitudes in a narrow range (see also Carrier et al., 2015;Got et al., 2017).As a result S(t) = (1 δ) N(t) S 0 where N(t) is the total number of earthquakes that have occurred up to the time t and S 0 is the initial, intact, area.The same reasoning may be used for E(t) (Amitrano & Helmstetter, 2006).It leads to write the continuity: is the proportion of rupture area created at the time t by N(t) earthquakes.Since δ is computed by fitting the displacement data by Equation 9knowing N(t), when the completeness magnitude is larger, N(t) is smaller but δ is larger.It is always possible to find a value of δ such that N(t)δ (and the continuity ψ(t) = e N(t)δ ) remains constant, if N(t) varies by a constant fraction during the pre-eruptive period.Continuity is thus a robust quantity with respect to changes in completeness magnitude.
In the incremental form, the continuity function is the solution of the damage evolution equation (by taking the time derivative of Equation 15): Journal of Geophysical Research: Solid Earth 10.1029/2023JB027595 GOT ET AL.

An Example of Time Variation for the Model State Variables
To illustrate the solution found in Equations 9-14 and the effect of various time distributions of earthquakes on the continuity and the state variables of the system, we choose to represent the earthquake number as a function of time by an inverse Omori-Utsu law of parameter p′ (see, e.g., Schmid & Grasso, 2012): where n 0 = n (t = 0), and p′ varies in the range [0 1.5] in order to represent with the same law, a constant earthquake rate (p′ = 0) and various increasing rates.
Figure 7 shows the family of solutions, parametered by p′, which gives the time variations of the model state variables.Equations 9-14 may therefore be used with a variety of continuity functions computed from the data (Equation 15); these equations constitute a generalization of the results formerly obtained by Got et al. (2017).
Figure 7 shows that N(t) increases and ψ(t) decreases more or less abruptly, depending on the value of p′; the displacement u(t) reproduces qualitatively the 3-stage features evidenced on Figures 3 and 4, which inspired the model.The magma overpressure drops more or less strongly, at the end of the process, according to the value of p′.This non-linear pressure drop is due to the non-linear increase of the reservoir volume; volume variation per unit time is equal to the magma flow (Equation 13 ).This increase in volume is due to damage, which induces a decrease in shear modulus.The characteristic time given by 1 τ 0 e N(t)δ + n(t)δ (Equations 3 and 15) shows the competition, in the pre-eruptive dynamics, between the magma feeding process (first term of 1 τ(t) , which dominates at low damage level, during the second stage) and the damage process (second term, which dominates at high damage level, during the third stage); this competition results in a maximum for τ(t).The magma-edifice interaction power (Equation 14) shows that a second maximum may appear for high enough values of p′; this interaction power increases with magma flow, but it drops with the edifice rock strength.

Results and Discussion
A careful analysis of the 34 1998-2017 inter-eruptive earthquake number time series shows that they often begin with a constant earthquake rate larger than the background seismicity rate; 22 of them (i.e., 63%) exhibit a preeruptive acceleration for at least 2 days (Figure 2 for 2004-2017, Figure S8 and Table S1 in Supporting Information S1).However, these earthquake time series show a remarkable variability, ranging from a linear variation to a very sudden acceleration.
Analyzing jointly earthquake and displacement time series, taking into account the proposed Kelvin damageable visco-elastic model (Equations 2-5), leads us to describe qualitatively the pre-eruptive dynamics in three main stages (Figure S9 in Supporting Information S1): -primary creep stage, during which the overpressure in the shallow-level reservoir increases close to the elastic limit -no earthquake occurs in the vicinity of the pressurized magma reservoir.During this stage the edifice remains linear elastic, and creep refers to the visco-elastic behavior of the system (reservoir + edifice) during this stage.It could reach an equilibrium if the linear elastic limit was not reached.-Secondary creep stage occurs when the overpressure exceeds the linear elastic limit and produces rupture area, earthquakes and surface displacement at a constant rate.-Tertiary creep stage occurs when the earthquake number and the surface displacement accelerate.
In the last two stages, creep is the brittle creep (see, e.g.Amitrano & Helmstetter, 2006;Brantut et al., 2013;Heap et al., 2009) and refers to the effect of the progressive damage of rocks in the edifice.
Understanding this first qualitative interpretation requires a thorough understanding of the progressive damage process: how damage increases with earthquakes, how earthquakes interact, how stress transfers and produces an inverse Omori-Utsu law, how can we interpret the parameters n 0 and p′ of this law, and what importance they have for the pre-eruptive dynamics and volcanic system state variables (Sections 4.1-4.6).

Understanding the Relation Between Earthquake Number and Continuity Function
Continuity ψ(t) = 1 D(t) is directly related to the damage parameter D(t).It is an essential function that critically controls the time evolution of material strength; with the characteristic time τ 0 , continuity controls the reservoir overpressure and surface displacements (Equations 9-12).Its form (Equation 15) must therefore be carefully assessed, and compared to similar quantities proposed by previous authors.Weibull (1939) described the material failure during tensile test experiments as a stochastic process controlled by the material defects and their distribution; the survival probability of the material is the product of the survival probabilities of the defects in a given volume.The probability product accounts for the serial arrangement of cracks loaded in tension (the weakest-link theory), whereas the integration of survival probabilities over the volume accounts for the parallel arrangement of cracks (see, e.g., Zok, 2017 for a short discussion, Batdorf & Crose, 1974 who give a proof for Weibull theory in the case of tensile stress states, and Weiss et al., 2014 who investigate compressive stress states and limitations of Weibull theory).He defined a failure function (damage function; see also Krajcinovic & Silva, 1982) at the sample scale from the probability to reach the strength σ m ; its cumulative distribution is: where , and σ is the applied stress and v the volume; m is in the range 1-20, and termed the homogeneity index (Tang, 1997;Tang & Kaiser, 1998;Weibull, 1939;Wilkins, 1980).The distribution of the continuous parts in the material 1 F(σ) is sometimes called the stretched exponential cumulative distribution in materials science and can be compared to the continuity (Equation 15).Coleman (1956Coleman ( , 1958) ) proposed a description of the material failure in terms of the progressive failure of a fiber bundle loaded in tension (the Fiber Bundle Model, or FBM), by expressing the rate of fiber failure using the breakdown rule (see also Krajcinovic, 1996;Nanjo, 2017;Rundle et al., 2003;Turcotte & Glasscoe, 2004;Turcotte et al., 2003), or evolution equation where n f (t) is the number of unbroken fibers at time t, and ν(σ) the hazard rate, a function of the fiber stress σ(t).
The quantity n f 0 where n f0 is the initial number of unbroken fibers may be considered as a discrete representation of the material continuity in the continuum damage approach.
Equations 16 and 19 can be considered as identical if we consider that the earthquake rate is a measure of the hazard rate.The result of the integration of Coleman's Equation 19 is directly comparable to the continuity function 1 F(σ) deduced from Weibull's Equation 18 if we consider that n m (σ) represents ν(σ).Therefore we can consider, as a first approach, the convergence of these reasonings as a validation of our expression of the continuity as a function of the earthquake number (Equation 15).This expression has to be linked with other existing relations and approaches used, especially in volcanology.The damage evolution or kinetic equation of Kachanov (1958Kachanov ( , 1986) may be written (where ϵ el is the linear elastic strain, A K and n K are positive constants), from which the Voight (1988) evolution equation can be derived by using the change of variable ψ = dt dΩ , with A V = A K ϵ el n K and β = n K + 2, positive constants.
Equation 21 is the basis of the Failure Forecast Method (FFM), used for eruption prediction.Equation 20 leads to where t c = (A K (n K + 1)ϵ el n K ) 1 is the rupture time.This solution for the continuity corresponds to Equation 15 ψ (t) = e N(t)δ when N(t) follows an inverse Omori law N(t) = n 0 t c ln(1 t t c ) , when n 0 t c δ = 1 1+n K = a.Notice that the Kachanov-Voight approach implies that critical size and time exist, whereas Coleman (1956) breakdown rule and Equation 15 do not.Das and Scholz (1981) and Main and Meredith (1989) 15when N(t) is the inverse Omori law.In the following sections we will show in which mechanical conditions this latter law appears.

Linking Crack Interaction and Continuity
At this point, it is of importance to understand the effect of the interaction between cracks in the form taken by damage and continuity (see, e.g., Amitrano et al., 1999;Costin, 1985Costin, , 1987;;Dahm & Becker, 1998;Reches & Lockner, 1994;Main 2000;Rivalta & Dahm, 2004), and in the earthquake number N(t).
Computing the effective Young's modulus E′ of a cracked elastic solid, Walsh (1965) has shown that where E is the Young's modulus and v is the volume of the solid, v is the volume of a spherical crack, and N the number of cracks.N v v is the crack density.
This reasoning is performed for one isolated crack then extended to N cracks by using the superposition theorem for linear elasticity.This is actually the case when the cracks are weakly interacting, that is, when the crack density N v v ≪ 1.Under this condition, one can write: where ψ wi is the continuity in the weakly interacting crack case (Budiansky & O'Connell, 1976, Self-Consistent Method).However, this reasoning cannot be generalized to the case where N v v tends toward 1, where cracks can strongly interact.
To take into account the interaction between cracks, Bruner (1976) proposed to introduce the cracks one at a time in an equivalent homogeneous elastic solid having the effective moduli due to the introduction of the previous cracks.He found that, without any condition on N v v : where ψ si is the continuity in the strongly interacting crack case.Equation 24 appears to be a first-order approximation of Equation 25(see also Cox & Meredith, 1993 for a discussion).Equation 25is similar to Equation 15, which may be derived from the evolution Equation 16.
Progressive damage and earthquake interaction are the physical justification of our Equation 15, which is directly related to Bruner (1976)'s crack interaction theory, to Coleman (1956Coleman ( , 1958))'s FBM load transfer schemes, and finally, to heterogeneity distribution (see, e.g., Weibull, 1939).This also shows that N(t) is the true control parameter of the continuity, and not directly the time as in Kachanov (1958).

Consequences for the Inverse Omori-Utsu Law
When failure occurs by progressive damage close to the limit equilibrium, the action, or applied force F is equal to the reaction, or strength, of the solid so that we can write: where σ 0 is the stress applied on the initially intact, finite area S 0 , and σ is the Kachanov effective stress, actually applied on the undamaged area S.This area S decreases with the damage function D: S = (1 D)S 0 = ψS 0 and the Kachanov effective stress increases when S decreases: the stress initially applied on the damaged area ΔS transfers on S at each rupture, increasing the probability of earthquakes and failure.
After Weibull (1939) in various materials, Charles (1958), Elliott (1958) and Mould and Southwick (1959) in glass, and Cruden (1970) (which contains an extended discussion of brittle creep in rocks), Cruden (1974), Wilkins (1980), Atkinson and Meredith (1981), Lajtai and Schmidtke (1986) and Lajtai et al. (1991) in rocks have established a widely used static fatigue or brittle creep law in various materials linking time t f and stress σ at rupture, that can be written (see, e.g., Anderson &Grew, 1977, andHeap et al., 2009 more recently, for a review): where t f0 is the characteristic time of the failure occurring at stress σ 0 , and α is found to be in the range 10-100 for intact rock samples (see, e.g., Amitrano & Helmstetter, 2006;Wilkins, 1980).This relation is sometimes referred as Charles (1958)'s law.Gupta (1982) has shown the equivalence of Charles (1958)'s power law and Wiederhorn and Bolz (1970)'s exponential law of slow crack growth.A parallel can also be drawn with Glen (1952)'s law for ice, and with Herschel-Bulkley's power law for non-newtonian visco-plastic fluids (see, e.g., Bonn et al., 2017;Nanjo, 2017).
Failure time t f can be estimated as the average time between two ruptures dt dN , which is the inverse of the earthquake rate or probability of occurrence n(t).In that case, Equation 28 can be written: Considering that the stress inducing earthquakes is the effective stress applied on undamaged areas (Equation 27), Equation 29 can be written: The earthquake rate n(t) therefore results from the adaptation of the effective stress field to the continuity (itself controlled by the distribution of material defects, i.e., the crack density), in response to an external stress forcing.
Equation 30 can be used to derive a general earthquake production law during the damage process: The left side of Equation 31 is controlled by the stress field time evolution, and the right side is controlled by the actual material crack distribution and interaction.

Weak Crack Interaction Case:
The Inverse Omori-Utsu Law With Low p′ ( p′ < 0.5) In the weakly interacting crack case, continuity is a simple linear function of N(t) (Equation 24) and Equation 32may be written: which is a first-order linear differential equation, the solution of which is: α+1) and N c = N (t = t c ).The corresponding earthquake rate n(t) = dN(t) dt is that is, an inverse Omori-Utsu law with a positive parameter p′ = α α+1 < 1.As a consequence, the inverse Omori-Utsu law with p′ < 1 observed in the pre-eruptive earthquake time series is the result of quasi-static damage (Equation 26) by weakly interacting cracks; crack interaction and finite-size or free-surface effect grow with p′.This reasoning and result can be compared with those obtained by Nanjo (2017).Reaching complete failure with weakly interacting cracks suggests subcritical crack growth, occurring at low stress threshold (the long-term strength) and low stress transfer with a few finite-size effect in weak, eventually fractured volumes, the fracturing process as a whole being a long-term process.This case is similar to the low-m case in Figures 3a-3e of Tang andKaiser (1998) or in Figures 4.10-4.11 of Tang andHudson (2011), to the compliant layer case in Gudmundsson (2005), to the loose-sand case in geomechanics (see, e.g., Nicot et al., 2023), or to the stable mode of Lyakhovsky and Shalev (2021).The end-member case where there is no crack interaction, that is p′ = 0, corresponds to the failure at constant displacement or strain described, for example, in Turcotte et al. (2003) and Gudmundsson (2012), which is possible in already fractured media; therefore n 0 is the event rate without crack interaction.Crack interaction and finite-size effect require a minimum rock strength to appear.
From Equation 34and p′ = α α+1 < 1, we can write that N c = n 0 t c 1 p′ and that may be written where ln p′ is the Tsallis q-logarithm with q = p′ (Tsallis, 1988), or where n 0 t c is the event number occurring without crack interaction, and ψ 0 (t) = 1 t t c is the continuity function in the same hypothesis (p′ = 0); the logarithm (or exponential) function accounts for the interactions, which is coherent with Bruner's incremental approach.Therefore, the inverse Omori-Utsu law may be written using the Tsallis nonadditive q-entropy (Tsallis, 1988) S p′ (for q = p′), if the total number of possible configurations or energy (actually strength) states is W = 1 ψ 0 , all states being equally probable:

Strong Crack Interaction Case:
The Inverse Omori-Utsu Law With High p′ ( p′ ≥ 0.5) When p′ tends to 1, n(t) and the interactions between cracks increase, which requires sufficient rock strength; the characteristic length of the stress field perturbation reaches the size of the body to be ruptured.This may be compared to the high-m case in Figures 3a'-3e' of Tang and Kaiser (1998) Nicot et al., 2023), or to the runaway mode of Lyakhovsky and Shalev (2021).In the strong crack interaction case, the continuity is expressed by Equations 15 or 25, and Equation 32 writes (see, e.g., Turcotte et al., 2003 for a similar reasoning): which is a first-order linear differential equation, the solution of which is: if N (t = 0) = 0 and t c = 1 n 0 αδ .This corresponds to a simple (p′ = 1) inverse Omori law n(t) = n 0 1 t tc , and may be written As for p′ ≠ 1, the logarithm function accounts for the interactions; Equation 42 corresponds to Equation 38when p′ = 1.Notice that 1 α = n 0 δt c = a, which means that in the inverse Omori law case, the stress field variation, controlled by α, is exactly adapted to the distribution of defects, controlled by n 0 δt c .Equation 42 may be written: where ψ Omori (t) = (1 t t c ) a (see Got et al., 2017).
Equation 43 may be compared to the Boltzmann law for entropy (whose Tsallis nonadditive q-entropy is a generalization), where the continuity is the inverse of the possible number of states (combinations of elements that are ruptured or not) W on the rupture area.Basaran and Nie (2004), Chan et al. (2012), Basaran et al. (2014) and Tucker et al. (2014) express the damage function using a decreasing exponential of entropy (see also Bryant et al., 2008 anda complete review in Osara &Bryant, 2019).Osara and Bryant (2019) and Idris et al. (2020) found that the normalized entropy generation during a fatigue degradation process was well approximated by the normalized number of loading cycles.
In Equation 28, α is a positive finite constant so that 0 < p′ < 1, and the earthquake number N(t) in Equation 34reaches a finite value N c .The limit case p′ → 1 implies a large α, a very small t f (Equation 28) that is a large n(t), and a small δ to have a physical t c = 1 n 0 αδ ; it means that the complete, progressive, quasi-static rupture of an intact, finite, area occurs by means of a very large number of very small, frequent and interacting earthquakes.At each time step, the applied stress is exactly equal to the rock mass strength (the limit equilibrium); in this case, crack interaction is the maximum possible at equilibrium.The inverse Omori law (p′ = 1) appears to be an upper bound, for a given n 0 t c , for earthquake production in strong and limited rock masses during quasi-static crack propagation.If the rupture occurs exactly at the peak rock mass strength, or if the rock mass has a plastic behavior, this corresponds to the constant load case described in Turcotte et al. (2003) and Gudmundsson (2012).If the rock mass weakens with strain, the maximum strength of the rock mass and the applied stress (and pressure) can decrease while remaining equal, maintaining conditions similar to the constant stress condition.
From this analysis it appears that the case where p′ > 1 occurs when the action (magma overpressure) is larger than the reaction (strength of the edifice), when the edifice rock mass weakens more strongly with time and strain than magma pressure.In this case the process is no longer quasi-static damage: the difference between action and reaction produces a strong acceleration and drives the rupture front, and the earthquake production.
Notice also that averaging the pre-eruptive earthquake number time series with various p′ (see, e.g., Chastin & Main, 2003;Collombet et al., 2003;Schmid et al., 2012) reveals the persistent physical feature (the free-surface or finite-size effect) and cancels the variable part due to changing rock mass strength and magma pressure.
From this discussion, we conclude that the expression ψ(t) = e N(t)δ (Equation 15) of the continuity function is physically well-grounded for rock masses loaded in tension, and may be used to model the pre-eruptive mechanical weakening of volcanic edifice rock masses, in particular of their elastic moduli.The solution of the Kachanov-Voight evolution equation is a particular case of this Equation 15when N(t) follows an inverse Omori law.

Consequences for the Understanding of the Secondary and Tertiary Creep Stages During Pre-Eruptive Periods
From the latter analysis we can complete our qualitative description of the earthquake and surface displacement time series in terms of creep stages or phases (Figure S9 in Supporting Information S1): -Secondary (brittle) creep stage occurs when the overpressure exceeds the elastic limit and produces sub-critical crack propagation, earthquakes and surface displacement at a constant rate.Crack interaction is therefore limited; no finite-size or free-surface effect appears.-Tertiary (brittle) creep stage occurs when the earthquake number and the surface displacement accelerate, that is when crack interaction becomes strong.Finite-size or free-surface effect appears, that is, condition for crack propagation is met at the scale of the edifice; crack interaction is maximum up to the surface, stress perturbation size and correlation length between cracks reach the edifice characteristic size, cracks can propagate rapidly up to the surface (see, e.g., Rivalta & Dahm, 2006 and the references therein) and an eruption may occur.This stage ends with the earthquake swarms and eruptions described in Duputel et al. (2019), who evidenced short and intense swarms before the short (1-5 hr) rapid (0.1-0.3 m/s) aseismic propagation of the rupture and magma immediately before the eruption.

Analysis of the Parameters of the Inverse Omori-Utsu Law
We fitted the 24 summit/proximal pre-eruptive earthquake time series using Equation 34, searching for n 0 and p′ (Figure S8 and Table S1 in Supporting Information S1).Most (20) of the eruptions were preceded by 1,000-4,000 earthquakes so no stacking was required to reveal an inverse Omori-Utsu law.The fits were performed by adjusting the secondary creep phase controlling n 0 , and the end of the tertiary creep phase that controls p′; they give n 0 and p′ for the homogeneous distribution of cracks equivalent to the real medium.
Results may be represented by two different linear laws (Figure 8): if p′ < 1, and if p′ > 1.The two lines intersect at p′ ≈ 1, for n 0 t c N c ≈ 0.1275.
From Equation 36 we deduce that, when p′ < 1, where a = n 0 t c δ; in the theoretical case where N c δ = 1 (no overlapping of earthquake sources), p′ = 1 a; p′N c δ represents the fraction of the rupture area created with crack interaction, and a = n 0 t c δ the one without crack interaction.
Equation 44 shows that the actual p′ is 1.15 times larger than expected from Equation 36; it results, for a given t c , in a larger N c (and area to be ruptured) than theoretically expected.This may be due to the overlapping of earthquakes sources, to non-zero residual stresses and/or to healing.
Figure 8 expresses the partition between (a) the stable crack propagation in weak, fractured, rock masses with a low stress transfer (p′ < 0.5 and large n 0 ), and (b) the final fast crack propagation in resistant rock masses with a strong crack interaction and stress transfer (high p′ values).This partition is inherent to the inverse Omori-Utsu law for a given N c (Equation 36); it recalls the two regimes evidenced by Maimon et al. (2012), or by, for example, Lyakhovsky and Shalev (2021) in hydraulic fracturing ("stable" or "runaway").
In the Supporting Information S1 Section B we computed, for each eruption, the quantities N c δ and a = n 0 t c δ.We show that the parameter a (Figure 9) is equivalent to the Monkman-Grant constant (Monkman & Grant, 1956) used in materials sciences.We show (Section B in Supporting Information S1) that this approach allows evidencing an efficient estimator for t c when the Monkman-Grant parameter a is a constant.However our computations show that N c δ increased (Figure S10 in Supporting Information S1) and a decreased (Figure 9) with time from 2004 to 2017, which limits the use of this approach in that case.After 2014 most of the earthquake rupture area was created just before the eruption, which may be an evidence of increasing instability.The parameters a and N c δ should be systematically computed.

Consequences for the Volcanic System State Variables
As shown in the previous sections, the stress transfer dynamics (characterized by p′) controls the seismicity and damage/continuity evolution, and continuity controls the state variables of the volcanic system, with the characteristic time of the magma transfer system τ 0 .Figure 7 shows that, when p′ increases, the state variables may strongly change before the final rupture and eruption.
To investigate the possible real behaviors, system state variables have been computed for 24 summit/proximal Piton de la Fournaise eruptions from 2003 to 2017 .The model parameters τ 0 , u el and δ are estimated by fitting the theoretical displacement u(t) to the displacement data, using the earthquake number N(t).The often high-quality of the fit of the displacement is directly related to the form taken by the continuity function; it may be understood as an a posteriori validation of this form.
The results of the system state variable computations  Information S1) may be represented using two end-members (Sections 4.6.1 and 4.6.2).

The Case of a Clear Acceleration of N(t)
In this case (Figures 10 and 11) p′ > 0.5, crack interaction, tertiary creep and finite size/free-surface effect are strong; it can be compared to the runaway faulting mode of Lyakhovsky and Shalev (2021), where the seismicity is concentrated along the rupture area.The continuity and the overpressure strongly drop down to the rupture time t c .The strong and eventually early drop in overpressure, and its relation with seismicity are consistent with experimental results obtained by, for example, Vinciguerra et al. (2004), Gehne et al. (2019) or Benson et al. (2020) and modeling (e.g., Figure 3 of Shalev & Lyakhovsky, 2013).State variables allow identifying specific times (Figure 12; see also Got et al., 2017) to characterize the evolution of the pre-eruptive damage/brittle creep process and magma-edifice coupling.Especially, the magma-edifice interaction power Π(t) systematically exhibits a second maximum (time t = t GI ) just before the final rupture and subsequent eruption (Section C in Supporting Information S1); this feature may be compared to the result obtained by Main and Meredith (1989) for stress intensity factors.The characteristic time τ(t) = ( ψ(t) τ 0 + n(t)δ) 1 , from Equation 3, drops below τ 0 for 50% of these cases, which means that rupture area production n(t)δ dominates over changes in the visco-elastic response of the edifice at this stage.After the second maximum of Π(t) the stress transfer and crack interaction is very strong; it corresponds to the first, seismic, pre-eruptive phase investigated by Duputel et al. (2019).The final drop of Π(t) and τ(t) show that the reservoir volume increase and the corresponding drop in pressure, due to fracturing (the continuity and shear modulus drop in the model), are faster than the viscous magma flow Q(t).
The time of these maxima can be related to the eruption time t c ; the relative delay t c t GI t c decreases with p′ for p′ > 0.6, which can be described by a linear law (Figure 13): so that an estimator of t c could be:    15) and from δ, estimated from the fit of the computed displacement u(t) (Equation 12) to the displacement data.For these values of p′, Π(t) reaches its second maximum and decreases before the eruption, and the normalized overpressure decreases below 0.5.The characteristic time τ(t) may decrease below τ 0 .
Estimations provided by tc values from the complete time series are within 10% of the true value in 75% of the cases (Figure 13), when t GI is available (60% of the cases).
The analysis of the continuity in the vicinity of the rupture time allows us to evidence that for very strong final crack interactions (p′ > 1.3), final rupture and eruption occur shortly after ψ(t = t ψ 25 ) ≈ 0.25 (Figure 11).For moderate to strong acceleration of N(t) and crack interaction (0.5 < p′ < 1.3), final rupture occurs for larger values of the normalized time delay   S1 in Supporting Information S1 for eruption parameters and caption of Figure 10 for more explanations).The continuity is a convex up, decreasing curve below 0.5, when it drops to 0. The power Π(t) reaches a second, sharp, maximum and decreases abruptly, and the normalized overpressure drops to 0 in a few days before the eruption.The characteristic time τ(t) decreases below τ 0 .Eruption number 23 (red lines) clearly shows the relation between the summit deformation and seismicity, during a pre-eruptive earthquake sequence comprising three successive earthquake swarms.

The Case Where N(t) Is Composed by Several Accelerations
This case corresponds to a piecewise rupture, that is, a series of rapid ruptures of resistant patches producing earthquake swarms that are not followed by eruptions, but rather by intrusions.This behavior is due to the hectometric-scale heterogeneity of the medium (and the corresponding stress perturbation and finite-size effects), which is smaller than the edifice scale.The characteristic time of the rupture of a resistant patch being shorter than that of the feeding system, the magma pressure drops strongly due to the rapid volume creation in the edifice; elastic potential energy accumulated before the rupture of the resistant patch is not sufficient for the magma to reach the surface.Magma pressure then increases again slowly.Weak heterogeneities or layers can deform and create volume laterally at low pressure.This leads to complex/chaotic patterns for N(t), ψ(t), and Π(t) (Figure S13 in Supporting Information S1).This may be due to layering, and may favor dyke arrest (see, e.g., Gudmundsson, 2005); at larger scales, this process may be involved in the creation of magma reservoirs, especially when a weak layer favors the displacement of a volcano flank.In this case t GI exists but provides a strong residual when fitting Equation 47 with t c data and contributes to enlarge the tails of the corresponding pdf (Figure 13).

The Case Where N(t) Decelerates With Time
This is the case after the eruptions occurring between the 2 April 2007 eruption and 2010 (Figure S14 in Supporting Information S1).In that case magma overpressure and Kachanov effective stress slowly decrease by an internal diffusion process in an increasing volume.It may correspond to the creation of a sill structure, for example.The stress perturbation diffuses over an increasing contact area and continuity increases until a new equilibrium is reached.The earthquake number follows an Omori law.

Conclusion
In this work we have proposed a simple semi-coupled model for magma-edifice coupling, for low-compressibility magmas, which allows the computation of pre-eruptive surface displacements, magma overpressure and flow, and magma-edifice interaction power as a function of the continuity of the volcanic edifice and of the characteristic time of the feeding system.Earthquake numbers and surface displacement time series are input data of the model.Model parameters are three: the linear elastic displacement, the characteristic time of the feeding system, and the damage per earthquake; they are estimated by fitting the surface displacement time series.Magma overpressure and flow, and magma-edifice interaction power are computed from the values of these parameters.
Continuity is the intact part of the rock between defects or cracks; it controls the rock elastic parameters.Its time evolution is dominated by the progressive damage and brittle fracture process during pre-eruptive periods.
Convergent reasonings lead us to propose that, at least in the tensile states of stress that we have investigated, continuity is a decreasing exponential of the earthquake number; earthquake number is therefore a logarithm of the continuity, which may be understood as a general form of the Omori-Utsu law similar to Boltzmann-Tsallis between the eruption time t c and the time t ψ 25 for which the continuity ψ(t) is equal to 0.25, as a function of p′ (blue crosses), and model (red line, Equation 48).(b) Sample probability of the relative residual between the eruption time t c and its estimation from t ψ 25 (Equation 49).
entropy.We show how damage and crack interaction produce the inverse Omori-Utsu law for earthquakes during pre-eruptive periods.This law is therefore a characteristic feature of the progressive failure of pressurized volcanic edifices.
The systematic analysis of 24 pre-eruptive earthquake and surface displacement time series preceding summit/ proximal eruptions at Piton de la Fournaise volcano shows that, though it exhibits variability, pre-eruptive earthquake number can indeed be modeled by an inverse Omori-Utsu law.Our study gives physical bases for interpreting the inverse Omori-Utsu parameters n 0 (earthquakes occurring without interaction) and p′ (proportion of interacting earthquakes) and allows identifying two cases: with weak crack interaction (p′ < 0.5).This is the case when fracturing occurs, and eventually magma propagates in a weak, already fractured, poorly cohesive medium.This corresponds to the constant strain boundary condition case of Turcotte et al. (2003) and Gudmundsson (2012), and the stable mode of Lyakhovsky and Shalev (2021).Finite-size/free-surface effect is very limited, and earthquake rate is almost constant; the magma-edifice interaction power Π(t) increases before the eruption, but it never shows a decrease after reaching a maximum.Final rupture and eruption occur when continuity is close to 0.25.with strong crack interaction (p′ ≥ 0.5), when fracturing occurs and eventually magma propagates into strong rock with some defects, so that crack interaction and propagation dominate the progressive damage process.This is the constant force case of Kachanov (1958), Coleman (1958), Turcotte et al. (2003) and Gudmundsson (2012) where the intact rock area decreases, Kachanov's effective stress increases and edifice finite-size effect controls the final rupture process up to the surface; this may also be related to the runaway mode of Lyakhovsky and Shalev (2021).Magma overpressure strongly decreases before the final rupture and eruption, and the power Π(t) reaches a second maximum before the eruption.The case p′ > 1 arises from forced crack interaction, that is, when magma overpressure is larger than the edifice strength.
Specific times have been identified and related to the eruption time t c and the parameter p′, knowing the entire time series.We used them to define estimators of the eruption time.We used the time for which continuity reaches the value of 0.25, which provided an estimator of t c for which 60% of the estimations are within 10% of the true value of t c .We also used the time for which the magma-edifice interaction power reaches its second maximum; when available (60% of the cases), this time provided estimations of t c within 10% of the true eruption time in 75% of the cases.Variability in the damage process induces deviation from the mean value; the result of this study is a probability for t c , knowing the entire time series, for a given value of specific times (t ψ 25 or t GI ).
The Monkman-Grant parameter a = n 0 δt c is the fraction of the total area to be ruptured, ruptured without crack interaction; it may be used to characterize the damage state of the edifice; its time evolution shows that damage and weakening increased from 2003 to 2017 at the summit of Piton de la Fournaise volcano.
We propose to continue this investigation to check the model extensively on other basaltic volcanoes.The approach is easy to implement as a monitoring module; it will be performed in the WebObs-IPGP environment.
Laboratory experimental studies should be used to investigate this model, especially the control of rock heterogeneity distribution on the seismicity and fluid overpressure decrease.We propose to call this approach "DAMAGE."

Figure 1 .
Figure 1.(a) Map of the Piton de la Fournaise volcano showing the OVPF permanent GPS network (black diamonds) and seismic network (blue triangles: stations operating between 2003 and 2011; red triangles: stations added to the network after 2011).Color dots represent the earthquakes recorded by the OVPF seismic network, during the pre-eruptive periods preceding the 24 2004-2017 eruptions used in this study.See the detailed map on Figure S1 in Supporting Information S1.Coordinates are given in the Gauss-Laborde La Reunion projection system.(b) Vertical cross-section of the volume studied in the Piton de la Fournaise volcano summit showing the earthquakes recorded during the studied pre-eruptive periods.Color indicates the normalized occurrence time.Black ellipse represents the magma reservoir inferred from GPS surface displacements inversion(Peltier et al., 2007).

Figure 2 .
Figure 2. Normalized number of earthquakes as a function of the normalized time, the unit time being the duration of the preeruptive period.Color scale indicates the time of the 24 summit/proximal eruptions studied, from the beginning of 2004 to the end of 2017.The yellow curve represents the average of the normalized earthquake numbers; this curve corresponds to an inverse Omori-Utsu law, with a parameter p′ = 1.1.

Figure 3 .
Figure 3. Examples of pre-eruptive dynamics showing the three inflation stages (see text for details).Green circles: earthquake number as a function of time.Blue crosses: surface horizontal displacement recorded at one summit GPS station (SNEG for Piton de la Fournaise, GFUM for Grimsvötn), as a function of time.Red vertical dashed line shows the eruption time.(a) Eruption of 20 July 2006 at Piton de la Fournaise volcano.(b) 2004-2011 pre-eruptive period at Grimsvötn volcano (modified fromGot et al., 2017).

Figure 4 .
Figure 4. 2006-2018 pre-eruptive period at Sierra Negra volcano (modified fromBell et al., 2021).The GPS station is located at the center of the caldera.

Figure 6 .
Figure6.Schematic stress-strain representation of the rupture process (see, e.g.,Li, 1987 for a similar representation) and the magma-edifice coupling in a basaltic volcano, with a special attention to post-peak processes.Pressure in the magma reservoir loads the volcanic edifice, but decreases when the reservoir volume increases; during a quasi-static deformation, the stressstrain curves of the edifice and the magma are almost superimposed.FigureS6in Supporting Information S1 shows the relation with the various deformation phases.(a) Blue line: stress-strain curve of a resistant, cohesive, volcanic edifice.(a1, green line): case where the magma pressure drops more sharply with strain than the edifice strength; the edifice is unloaded during the post-peak deformation, which remains a quasi-static incremental process.The magenta lines correspond to the visco-elastic re-loading of the edifice by the magma reservoir pressure.(a2, red line): case where the magma pressure drops less with strain than the edifice strength.An instability initiates (beginning of the red line), which leads to the dynamic macro-rupture of the edifice and dyke propagation (see, e.g.,Duputel et al., 2019).(b) Black line: stress-strain curve of a weak, fractured, volcanic edifice.In that case the probability that magma pressure drops more sharply with strain than the edifice strength (case b1, orange line) is high, and the post-peak deformation of the edifice remains quasi-static.

Figure 8 .
Figure 8. Inverse Omori-Utsu exponent p′ as a function of n 0 t c N c .Color scale represents N c = N (t = t c ). Straight lines correspond to the models represented by Equation 44, black line, and Equation 45, red line.

Figure 9 .
Figure 9. (a) Earthquake number during secondary creep (no crack interaction) n 0 t c as a function of δ.Color scale represents the time in days from 02/05/2004 (first eruption monitored using GPS measurements and the OVPF seismic network).The red curve shows the model n 0 t c = C δ (Monkman-Grant relation), for C = 0.9.(b) Parameter a = n 0 t c δ as a function of time in days from 02/05/2004 (same color scale as in panel (a)).The red vertical line figures the 30 March 2007 eruption.

Figure 10 .
Figure 10.Model state variables as a function of time, computed for eruptions showing a moderate acceleration in the earthquake number N(t) (0.58 < p′ < 0.7).Numbers in the legend cartoon indicates the eruption number (see Table S1 in Supporting Information S1 for eruption parameters).(a) Cumulative number of earthquakes N(t).(b) Thin line: horizontal displacement recorded at the summit SNEG GPS station; thick line: modeled displacement u(t).(c) Continuity ψ(t).(d) Normalized characteristic time τ(t) τ 0 .(e) Normalized overpressure ΔP(t) P .(f) Normalized power Π Π 0 .The continuity is computed directly from N(t) (Equation15) and from δ, estimated from the fit of the computed displacement u(t) (Equation12) to the displacement data.For these values of p′, Π(t) reaches its second maximum and decreases before the eruption, and the normalized overpressure decreases below 0.5.The characteristic time τ(t) may decrease below τ 0 .

Figure 11 .
Figure11.Model state variables as a function of time, computed for eruptions showing a strong acceleration in the earthquake number N(t) (0.90 < p′ < 1.62, strong crack interaction; see TableS1in Supporting Information S1 for eruption parameters and caption of Figure10for more explanations).The continuity is a convex up, decreasing curve below 0.5, when it drops to 0. The power Π(t) reaches a second, sharp, maximum and decreases abruptly, and the normalized overpressure drops to 0 in a few days before the eruption.The characteristic time τ(t) decreases below τ 0 .Eruption number 23 (red lines) clearly shows the relation between the summit deformation and seismicity, during a pre-eruptive earthquake sequence comprising three successive earthquake swarms.

Figure 12 .
Figure 12.Model state variables as a function of time, computed for the 20 July 2006 eruption, showing the specific times used in the discussion.See Table S1 in Supporting Information S1 for model parameters.(a) Cumulative number of earthquakes N(t).(b) Red line: horizontal displacement recorded at the summit SNEG GPS station; blue line: modeled displacement u(t).(c) Continuity ψ(t).The green dashed vertical line shows the t ψ 25 time.(d) Normalized characteristic time
ΔP 50 time.(f)Normalizedpower Π Π 0 .The green dashed vertical line shows the t GI time.The definition of Π Π 0 (Equation14) shows that t GI = t ΔP 50 .Model state variables as a function of time, computed for eruptions showing no or a very limited acceleration in the Table S1 in Supporting Information S1 for model parameters.(a)Cumulativenumber of earthquakes N(t).(b)Red line: horizontal displacement recorded at the summit SNEG GPS station; blue line: modeled displacement u(t).(c)Continuityψ(t).The green dashed vertical line shows the t ψ 25 time.(d)Normalizedcharacteristictimeτ(t)τ0 .(e)NormalizedoverpressureΔP(t)P.The green dashed vertical line shows the t earthquake number N(t) ( p′ < 0.35, weak crack interaction; see TableS1in Supporting Information S1 for eruption parameters and caption of Figure10for more explanations).Continuity is a concave up, decreasing curve down to 0.25; it does not drop strongly below.The characteristic time τ(t) is always larger than τ 0 .The normalized power Π(t) increases before the eruption, but exceptionally reaches its second maximum.