# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2023JB027595

### Special Section:

Advances in understanding volcanic processes

### Key Points:

* Rock mass continuity decreases exponentially with earthquake number and controls magma pressure and flow

* Pre-eruptive earthquake time series are controlled by heterogeneity distribution and crack interaction for quasi-static deformation

* When crack interact, magma-edifice interaction power exhibits a maximum before the eruption, used to estimate the eruption time at PdF

### Supporting Information:

Supporting Information may be found in the online version of this article.

### Correspondence to:

J.-L. Got,
Jean-Luc.Got@univ-smb.fr

### Citation:

Got, J.-L., Peltier, A., Marsan, D., Ferrazzini, V., Brothelande, E., & Carrier, A. (2024). Pre-eruptive damage, weakening and magma-edifice coupling at Piton de la Fournaise volcano. *Journal of Geophysical Research: Solid Earth*, 129, e2023JB027595. https://doi.org/10.1029/2023JB027595

Received 10 AUG 2023
Accepted 19 FEB 2024

### Author Contributions:

**Conceptualization:** Jean-Luc Got
**Data curation:** Aline Peltier, Valerie Ferrazzini, Elodie Brothelande

**Formal analysis:** Jean-Luc Got

**Funding acquisition:** Jean-Luc Got

**Investigation:** Jean-Luc Got

**Methodology:** Jean-Luc Got

**Project administration:** Jean-Luc Got

**Resources:** Aline Peltier, Valerie Ferrazzini, Elodie Brothelande

**Software:** Jean-Luc Got, Aurore Carrier

**Supervision:** Jean-Luc Got

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_1.png" /></figure>

## Pre-Eruptive Damage, Weakening and Magma-Edifice Coupling at Piton de la Fournaise Volcano

Jean-Luc Got<fnref n="1" /> ID, Aline Peltier<fnref n="2" /> ID, David Marsan<fnref n="1" />, Valerie Ferrazzini<fnref n="2" /> ID, Elodie Brothelande<fnref n="3" />, and Aurore Carrier<fnref n="1" />

**Abstract** 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.

**Plain Language Summary** Eruptions in basaltic volcanoes are often preceded by increasing seismicity and surface deformation. In this article we study the processes at work in the volcanic edifice rock mass. A rock mass, deformed and fractured, weakens: this process is called damage. Rock elastic characteristics decrease with damage. The way damage evolves under the effect of magma pressure contributes to determine the time evolution of the pre-eruptive deformation. We show that continuity, a measure of damage, is a decreasing exponential of the cumulative number of earthquakes; this relationship is comparable to the definition of Boltzmann entropy. A model of initially elastic edifice, damaged under the pressure of a viscous, incompressible, magma in a reservoir has been developed. It is shown that, when fractures interact, continuity and magma pressure strongly drops, magma flux increases and the power of the magma-edifice interaction reaches a maximum before the eruption. A relation between the time of this maximum and the eruption time has been established and checked using 24 summit eruptions at Piton de la Fournaise; it provided an estimation within 10% of the eruption time in 60%–75% of the cases. Checking this relation with very numerous eruptions on basaltic volcanoes will allow knowing its reliability.

## 1. 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

<footnote n="1">¹Univ. Grenoble Alpes, Univ. Savoie Mont-Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, Grenoble, France,</footnote>

<footnote n="2">²Observatoire Volcanologique du Piton de la Fournaise, IPGP, Sorbonne Paris Cité, UMR, CNRS, Paris, France,</footnote>

<footnote n="3">³Department of Marine Geosciences, Rosenstiel School of Marine and Atmospheric Science, University of Miami, Miami, FL, USA</footnote>
Validation: Jean-Luc Got

Visualization: Jean-Luc Got

Writing – original draft: Jean-Luc Got

Writing – review & editing: Jean-Luc Got, Aline Peltier, David Marsan

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. Pre-eruptive 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 to 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.
<figure><img src="image_2.png" /><figcaption>**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).</figcaption></figure>

## 2. 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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normalized number of earthquakes as a function of the normalized time, the unit time being the duration of the pre-eruptive 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$.</figcaption></figure>

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.

# 3. Model

## 3.1. 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
<figure><img src="image_4.png" /><figcaption>**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 from Got et al., 2017).</figcaption></figure>

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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** 2006–2018 pre-eruptive period at Sierra Negra volcano (modified from Bell et al., 2021). The GPS station is located at the center of the caldera.</figcaption></figure>

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 $\tau_0$ (controlled by the magma viscosity and feeding system geometry), and the proportion $\delta$ 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
<figure><img src="image_6.png" /><figcaption>**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 $\mu$. In the text, the feeding system will be represented by its characteristic time $\tau_0 = \frac{8\mu H_c a_c^2}{G_0 a_c^2}$, where $G_0 = G(t = 0)$. Stars represent earthquakes.</figcaption></figure>

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 gas-permeable 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):

$$ \frac{\pi a_c^4}{8\mu H_c} (P - \Delta P(t)) = \pi a_r^3 \frac{d}{dt} \left( \frac{\Delta P(t)}{G(t)} \right) \quad (1) $$

$\Delta 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):

$$ \frac{d\Delta P(t)}{dt} + \frac{1}{\tau(t)}\Delta P(t) = \frac{P}{\tau_0}\psi(t) \quad (2) $$

where

$$ \frac{1}{\tau(t)} = \frac{1}{\tau_0}\psi(t) - \frac{d}{dt}\ln \psi(t) \quad (3) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** 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 stress-strain curves of the edifice and the magma are almost superimposed. Figure S6 in Supporting Information S1 shows the relation with the various deformation phases.</figcaption></figure>

(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.

The characteristic time $\tau(t)$ results from both the feeding (characteristic time $\tau_0$) and the damage process, represented by the Kachanov's continuity $\psi(t) = S(t)/S_0 = E(t)/E_0 = G(t)/G_0 = 1 - D(t)$, 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)$, $E_0 = E(t=0)$, $G_0 = G(t=0)$, $D(t=0) = 0$ and $\psi(t=0) = 1$. Magma reservoir overpressure, surface displacement $u(t)$ and strain $\epsilon(t)$ are related by (Got et al., 2017):

$$ \frac{\Delta P(t)}{P} = \psi(t) \frac{u(t)}{u_{el}} = \psi(t) \frac{\epsilon(t)}{\epsilon_{el}} \quad (4) $$

(where $u_{el}$ and $\epsilon_{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

$$ \tau_0 \frac{d\epsilon(t)}{dt} + \psi(t)\epsilon(t) = \epsilon_{el} \quad (5) $$

which is the strain form of the equation governing a damageable Kelvin model (Got et al., 2017), with $\tau_0 = \eta/E_0$, where $\eta$ is the viscosity in the Kelvin model.

## 3.2. A Solution for the Time Evolution of Normalized Overpressure, Surface Displacement, Strain, Magma Flow and Power Only Depending on $\tau_0$ and $\psi(t)$

The solution of the homogeneous equation associated to Equation 2 is

$$ \Delta P(t) = C e^{-\int \frac{dt}{\tau_0}} \quad (6) $$

that is,

$$ \Delta P(t) = C \psi(t) e^{-\frac{1}{\tau_0} \int \psi(t) dt} \quad (7) $$

Determining the constant C by the Euler-Lagrange method, we find

$$ C(t) = \frac{P}{\tau_0} \int_0^t e^{\frac{1}{\tau_0} \int \psi(t) dt} dt \quad (8) $$

If we write

$$ I_1(t) = \frac{1}{\tau_0} \int \psi(t) dt \quad (9) $$

the magma reservoir overpressure writes

$$ \frac{\Delta P(t)}{P} = \frac{\psi(t)}{\tau_0} e^{-I_1(t)} \int_0^t e^{I_1(t)} dt \quad (10) $$

The displacement $u(t)$ and strain $\epsilon(t)$ are

$$ \frac{u(t)}{u_{el}} = \frac{1}{\tau_0} e^{-I_1(t)} \int_0^t e^{I_1(t)} dt = \frac{\epsilon(t)}{\epsilon_{el}} \quad (11) $$

Notice that Equation 11 may be rewritten, in a compact form:
$$ \frac{u(t)}{u_{el}} = \frac{\epsilon(t)}{\epsilon_{el}} = \left( \tau_0 \frac{d}{dt} \ln \left( \int_0^t e^{I_1(t)} dt \right) \right)^{-1} \qquad (12) $$

As quoted by Got et al. (2017), magma flow may be derived from Poiseuille's law under the form:

$$ Q(t) = Q_0 \left( 1 - \frac{\Delta P(t)}{P} \right) \qquad (13) $$

so that we can define (Got et al., 2017) a normalized magma-edifice interaction power:

$$ \frac{\Pi(t)}{\Pi_0} = \frac{\Delta P(t)}{P} \frac{Q(t)}{Q_0} = \frac{\Delta P(t)}{P} \left( 1 - \frac{\Delta P(t)}{P} \right) \qquad (14) $$

where $\Pi_0 = PQ_0$.

Equations 9–12 constitute a general analytical solution of Equation 2. This solution allows expressing the system state variables ($\frac{\Delta P}{P}$, $\frac{Q}{Q_0}$, $\frac{\Pi}{\Pi_0}$) as a function of the characteristic time $\tau_0$ of the feeding system and the continuity $\psi(t)$ only. The linear elastic displacement $u_{el}$ is a scaling parameter. Notice that the characteristic time $\tau_0$ controls the viscous component of the model and $\psi$ the brittle or plastic component. In the following section, we will see that $\psi(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 $\tau_0$, $\delta$ 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).

### 3.3. Expression of the Continuity $\psi(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 - \delta)S(t)$ where $\delta = \frac{dS}{S}$ is the relative incremental decrease of the intact area $S$, for each earthquake; $\delta$ 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 - \delta)^{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:

$$ \psi(t) = \frac{S(t)}{S_0} = (1 - \delta)^{N(t)} \approx e^{-N(t)\delta} \qquad (15) $$

if $\delta \ll 1$; $N(t)\delta$ is the proportion of rupture area created at the time $t$ by $N(t)$ earthquakes. Since $\delta$ is computed by fitting the displacement data by Equation 9 knowing $N(t)$, when the completeness magnitude is larger, $N(t)$ is smaller but $\delta$ is larger. It is always possible to find a value of $\delta$ such that $N(t)\delta$ (and the continuity $\psi(t) = e^{-N(t)\delta}$) 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):

$$ \frac{d\psi(t)}{dt} = -n(t)\psi(t)\delta \qquad (16) $$
where $n(t) = \frac{dN(t)}{dt}$ is the earthquake rate; this evolution equation will be compared to the damage evolution equations of Coleman (1956, 1958) and Kachanov (1958, 1986).

### 3.4. 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):

$$n(t) = \frac{n_0}{\left(1 - \frac{t}{t_c}\right)^{p'}} \quad (17)$$

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 $\psi(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 $Q(t) = Q_0(1 - \frac{\Delta P(t)}{P})$). This increase in volume is due to damage, which induces a decrease in shear modulus. The characteristic time given by $\frac{1}{\tau(t)} = \frac{1}{\tau_0}\psi(t) - \frac{d}{dt}\ln \psi(t) = \frac{1}{\tau_0}e^{-N(t)\delta} + n(t)\delta$ (Equations 3 and 15) shows the competition, in the pre-eruptive dynamics, between the magma feeding process (first term of $\frac{1}{\tau(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 $\tau(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.

### 4. 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 pre-eruptive 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.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Model state variables as a function of the normalized time $t_c$, computed for $n_0 = 5$ earthquakes/day, $\delta = 0.0001$, $\tau_0 = 0.05 t_c$, for various values of the Omori-Utsu parameter $p'$, from 0 to 1.5 (colors). (a) Cumulative number of earthquakes $N(t)$. (b) Displacement $u(t)$. (c) Continuity $\psi(t)$. (d) Normalized characteristic time $\tau_c/\tau_0$. (e) Normalized overpressure $\Delta P/P$. (f) Normalized power $\Pi_0 p' = 0$ corresponds to $n(t) = n_0 = \text{constant}$, that is, $\psi(t) = e^{-n_0 \delta t}$, for which $n(t)$ and $\psi(t)$ have no characteristic time. $p' = 1$ corresponds to the inverse Omori law; it allows retrieving Got et al. (2017)'s solution.</figcaption></figure>

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).

## 4.1. Understanding the Relation Between Earthquake Number and Continuity Function

Continuity $\psi(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 $\tau_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 $\sigma_m$; its cumulative distribution is:

$$F(\sigma) = 1 - e^{-\int n_m(\sigma)dv} \quad (18)$$

where $n_m(\sigma) = (\sigma/\sigma_m)^m$, and $\sigma$ 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(\sigma)$ is sometimes called the stretched exponential cumulative distribution in materials science and can be compared to the continuity (Equation 15).

Coleman (1956, 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

$$\frac{dn_f(t)}{dt} = -\nu(\sigma)n_f(t) \quad (19)$$

where $n_f(t)$ is the number of unbroken fibers at time $t$, and $\nu(\sigma)$ the hazard rate, a function of the fiber stress $\sigma(t)$. The quantity $\frac{n_f(t)}{n_{f0}}$ 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(\sigma)$ deduced from Weibull's Equation 18 if we consider that $n_m(\sigma)$ represents $\nu(\sigma)$. 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 (1958, 1986) may be written

$$\frac{d\psi(t)}{dt} = -A_K \left(\frac{\epsilon_{el}}{\psi(t)}\right)^{n_K} \quad (20)$$

(where $\epsilon_{el}$ is the linear elastic strain, $A_K$ and $n_K$ are positive constants), from which the Voight (1988) evolution equation

$$\frac{d^2\Omega}{dt^2} = A_V \left(\frac{d\Omega}{dt}\right)^\beta \quad (21)$$

can be derived by using the change of variable $\psi = \frac{dt}{d\Omega}$, with $A_V = A_K\epsilon_{el}^{n_K}$ and $\beta = n_K + 2$, positive constants. Equation 21 is the basis of the Failure Forecast Method (FFM), used for eruption prediction. Equation 20 leads to

$$\psi(t) = \left(1 - \frac{t}{t_c}\right)^{\frac{1}{1+n_K}} \quad (22)$$

where $t_c = (A_K(n_K + 1)\epsilon_{el}^{n_K})^{-1}$ is the rupture time. This solution for the continuity corresponds to Equation 15 $\psi(t) = e^{-N(t)\delta}$ when $N(t)$ follows an inverse Omori law $N(t) = -n_0 t_c \ln(1 - \frac{t}{t_c})$, when $n_0 t_c \delta = \frac{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) also found a result similar to Equation 22 for crack length growth.
As a consequence, Weibull (1939)'s and Coleman (1956)'s FBM approaches, well established for tensile stress states, converge with Equation 15; Kachanov (1958)'s continuity (represented by Voight 1988's FFM in volca- nology) is a specific case of Weibull and Coleman's approaches, and of Equation 15 when N(t) is the inverse Omori law. In the following sections we will show in which mechanical conditions this latter law appears.

## 4.2. 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, 1985, 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

$$ \frac{E'}{E} = \frac{1}{1 + N \frac{\bar{v}}{v}} = \psi_{Walsh} \qquad (23) $$

where E is the Young's modulus and v is the volume of the solid, $\bar{v}$ is the volume of a spherical crack, and N the number of cracks. $N \frac{\bar{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 \frac{\bar{v}}{v} \ll 1$. Under this condition, one can write:

$$ \frac{E'}{E} = 1 - N \frac{\bar{v}}{v} = \psi_{wi} \qquad (24) $$

where $\psi_{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 \frac{\bar{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 \frac{\bar{v}}{v}$:

$$ \frac{E'}{E} \approx e^{-N \frac{\bar{v}}{v}} = \psi_{si} \qquad (25) $$

where $\psi_{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 25 is 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 (1956, 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).

## 4.3. 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:

$$ F = \sigma_0 S_0 = \sigma S \qquad (26) $$

where $\sigma_0$ is the stress applied on the initially intact, finite area $S_0$, and $\sigma$ 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 = \psi S_0$ and the Kachanov effective stress
$$ \sigma = \frac{\sigma_0}{\psi} \qquad (27) $$

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 $\sigma$ at rupture, that can be written (see, e.g., Anderson & Grew, 1977, and Heap et al., 2009 more recently, for a review):

$$ \frac{t_f}{t_{f0}} = \left(\frac{\sigma}{\sigma_0}\right)^{-\alpha} \qquad (28) $$

where $t_{f0}$ is the characteristic time of the failure occurring at stress $\sigma_0$, and $\alpha$ 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 $\frac{dt}{dN}$, which is the inverse of the earthquake rate or probability of occurrence $n(t)$. In that case, Equation 28 can be written:

$$ \frac{n(t)}{n_0} = \left(\frac{\sigma}{\sigma_0}\right)^{\alpha} \qquad (29) $$

Considering that the stress inducing earthquakes is the effective stress applied on undamaged areas (Equation 27), Equation 29 can be written:

$$ \frac{n(t)}{n_0} = \left(\frac{1}{\psi}\right)^{\alpha} \qquad (30) $$

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:

$$ \left(\frac{1}{n_0} \frac{dN(t)}{dt}\right)^{-\frac{1}{\alpha}} = \psi(N(t)) \qquad (31) $$

or

$$ \frac{dN(t)}{dt} = n_0(\psi(N(t)))^{-\alpha} \qquad (32) $$

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.

### 4.3.1. 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 32 may be written:

$$ \frac{dN(t)}{dt} = n_0 \left(1 - \frac{N(t)}{N_c}\right)^{-\alpha} \qquad (33) $$

which is a first-order linear differential equation, the solution of which is:
$$N(t) = N_c \left( 1 - \left( 1 - \frac{t}{t_c} \right)^{\frac{1}{a+1}} \right) \qquad (34)$$

if $N(t = 0) = 0$, $t_c = \frac{N_c}{n_0(a+1)}$ and $N_c = N(t = t_c)$. The corresponding earthquake rate $n(t) = \frac{dN(t)}{dt}$ is

$$n(t) = n_0 \left( 1 - \frac{t}{t_c} \right)^{-\frac{a}{a+1}} \qquad (35)$$

that is, an inverse Omori-Utsu law with a positive parameter $p' = \frac{a}{a+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 and Kaiser (1998) or in Figures 4.10–4.11 of Tang and Hudson (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 34 and $p' = \frac{a}{a+1} < 1$, we can write that $N_c = \frac{n_0 t_c}{1-p'}$ and

$$N(t) = \frac{n_0 t_c}{1 - p'} \left( 1 - \left( 1 - \frac{t}{t_c} \right)^{1 - p'} \right) \qquad (36)$$

that may be written

$$N(t) = -n_0 t_c \ln p' \left( 1 - \frac{t}{t_c} \right) \qquad (37)$$

where $\ln p'$ is the Tsallis q-logarithm with $q = p'$ (Tsallis, 1988), or

$$\frac{N(t)}{n_0 t_c} = -\ln p' \psi_0(t) \qquad (38)$$

where $n_0 t_c$ is the event number occurring without crack interaction, and $\psi_0(t) = 1 - \frac{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 = \frac{1}{\psi_0}$, all states being equally probable:

$$\frac{N(t)}{n_0 t_c} = S_{p'}(t) \qquad (39)$$

### 4.3.2. Strong Crack Interaction Case: The Inverse Omori-Utsu Law With High $p' (p' \ge 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) or in Figures 4.10–4.11 of Tang and Hudson (2011), to the stiff layer case in Gudmundsson (2005), to the dense-sand case in geomechanics (see, e.g.,
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):

$$ \frac{dN(t)}{dt} = n_0 e^{\alpha \delta N(t)} \qquad (40) $$

which is a first-order linear differential equation, the solution of which is:

$$ N(t) = -n_0 t_c \ln\left(1 - \frac{t}{t_c}\right) \qquad (41) $$

if $N(t=0) = 0$ and $t_c = \frac{1}{n_0 a \delta}$. This corresponds to a simple ($p' = 1$) inverse Omori law $n(t) = \frac{n_0}{1 - \frac{t}{t_c}}$, and may be written

$$ \frac{N(t)}{n_0 t_c} = -\ln \psi_0(t) \qquad (42) $$

As for $p' \neq 1$, the logarithm function accounts for the interactions; Equation 42 corresponds to Equation 38 when $p' = 1$. Notice that $\frac{1}{\alpha} = n_0 \delta t_c = a$, which means that in the inverse Omori law case, the stress field variation, controlled by $\alpha$, is exactly adapted to the distribution of defects, controlled by $n_0 \delta t_c$. Equation 42 may be written:

$$ N(t)\delta = -\ln \psi_{Omori}(t) \qquad (43) $$

where $\psi_{Omori}(t) = \left(1 - \frac{t}{t_c}\right)^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 and a 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, $\alpha$ is a positive finite constant so that $0 < p' < 1$, and the earthquake number $N(t)$ in Equation 34 reaches a finite value $N_c$. The limit case $p' \to 1$ implies a large $\alpha$, a very small $t_f$ (Equation 28) that is a large $n(t)$, and a small $\delta$ to have a physical $t_c = \frac{1}{n_0 \alpha \delta}$; 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 $\psi(t) = e^{-N(t)\delta}$ (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 15 when $N(t)$ follows an inverse Omori law.

## 4.4. 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.

## 4.5. 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):

$$p' = 1.15 - 1.15 \frac{n_0 t_c}{N_c} \qquad (44)$$

if $p' < 1$, and

$$p' = 1.75 - 5.7 \frac{n_0 t_c}{N_c} \qquad (45)$$

if $p' > 1$. The two lines intersect at $p' \approx 1$, for $\frac{n_0 t_c}{N_c} \approx 0.1275$.

From Equation 36 we deduce that, when $p' < 1$, $N_c = N(t = t_c) = \frac{n_0 t_c}{1-p'}$ so that $p' = 1 - \frac{n_0 t_c}{N_c} = 1 - \frac{a}{N_c \delta}$ where $a = n_0 t_c \delta$; in the theoretical case where $N_c \delta = 1$ (no overlapping of earthquake sources), $p' = 1 - a$; $p' N_c \delta$ represents the fraction of the rupture area created with crack interaction, and $a = n_0 t_c \delta$ 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 earthquake 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.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Inverse Omori-Utsu exponent $p'$ as a function of $n_0 t_c / N$. 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.</figcaption></figure>

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\delta$ and $a = n_0 t_c \delta$. 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\delta$ 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\delta$ should be systematically computed.

## 4.6. 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 $\tau_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 (Figures 10–14, Figures S12–S13 in Supporting Information S1). The model parameters $\tau_0$, $u_{el}$ and $\delta$ 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 (Figures 10–14 and Figures S12–S13 in Supporting Information S1) may be represented using two end-members (Sections 4.6.1 and 4.6.2).
<figure><img src="image_10.png" /><figcaption>Figure 9. (a) Earthquake number during secondary creep (no crack interaction) $n_0 t_c$ as a function of $\delta$. 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 = \frac{C}{\delta}$ (Monkman-Grant relation), for $C = 0.9$. (b) Parameter $a = n_0 t_c \delta$ 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.</figcaption></figure>

### 4.6.1. 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 $\Pi(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 $\tau(t) = (\frac{\psi(t)}{\tau_0} + n(t)\delta)^{-1}$, from Equation 3, drops below $\tau_0$ for 50% of these cases, which means that rupture area production $n(t)\delta$ dominates over changes in the visco-elastic response of the edifice at this stage. After the second maximum of $\Pi(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 $\Pi(t)$ and $\tau(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 $\frac{t_c - t_{GI}}{t_c}$ decreases with $p'$ for $p' > 0.6$, which can be described by a linear law (Figure 13):

$$ \frac{t_c - t_{GI}}{t_c} = -0.25p' + 0.43 \qquad (46) $$

so that an estimator of $t_c$ could be:

$$ \hat{t}_c = \frac{t_{GI}}{0.25p' + 0.57} \qquad (47) $$
<figure><img src="image_11.png" /><figcaption>**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)/τ₀. (e) Normalized overpressure ΔP(t)/P. (f) Normalized power Π/Π₀. The continuity is computed directly from *N*(*t*) (Equation 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 τ₀.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** 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 Table 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 τ₀. 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.</figcaption></figure>

Estimations provided by $\hat{t}_c$ 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 $\psi(t = t_{\psi25}) \approx 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 $\frac{t_c - t_{\psi25}}{t_c}$ (Figure 15).
<figure><img src="image_13.png" /><figcaption>**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*<sub>ψ25</sub> time. (d) Normalized characteristic time *τ*(*t*)/τ<sub>0</sub>. (e) Normalized overpressure Δ*P*(*t*)/P. The green dashed vertical line shows the *t*<sub>ΔP50</sub> time. (f) Normalized power Π/Π<sub>0</sub>. The green dashed vertical line shows the *t*<sub>GI</sub> time. The definition of Π/Π<sub>0</sub> (Equation 14) shows that *t*<sub>GI</sub> = *t*<sub>ΔP50</sub>.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13. (a) Normalized time delay $\frac{t_c-t_{G1}}{t_c}$ as a function of $p'$ (blue crosses), and model (red line, Equation 46), for $p' > 0.6$. No second maximum of $\Pi(t)$ was found for values of $p'$ lower than 0.6. (b) Sample probability of the relative residual of the estimation of $t_c$ from $t_{G1}$ (Equation 47).</figcaption></figure>

### 4.6.2. The Case Where $N(t)$ Is Quasi-Linear

In this case (Figure 14) $p' < 0.5$, and the progressive damage and crack growth occurs with weak crack interaction, stress transfer and free-surface effect. This process is not or weakly controlled by a characteristic length and the continuity is a decreasing exponential of the time, with no critical time. The overpressure drops slowly with time. This process recalls the stable mode of Lyakhovsky and Shalev (2021), when pressure and seismicity propagate by diffusing slowly in a fractured volume around an injection well. Figure 3 of Shalev and Lyakhovsky (2013) also shows how the progressive failure of the rock mass around a wellbore (due to repeated hydraulic tests) decreases the peak pressure and the post-peak pressure weakening. Characteristic time $\tau(t)$ and magma-edifice interaction power $\Pi(t)$ do not exhibit drops or maxima before the final rupture and eruption. However $\Pi(t)$ shows systematic increases before the eruption that gives information on the pre-eruptive process dynamics.

For linear $N(t)$, final rupture occurs for $\psi(t = t_{\psi25}) \approx 0.25$. Together with the values of $\frac{t_c-t_{\psi25}}{t_c}$ obtained for $p' > 0.5$, this leads to propose a Gaussian model for the normalized time delay $\frac{t_c-t_{\psi25}}{t_c}$ (Figure 15):

$$ \frac{t_c - t_{\psi25}}{t_c} = \Delta t_{\max} e^{-\left(\frac{p' - p_c}{\sigma_{p'}}\right)^2} \quad (48) $$

so that another estimator of $t_c$ could be:

$$ \hat{t}_c = \frac{t_{\psi25}}{1 - \Delta t_{\max} e^{-\left(\frac{p' - p_c}{\sigma_{p'}}\right)^2}} \quad (49) $$

with $\Delta t_{\max} \approx 0.33$, $p'_{\min} \approx 0.9$ and $\sigma_{p'} \approx 0.5$, which provides an approximate value of $t_c$ knowing $t_{\psi25}$. 60% of the values provided by this estimator, from the complete time series, are within 10% of the true $t_c$ value (Figure 15).

This symmetric model (Equation 48) shows that the very strong final acceleration and crack interaction that occur for $p' > 1.3$ is not strongly different from the weak interaction that occurs for $p' < 0.5$, from the point of view of the rupture time. This shows that a short, strong, final seismic crisis plays a minor role in the complete rupture process. This $\psi$-based estimator (Equation 49) does not require a specific interval for $p'$ values, unlike the $t_{GR}$-based model (Equation 47); it can be considered as more general. Other estimators will be discussed in further works.

Two other marginal, but however interesting, cases exist.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Model state variables as a function of time, computed for eruptions showing no or a very limited acceleration in the earthquake number *N*(*t*) (*p*′ < 0.35, weak crack interaction; see Table S1 in Supporting Information S1 for eruption parameters and caption of Figure 10 for 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 τ₀. The normalized power Π(*t*) increases before the eruption, but exceptionally reaches its second maximum.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) Normalized time delay $\frac{t_c-t_{\psi 25}}{t_c}$ between the eruption time $t_c$ and the time $t_{\psi 25}$ for which the continuity $\psi(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_{\psi 25}$ (Equation 49).</figcaption></figure>

### 4.6.3. 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).

### 4.6.4. 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.

## 5. 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
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 $\Pi(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' \ge 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 $\Pi(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_{\psi_{25}}$ or $t_{GI}$).

The Monkman-Grant parameter $a = n_0\delta 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."

## Data Availability Statement

Data used for this study were acquired by the Observatoire Volcanologique du Piton de la Fournaise—Institut de Physique du Globe de Paris (Institut de Physique du Globe de Paris, 2021). Data from the permanent seismic and Global Navigation Satellite System (GNSS) monitoring network can be downloaded from the IPGP Data Center (http://volobsis.ipgp.fr/).

* Data used in this study can be downloaded at https://doi.org/10.5281/zenodo.8207390 (Got et al., 2023).

* Software used in this study can be downloaded at https://doi.org/10.5281/zenodo.8207390 (Got et al., 2023).

## Acknowledgments

The authors would like to thank the Editor, Associate Editor and Reviewers for their constructive comments, which helped to improve the article.

## References

Amitrano, D., Grasso, J.-R., & Hantz, D. (1999). From diffuse to localised damage through elastic interaction. *Geophysical Research Letters*, 26(14), 2109–2112. https://doi.org/10.1029/1999gl900388

Amitrano, D., & Helmstetter, A. (2006). Brittle creep, damage, and time to failure in rocks. *Journal of Geophysical Research*, 111(B11). https://doi.org/10.1029/2005JB004252
Anderson, O., & Grew, P. (1977). Stress corrosion theory of crack propagation with applications to geophysics. *Reviews of Geophysics and Space Physics, 15*(1), 77–104. https://doi.org/10.1029/rg015i001p00077

Atkinson, B., & Meredith, P. (1981). Stress corrosion cracking of quartz: A note on the influence of chemical environment. *Tectonophysics*, 77(1-2), T1-T11. https://doi.org/10.1016/0040-1951(81)90157-8

Basaran, C., & Nie, S. (2004). An irreversible thermodynamics theory for damage mechanics of solids. *International Journal of Damage Mechanics*, 13(3), 205–223. https://doi.org/10.1177/1056789504041058

Basaran, C., Nie, S., Gomez, J., Gunel, E., Li, S., Lin, M., et al. (2014). Thermodynamic theory for damage evolution in solids. *Handbook of Damage Mechanics*, 13(3), 205–223.

Batdorf, S., & Crose, J. (1974). Statistical theory for the fracture of brittle structures subjected to nonuniform polyaxial stresses. *Journal of Applied Mechanics*, 41(2), 459–464. https://doi.org/10.1115/1.3423310

Bell, A., Greenhough, J., Heap, M., & Main, I. (2011). Challenges for forecasting based on accelerating rates of earthquakes at volcanoes and laboratory analogues. *Geophysical Journal International*, 185(2), 718–723. https://doi.org/10.1111/j.1365-246x.2011.04982.x

Bell, A., Hernandez, S., La Femina, P., & Ruiz, M. (2021). Uplift and seismicity driven by magmatic inflation at Sierra Negra volcano, Galapagos Islands. *Journal of Geophysical Research*, 126(7). https://doi.org/10.1029/2021JB022244

Bell, A., Naylor, M., Heap, M., & Main, I. (2011). Forecasting volcanic eruptions and other material failure phenomena: An evaluation of the failure forecast method. *Geophysical Research Letters*, 38(15), L15304. https://doi.org/10.1029/2011GL048155

Benson, P., Austria, D., Gehne, S., Butcher, E., Harnett, C., Fazio, M., et al. (2020). Laboratory simulations of fluid-induced seismicity, hydraulic fracture, and fluid flow. *Geomechanics for Energy and the Environment*, 24, 100169. https://doi.org/10.1016/j.gete.2019.100169

Bonn, D., Morton, M., Berthier, L., Divoux, T., & Manneville, S. (2017). Yield stress materials in soft condensed matter. *Reviews of Modern Physics*, 89, 35005.

Boué, A., Lesage, P., Cortés, G., Valette, B., Reyes-Davila, G., Arambula-Mendoza, R., & Budi-Santoso, A. (2016). Performance of the 'material Failure Forecast Method' in real-time situations: A Bayesian approach applied to effusive and explosive eruptions. *Journal of Volcanology and Geothermal Research*, 327, 622–633. https://doi.org/10.1016/j.jvolgeores.2016.10.002

Brantut, N., Heap, M., Meredith, P., & Baud, P. (2013). Time-dependent cracking and brittle creep in crustal rocks: A review. *Journal of Structural Geology*, 52, 17–43. https://doi.org/10.1016/j.jsg.2013.03.007

Bruner, W. (1976). Comment on “Seismic velocities in dry and saturated cracked solids” by O’Connell & Budiansky. *Journal of Geophysical Research*, 81(14), 2573–2576. https://doi.org/10.1029/jb081i014p02573

Bryant, M., Khonsari, M., & Ling, F. (2008). On the thermodynamics of degradation. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 464(2096), 2001–2014. https://doi.org/10.1098/rspa.2007.0371

Budiansky, B., & O’Connell, R. (1976). Elastic moduli of a cracked solid. *International Journal of Solids and Structures*, 12(2), 81–97. https://doi.org/10.1016/0020-7683(76)90044-5

Carrier, A., Got, J.-L., Peltier, A., Ferrazzini, V., Staudacher, T., Kowalski, P., & Boissier, P. (2015). A damage model for volcanic edifices: Implications for edifice strength, magma pressure, and eruptive processes. *Journal of Geophysical Research*, 120(1), 567–583. https://doi.org/10.1002/2014jb011485

Chan, D., Subbarayan, G., & Nguyen, L. (2012). Maximum-entropy principle for modeling damage and fracture in solder. *Journal of Electronic Materials*, 41(2), 398–411. https://doi.org/10.1007/s11664-011-1804-9

Charles, R. (1958). Static fatigue of glass. *II. Journal of Applied Physics*, 29(11), 1555–1560. https://doi.org/10.1063/1.1722992

Chastin, S., & Main, I. (2003). Statistical analysis of daily seismic event rate as a precursor to volcanic eruptions. *Geophysical Research Letters*, 30(13), 1671–1674. https://doi.org/10.1029/2003gl016900

Coleman, B. (1956). Time dependence of mechanical breakdown phenomena. *Journal of Applied Physics*, 27(8), 862–866. https://doi.org/10.1063/1.1722504

Coleman, B. (1958). Statistics and time dependence of mechanical breakdown in fibers. *Journal of Applied Physics*, 29(6), 968–983. https://doi.org/10.1063/1.1723343

Collombet, M., Grasso, J.-R., & Ferrazzini, V. (2003). Seismicity rate before eruptions on Piton de la Fournaise volcano: Implications for eruption dynamics. *Geophysical Research Letters*, 30(21), 2099–3003. https://doi.org/10.1029/2003gl017494

Costin, L. (1985). Damage mechanics in the post-failure regime. *Mechanics of Materials*, 4(2), 149–160. https://doi.org/10.1016/0167-6636(85)90013-4

Costin, L. (1987). Time-dependent deformation and failure. *Fracture Mechanics of Rock*, 167–216.

Cox, S., & Meredith, P. (1993). Microcrack formation and material softening in rock measured by monitoring acoustic emissions. *International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts*, 30(1), 11–24. https://doi.org/10.1016/0148-9062(93)90172-a

Cruden, D. (1970). A theory of brittle creep in rock under uniaxial compression. *Journal of Geophysical Research*, 75(17), 3431–3442. https://doi.org/10.1029/jb075i017p03431

Cruden, D. (1974). The static fatigue of brittle rock under uniaxial compression. *International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts*, 11(2), 67–73. https://doi.org/10.1016/0148-9062(74)92650-3

Dahm, T., & Becker, T. (1998). On the elastic and viscous properties of media containing strongly interacting in-plane cracks. *Pure and Applied Geophysics*, 151, 1–16. https://doi.org/10.1007/s000240050102

Das, S., & Scholz, C. (1981). Theory of time-dependent rupture in the earth. *Journal of Geophysical Research*, 86(B7), 6039–6051. https://doi.org/10.1029/2013JB002766

Delaney, P., & McTigue, D. (1994). Volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kilauea Volcano. *Bulletin of Volcanology*, 56(6–7), 417–424. https://doi.org/10.1007/bf00302823

Dieterich, J. (1994). A constitutive law for rate of earthquake production and its application to earthquake clustering. *Journal of Geophysical Research*, 99(B2), 2601–2618. https://doi.org/10.1029/93jb02581

Duputel, Z., Ferrazzini, V., Lengliné, O., Michon, L., Fontaine, F., & Massin, F. (2021). Seismicity of La Réunion Island. *Comptes Rendus Geoscience*, 353(S1), 237–255. https://doi.org/10.5802/crgeos.77

Duputel, Z., Lengliné, O., & Ferrazzini, V. (2019). Constraining spatiotemporal characteristics of magma migration at Piton De La Fournaise volcano from pre-eruptive seismicity. *Geophysical Research Letters*, 46(1), 119–127. https://doi.org/10.1029/2018GL080895

Dzurisin, D. (2007). *Volcano deformation: Geodetic monitoring techniques*. Springer.

Elliott, H. (1958). Stress rupture in glass. *Journal of Applied Physics*, 29(2), 224–225. https://doi.org/10.1063/1.1723074

Gehne, S., Benson, P., Koor, N., Dobson, K. J., Enfield, M., & Barber, A. (2019). Seismo-mechanical response of anisotropic rocks under simulated hydraulic fracture conditions: New experimental insights. *Journal of Geophysical Research*, 124(9), 9562–9579. https://doi.org/10.1029/2019jb017342

Glen, J. W. (1952). Experiments on the deformation of ice. *Journal of Glaciology*, 2(12), 111–114.
Got, J.-L., Carrier, A., Marsan, A., Jouanne, V., Vogfjord, K., & Villemin, T. (2017). An analysis of the nonlinear magma-edifice coupling at Grímsvötn volcano (Iceland). *Journal of Geophysical Research: Solid Earth*, 122(2), 826–843. https://doi.org/10.1002/2016JB012905

Got, J.-L., Peltier, A., Marsan, D., Ferrazzini, V., Brothelande, E., & Carrier, A. (2023). Pre-eruptive damage, weakening and magma-edifice coupling at Piton de la Fournaise volcano [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.8207390

Gudmundsson, A. (2005). The effects of layering and local stresses in composite volcanoes on dyke emplacement and volcanic hazards. *C.R. Geoscience*, 337(13), 1216–1222. https://doi.org/10.1016/j.crte.2005.07.001

Gudmundsson, A. (2012). Strengths and strain energies of volcanic edifices: Implication for eruptions, collapse calderas, and landslides. *Natural Hazards and Earth System Sciences*, 12(7), 2241–2258. https://doi.org/10.5194/nhess-12-2241-2012

Gupta, P. (1982). Relation between power and exponential laws of slow crack growth. *Journal of the American Ceramic Society*, 65(10), C-163. https://doi.org/10.1111/j.1151-2916.1982.tb10350.x

Heap, M., Baud, P., Meredith, P., Bell, A., & Main, I. (2009). Time-dependent brittle creep in Darley Dale sandstone. *Journal of Geophysical Research*, 114(B07203). https://doi.org/10.1029/2008jb006212

Heap, M., Villeneuve, M., Albino, F., Farquharson, J., Brothelande, E., Amelung, F., et al. (2020). Towards more realistic values of elastic moduli for volcano modelling. *Journal of Volcanology and Geothermal Research*, 390, 106684. https://doi.org/10.1016/j.jvolgeores.2019.106684

Heimisson, E., & Segall, P. (2020). Physically consistent modeling of dike induced deformation and seismicity: Application to the 2014 Bar-darbunga dike, Iceland. *Journal of Geophysical Research*, 125(2). https://doi.org/10.1029/2019JB018141

Herring, T., King, R., Floyd, M., & McClusky, S. (2018). *Introduction to GAMIT/GLOBK*. Release 10.7. Memorandum of MIT.

Idris, R., Abdullah, S., Thamburaja, P., & Zaidi Omar, M. (2020). Prediction of fatigue crack growth rate based on entropy generation. *Entropy*, 22(1), 9. https://doi.org/10.3390/e22010009

Institut de Physique du Globe de Paris, I. (2021). *Data collection of the volcanological observatory of Piton de la Fournaise*. Institut de Physique du Globe de Paris. https://doi.org/10.18715/REUNION.OVPF

Kachanov, L. M. (1958). On creep rupture time. *Isv. Akad. Nauk. SSR. Otd Tehkh. Nauk*, 8, 26–31.

Kachanov, L. M. (1986). *Introduction to continuum damage mechanics* (pp. 1–135). Martinus Nijhoff Publishers.

Kilburn, C. (2018). Forecasting volcanic eruptions: Beyond the Failure Forecast Method. *Frontiers in Earth Science*, 6, 133. https://doi.org/10.3389/feart.2018.00133

Kilburn, C., & Voight, B. (1998). Slow rock fracture as eruption precursor at Soufrière Hills volcano, Montserrat. *Geophysical Research Letters*, 25, 19–3668. https://doi.org/10.1029/98g101609

Kozono, T. (2021). The dynamics of dual-magma-chamber system during volcanic eruptions inferred from physical modeling. *Earth Planets and Space*, 73(1), 103. https://doi.org/10.1186/s40623-021-01421-4

Krajcinovic, D. (1996). *Damage mechanics*. Elsevier.

Krajcinovic, D., & Silva, M. (1982). Statistical aspects of the continuous damage theory. *International Journal of Solids and Structures*, 18(7), 551–562. https://doi.org/10.1016/0020-7683(82)90039-7

Lajtai, E., & Schmidtke, R. (1986). Delayed failure in rock loaded in uniaxial compression. *Rock Mechanics and Rock Engineering*, 19(1), 11–25. https://doi.org/10.1007/bf01043444

Lajtai, E., Scott Duncan, E., & Carter, B. (1991). The effect of strain rate on rock strength. *Rock Mechanics and Rock Engineering*, 24(2), 99–109. https://doi.org/10.1007/bf01032501

Lengliné, O., Marsan, D., Got, J., Pinel, V., Ferrazzini, V., & Okubo, P. (2008). Seismicity and deformation induced by magma accumulation at three basaltic volcanoes. *Journal of Geophysical Research*, 113(B12), B12305. https://doi.org/10.1029/2008jb005937

Li, V. (1987). Mechanics of shear rupture applied to earthquake zones. In B. Atkinson (Ed.), *Fracture mechanics of rock* (pp. 351–428). Academic Press.

Lyakhovsky, V., & Shalev, E. (2021). Runaway versus stable fracturing during hydraulic stimulation: Insights from the damage rheology modeling. *Rock Mechanics and Rock Engineering*, 54(10), 5449–5464. https://doi.org/10.1007/s00603-021-02395-1

Maimon, O., Lyakhovsky, V., Melnik, O., & Navon, O. (2012). The propagation of a dyke driven by gas-saturated magma. *Geophysical Journal International*, 189(2), 956–966. https://doi.org/10.1111/j.1365-246x.2012.05342.x

Main, I. (1999). Applicability of time-to-failure analysis to accelerated strain before earthquakes and volcanic eruptions. *Geophysical Journal International*, 139(3), F1–F6. https://doi.org/10.1046/j.1365-246x.1999.00004.x

Main, I. (2000). A damage mechanics model for power-law creep and earthquake aftershock and foreshock sequences. *Geophysical Journal International*, 142(1), 151–161. https://doi.org/10.1046/j.1365-246x.2000.000136.x

Main, I., & Meredith, P. (1989). Classification of earthquake precursors from a fracture mechanics model. *Tectonophysics*, 167(2-4), 273–283. https://doi.org/10.1016/0040-1951(89)90078-4

McGuire, W., & Kilburn, C. (1997). Forecasting volcanic events: Some contemporary issues. *Geologische Rundschau*, 86(2), 439–445. https://doi.org/10.1007/s005310050152

Monkman, F., & Grant, N. (1956). An empirical relationship between rupture life and minimum creep rate in creep rupture tests. *Proc of the ASTM*, 56, 593–605.

Mould, R., & Southwick, R. (1959). Strength and static fatigue of abraded glass under controlled ambient conditions: II, effect of various abrasions and the universal fatigue curve. *Journal of the American Ceramic Society*, 42(12), 582–592. https://doi.org/10.1111/j.1151-2916.1959.tb13578.x

Nanjo, K. Z. (2017). A fiber-bundle model for the continuum deformation of brittle material. *International Journal of Fracture*, 204(2), 225–237. https://doi.org/10.1007/s10704-016-0175-x

Nicot, F., Wang, X., Wautier, A., Wan, R., & Darve, F. (2023). Shear banding as a dissipative structure from a thermodynamic viewpoint. *Journal of the Mechanics and Physics of Solids*, 179, 105394. https://doi.org/10.1016/j.jmps.2023.105394

Osara, J., & Bryant, M. (2019). Thermodynamics of fatigue: Degradation-entropy generation methodology for system and process characterization and failure analysis. *Entropy*, 21(7), 685. https://doi.org/10.3390/e21070685

Peltier, A., Staudacher, T., & Bachlèlery, P. (2007). Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data. *Journal of Geophysical Research*, 112(B3). https://doi.org/10.1029/2006jb004379

Reches, Z., & Lockner, D. (1994). Nucleation and growth of faults in brittle rocks. *Journal of Geophysical Research*, 99(B9), 18159–18173. https://doi.org/10.1029/94jb00115

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., & Bascou, P. (2014). A two-magma chamber model as a source of deformation at Grímsvötn Volcano, Iceland. *Journal of Geophysical Research*, 119(6), 4666–4683. https://doi.org/10.1002/2013jb010569

Rivalta, E., & Dahm, T. (2004). Dyke emplacement in fractured media: Application to the 2000 intrusion at Izu islands, Japan. *Geophysical Journal International*, 157(1), 283–292. https://doi.org/10.1111/j.1365-246x.2004.02193.x
Rivalta, E., & Dahm, T. (2006). Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface. *Geophysical Journal International*, **166**(3), 1424–1439. https://doi.org/10.1111/j.1365-246x.2006.02962.x

Rivalta, E., & Segall, P. (2008). Magma compressibility and the missing source for some dike intrusions. *Geophysical Research Letters*, **35**(4), L04306. https://doi.org/10.1029/2007g01032521

Roman, D., & Cashman, K. (2006). The origin of volcano-tectonic earthquake swarms. *Geology*, **34**(6), 457–460. https://doi.org/10.1130/G22269.1

Rundle, J., Turcotte, D., Shcherbakov, R., Klein, W., & Sammis, C. (2003). Statistical physics approach to understanding the multiscale dynamics of earthquake fault systems. *Reviews of Geophysics*, **41**(4). https://doi.org/10.1029/2003RG000135

Schmid, A., & Grasso, J. (2012). Omori law for eruption foreshocks and aftershocks. *Journal of Geophysical Research*, **117**(B07302). https://doi.org/10.1029/2011jb008975

Schmid, A., Grasso, J.-R., Clarke, D., Ferrazzini, V., Bachlèry, P., & Staudacher, T. (2012). Eruption forerunners from multiparameter monitoring and application for eruptions time predictability (Piton de la Fournaise). *Journal of Geophysical Research*, **117**(B11), B11203. https://doi.org/10.1029/2012JB009167

Segall, P., Llenos, A., Yun, S., Bradley, A., & Syracuse, E. (2013). Time-dependent dike propagation from joint inversion of seismicity and deformation data. *Journal of Geophysical Research*, **118**(11), 5785–5804. https://doi.org/10.1002/2013jb010251

Shalev, E., & Lyakhovsky, V. (2013). The processes controlling damage zone propagation induced by wellbore fluid injection. *Geophysical Journal International*, **193**(1), 209–219. https://doi.org/10.1093/gji/ggt002

Spera, F. (2000). Physical properties of magma. In H. Sigurdsson (Ed.), *Encyclopedia of volcanoes* (pp. 171–190). Academic Press.

Tang, C. (1997). Numerical simulation of progressive rock failure and associated seismicity. *International Journal of Rock Mechanics and Mining Sciences*, **34**(2), 249–261. https://doi.org/10.1016/s0148-9062(96)00039-3

Tang, C., & Hudson, J. (2011). *Rock failure mechanisms explained and illustrated*. CRC Press.

Tang, C., & Kaiser, P. (1998). Numerical simulation of cumulative damage and seismic energy release during brittle rock failure – Part I: Fundamentals. *International Journal of Rock Mechanics and Mining Sciences*, **35**(2), 113–121. https://doi.org/10.1016/s0148-9062(97)00009-0

Tsallis, C. (1988). Possible generalization of Boltzmann-Gibbs statistics. *Journal of Statistical Physics*, **52**(1–2), 479–487. https://doi.org/10.1007/bf01016429

Tucker, J., Chan, D., Subbarayan, G., & Handwerker, C. (2014). Maximum entropy fracture model and its use for predicting cyclic hysteresis in Sn3.8Ag0.7Cu and Sn3.0Ag0.5 solder alloys. *Microelectronics Reliability*, **54**(11), 2513–2522. https://doi.org/10.1016/j.microrel.2014.04.012

Turcotte, D., & Glasscoe, M. (2004). A damage model for the continuum rheology of the upper continental crust. *Tectonophysics*, **383**(1–2), 71–80. https://doi.org/10.1016/j.tecto.2004.02.011

Turcotte, D., Newman, W., & Shcherbakov, R. (2003). Micro and macroscopic models of rock fracture. *Geophysical Journal International*, **152**(3), 718–728. https://doi.org/10.1046/j.1365-246x.2003.01884.x

Vasseur, J., Wadsworth, F., Heap, M., Main, I., & Dingwell, D. (2017). Does an inter-flaw length control the accuracy of rupture forecasting in geological materials? *Earth and Planetary Science Letters*, **475**, 181–189. https://doi.org/10.1016/j.epsl.2017.07.011

Vinciguerra, S., Meredith, P., & Hazzard, J. (2004). Experimental and modeling study of fluid pressure-driven fractures in Darley Dale sandstone. *Geophysical Research Letters*, **31**(9), L09609. https://doi.org/10.1029/2004gl019638

Voight, B. (1988). A method for prediction of volcanic eruptions. *Nature*, **332**(6160), 125–130. https://doi.org/10.1038/332125a0

Walsh, J. (1965). The effect of cracks on the uniaxial elastic compression of rock. *Journal of Geophysical Research*, **70**(2), 399–411. https://doi.org/10.1029/jz070i002p00399

Weibull, W. (1939). *A statistical theory of the strength of materials* (Vol. 151, pp. 1–45). Royal Institution for Engineering and Research.

Weiss, J., Girard, L., Gimbert, F., Amitrano, D., & Vandembroucq, D. (2014). Finite size effects on compressive strength. *Proceedings of the National Academy of Sciences of United States of America*, **111**(17), 6231–6236. https://doi.org/10.1073/pnas.1403500111

Wiederhorn, S., & Bolz, L. (1970). Stress corrosion and static fatigue of glass. *Journal of the American Ceramic Society*, **50**(10), 543–548. https://doi.org/10.1111/j.1151-2916.1970.tb15962.x

Wilkins, B. (1980). Slow crack growth and delayed failure of granite. *International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts*, **17**(6), 365–369. https://doi.org/10.1016/0148-9062(80)90520-3

Zok, F. (2017). On weakest link theory and Weibull statistics. *Journal of the American Ceramic Society*, **100**(4), 1265–1268. https://doi.org/10.1111/jace.14665