# Damage in cohesive granular materials: simulations and geophysical implications

Vincent Canel, Michel Campillo, Xiaoping Jia, Ioan R. Ionescu

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

## To cite this version:

Vincent Canel, Michel Campillo, Xiaoping Jia, Ioan R. Ionescu. Damage in cohesive granular materials: simulations and geophysical implications. Comptes Rendus. Géoscience, 2023, 355 (S3), pp.1-21. 10.5802/crgeos.203. hal-03997648

HAL Id: hal-03997648

https://hal.science/hal-03997648

Submitted on 22 May 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Vincent Canel, Michel Campillo, Xiaoping Jia and Ioan R. Ionescu

**Damage in cohesive granular materials: simulations and geophysical implications**

Published online: 13 April 2023

https://doi.org/10.5802/crgeos.203

**Part of Special Issue:** New Developments in Passive Seismic Imaging and Monitoring

**Guest editors:** Michel Campillo (University Grenoble-Alpes, Institut des Sciences de la Terre, Grenoble, France),

Andrew Curtis (School of GeoSciences, University of Edinburgh, Scotland),

Anne Obermann (Swiss Seismological Service, ETH, Zurich, Switzerland) and Nikolai

Shapiro (CNRS, Institut des Sciences de la Terre, Grenoble, France)

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Géoscience — Sciences de la Planète sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1778-7025
1984, Schaff and Beroza, 2004, Langlois and Jia, 2014). Noise-based passive monitoring allows to envision long-term continuous monitoring of the seismic velocity [Brenguier et al., 2008, 2014, Sens-Schönfelder and Wegler, 2006, Rivet et al., 2014]. Examples of these changes are now well documented, and it is possible to detect them continuously with unprecedented precision that can reach a few $10^{-5}$ for relative velocity changes [Wang et al., 2017] and very short time resolution [Mao et al., 2019b, Sens-Schönfelder and Eulenfeld, 2019].

One difficulty with these measurements is to distinguish between changes related to internally induced deformations and changes related to external forcing, such as precipitation [Sens-Schönfelder and Wegler, 2006, Barajas et al., 2021], temperature [Meier et al., 2010], or tides [Sens-Schönfelder and Eulenfeld, 2019]. As the causes of external forcing are fairly well known, correction strategies have been proposed [Wang et al., 2017]. Assessing these changes in the vicinity of rupture zones will improve through the development of sensitivity kernels for spatial inversion methods [Margerin et al., 2016, Obermann et al., 2016, van Dinther et al., 2021, Barajas et al., 2021, Zhang et al., 2021].

With the striking increase achieved in quantity and quality of seismic observations and the novel methods now possible, a way is opening to revisit the vision of earthquake processes. At the same time, geodetic measures reach impressive precision. At time scales from days to years, evidence for changes in elastic moduli in response to earthquakes [Brenguier et al., 2008] and transient tectonic deformations [Rivet et al., 2014] has been shown in the form of a velocity drop followed by a slow recovery similar to laboratory observations [TenCate et al., 2000]. These changes are also observed at the time scale of days in a mining environment [Olivier et al., 2015] and in laboratory experiments of stick-slip and slow slip [Scudereri et al., 2016], where a drop in the elastic moduli starts even before the slip is observed.

Such effects have been largely studied in rocks with laboratory experiments. They have revealed that the application of a uniaxial stress involves elastic wave velocity anisotropy [Nur and Simmons, 1969] that could be linked to the elastic nonlinear behavior of rocks [Johnson et al., 1996, Johnson and Rasolofosaon, 1996, Pasqualini et al., 2007]. Direct comparisons between rock degradation (linked to

a damage quantified by the distribution of microcracks) and wave velocity measurements could be highlighted experimentally [Hamiel et al., 2009]. Furthermore, theoretical considerations based on the inversion of the measured wave velocities and resulting microcrack density tensors allowed to describe the microcracks evolution and anisotropy [Sayers and Kachanov, 1991, 1995, Schubnel and Guéguen, 2003, Stanchits et al., 2006, Hall et al., 2008]. For the following numerical simulations, we consider granular materials rather than real rocks as they have been studied as synthetic rocks in previous experiments [Langlois and Jia, 2014, Canel et al., 2020].

An important issue is to explain how the very small macro-scale deformation that is associated with tectonic deformation can induce an *observable velocity drop*. In this study, we focus on the drop effect without consideration of the slow logarithmic type relaxation effect [TenCate et al., 2000] observed in the earth after a strong drop [Brenguier et al., 2008].

Here we use a simple cemented granular media model [Dvorkin et al., 1994, Langlois and Jia, 2014, Hemmerle et al., 2016] to represent the behavior of rocks around faults that are intensely fractured and cemented by precipitation due to fluid transfer. Our interest in this simulation is to link it to the observed temporal variations of seismic wave propagation velocities. More precisely, in the following, we use numerical simulation models to analyze the relation between the wave velocity drop and the damage evolution. We focus on numerical simulations in an edometric experimental setup, such that some comparisons with real experiments can be done. The edometric experiment is a typical experiment in soils and granular materials [Evesque, 2000, Sawicki and Swidzinski, 1995, Langlois and Jia, 2014], which involve macroscopic deformation of up to 15%. In the first stage, the goal is to deform and damage a medium in which we will later study wave propagation [Langlois and Jia, 2014]. Only small overall deformation (up to 1.5%) will be considered here. However, as even for small macroscopic strain the microscopic bond deformation is large and very heterogeneous, large deformation modeling will be considered at the microscopic scale.

The numerical granular material considered in this paper is made from a dense elastic bead packing that includes some cement. The numerical simulations can highlight the microscopic hetero-
geneities. For this we used the finite element method, instead of the discrete element method [Radjai and Dubois, 2011], and we meshed both the beads and the bonds. As far as we know, this approach has not been used in granular modeling. In contrast with the discrete element method approach, which is widely used in granular physics, the finite element method approach provides a fine description of the damage mechanism at bonded contacts and of the elastic wave propagation through the bonds and the particles (beads).

To limit the computational load, we consider only two-dimensional (2D) numerical simulations. The mechanical modeling of the bonds relies on a simple elasto-plastic model with damage, while the beads are assumed to be elastic and isotropic. Following the edometric experimental setup, we first detail the adopted geometry for the cemented granular material, and then we distinguish two processes in the following sections. The first process is quasi-static loading of the heterogeneous macroscopic sample that involves damage to the bonds. This first modeling aims at creation of a realistic state of deformation at the microscopic scale. The second modeling is the dynamic propagation of a small perturbation in a linear elastic regime for a given state of damage of the material. This second process can be considered as the dynamic version of the “tangent problem” associated to the quasi-static evolution. It represents the actual situation encountered in geophysical monitoring of active faults: an elastic, weak amplitude wave propagates in a medium where damage results from a slow internal deformation. Note that the macroscopic damage of the sample, which is defined from the variation of the macroscopic elastic coefficients, cannot be found from the strain-stress curve obtained from quasistatic unidirectional loading. For this, we would need several unloading/reloading experiments. The macroscopic damage is therefore probed by the wave propagation, as in the seismological measurements of the wave speed.

## 2. Quasistatic loading

### 2.1. Macroscopic setting

Let us first describe the numerical edometric experiment. The “sample” (or cell) of the cemented granular material initially occupies the 2D domain

$$ \mathcal{D}_0^{\text{macro}} = (0, L) \times (0, l_0) \text{ (with } L \text{ the constant length and } l_0 \text{ the initial width of the sample, see Figure 1) in the plane strain configuration. On its boundary, the sample is in contact without friction with four rigid walls. Three of them are at rest while the upper one moves slowly in the vertical direction. The process is strain driven: the sample width } l \text{ and then the macroscopic displacement } \Delta l = |l - l_0| \text{ evolve linearly with time. The microscopic deformation process is heterogeneously activated by the macroscopic deformation } \varepsilon_{\text{macro}} = \Delta l / l_0, \text{ which increases slowly from 0 to the final strain } \varepsilon_{\text{macro}}^f = 1.58\%. \text{ During the loading process, to find the macroscopic stress } \sigma_{\text{macro}} = F_{\text{macro}}/L, \text{ we compute the resultant force } F_{\text{macro}} \text{ that acts on the upper plate. In this configuration, the macroscopic Hooke's law at the macroscopic level is written as } \sigma_{\text{macro}} = M_{\text{macro}}\varepsilon_{\text{macro}}, \text{ where } M_{\text{macro}} \text{ is the macroscopic uniaxial strain modulus or P-wave modulus } (M = \lambda + 2\mu = E(1-\nu)/((1+\nu)(1-2\nu))) \text{ with } \lambda \text{ and } \mu \text{ the macroscopic first and second Lamé coefficients and } E \text{ and } \nu \text{ the macroscopic Young and Poisson moduli respectively).} $$

### 2.2. Microscopic setting

The cemented granular material, which occupies at the moment $t \in [0, T^{\text{quasi-static}}]$ the domain $\mathcal{D}_t \subset \mathcal{D}_0^{\text{macro}}$ is composed of grains, denoted by $\mathcal{G}_t$, and bonds, denoted by $\mathcal{B}_t$, i.e. $\tilde{\mathcal{D}}_t = \tilde{\mathcal{G}}_t \cup \tilde{\mathcal{B}}_t$ and $\mathcal{G}_t \cap \mathcal{B}_t = \emptyset$. For the initial geometric settings of grains $\mathcal{G}_0$ we used 382 disks (to model glass beads) of radii 300 µm to 600 µm. The beads are connected by the cement $\mathcal{B}_0$ composed of 956 bonds and gathered in a cell $\mathcal{D}_0^{\text{macro}}$ of size 3 cm x 1 cm. The word “bond” stands for the deformable region which links two beads. The packing fraction (i.e., the ratio between the beads area and the total area of the sample) is initially 77.2% and increases slowly during the compression, to reach 78.4% at the final strain $\varepsilon_{\text{macro}}^f$. The ratio between the bonds area and the beads area is 6.32%, while the ratio between the bond area and the total initial area of the sample is 4.88%. Even if the bonds area is small compared to the beads area, the bonds are modeled as a set of continuum 2D bodies.

After creation of the packing, we checked that the bonds are isotropically generated, so as not to induce any geometric anisotropy from the fabric that might hide or perturb possible anisotropic effects due to the future quasi-static loading. Figure 2b shows the polar
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Representation of the quasi-static edometric numerical experiment.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Geometric statistics of the initial isotropic packing. (a) Bond length as a function of the orientation angle of the corresponding bond. (b) Polar probability distribution of the orientation angles of the bonds. (c) Probability distribution of the lengths of the bonds. (d) Definition of the orientation angle $\theta$ of a bond.</figcaption></figure>

probability distribution of the orientation angles of the bonds (i.e., angle between the (*O*x) axis and the axis that connects the center of the beads, as schematized in Figure 2d). Half the polar diagram is shown, as it is $\pi$-periodic: each bond is indeed indexed on the basis of both its associated bead, with an angle $\theta$, and the other associated bead, with an angle $\theta + \pi$. The distribution is globally isotropic (a sample with more bonds would exhibit a smoother distribution).

Figure 2c shows the distribution of the length of the bonds, defined as the distance between the centers of the beads concerned and the sum of their radii; this then represents a length along the bond

symmetrical axis. The minimal values correspond to the closest beads and tend towards zero, as there is no minimal limit for the generation of bonds in the sample considered. The majority of the bonds have a length less than 0.05 mm.

Further geometric anisotropy can come from correlations between the lengths and the orientations of the bonds. Indeed, as we expect that short bonds are more affected by the macroscopic deformation, we verified that there is no privileged direction for them. Figure 2a shows the lengths of the bonds as functions of the orientation angles of the corresponding bonds in a polar plot. The isotropic distribution shows that
<table><caption>**Table 1.** Parameters for simulations</caption><thead><tr><th></th><th>ρ (kg/m³)</th><th>E₀ (GPa)</th><th>ν₀ (GPa)</th><th>M₀ (GPa)</th><th>cₚ₀ (m/s)</th><th>c₅₀ (m/s)</th></tr></thead><tbody><tr><td>Beads</td><td>2600</td><td>60</td><td>0.24</td><td>70.72</td><td>5200</td><td>3100</td></tr><tr><td>Bonds</td><td>1200</td><td>3.5</td><td>0.33</td><td>5.18</td><td>2100</td><td>1000</td></tr></tbody></table>

The index 0 refers to the state without damage.

there is no correlation between length and orientation angle.

## 2.3. *Mechanical modeling*

For the mechanical modeling of the bond material, which is detailed in Appendix, we chose the model of Lemaitre and Chaboche [Lemaitre and Chaboche, 1994] coupling isotropic ductile plastic damage and with an elasto-plastic law in the framework of an Eulerian (large deformations) description [see for instance Belytschko et al., 2000].

For this, we considered the additive decomposition of the rate deformation tensor in the elastic and plastic rates of deformation. For the elastic range, we considered an isotropic hypo-elastic law (i.e., the large strain generalization of Hooke's law with the Lamé coefficients $\lambda$ and $\mu$ written in terms of the Jauermann rate of the Cauchy stress tensor and the elastic rates of deformation). The plastic rate of deformation is related to the Cauchy stress tensor through the flow rule associated to the classical Von-Mises yield criterion (without hardening). The weakening effect is characterized by a damage (phenomenological) parameter $d \in [0, 1]$ and its evolution law is related to the cumulated plastic strain $\epsilon_p$.

The model parameters used here have not been characterized or calibrated experimentally. This important task is beyond the scope of present numerical study. However, the mechanical parameters have been chosen to correspond to a very ductile cement, as tetradecane or eicosane, that has been used in laboratory experiments showing seismic velocity reduction during compression [Langlois and Jia, 2014, Canel et al., 2020]. The elastic coefficients of the undamaged material are given in Table 1 and the yield limit is $\kappa_0 = 135$ MPa. For simplicity, the dependence of damage rate on the cumulated plastic strain is given through a piecewise linear function. When the cumulated plastic strain reaches an activation level

$\epsilon_{p,activation} = 2\%$, we consider the beginning of damage process, while the maximal level of damage is $d_{max} = 0.8$ corresponding to $\epsilon_{p,max} = 20\%$.

The beads material was supposed to be purely elastic (see Table 1 for the Lamé coefficients corresponding to glass), with no damage effects or plastic strains. This choice could appear too simple but corresponds to a granular material made from glass beads used in laboratory experiments for which no bead crushing has been observed [see for instance Langlois and Jia, 2014, Canel et al., 2020]. Moreover, this choice allows us to focus on the bond material, and to analyze the crucial role played by the bonds damage in the wave velocity drop. This is due to damage localization in the bonds, which are pointed out in the next numerical results obtained with the choice of damage parameters presented above. The presence of a maximal level of damage in the model prevents total bonds failure (that would be very difficult to handle in an elasto-plastic FE computation of a large number of bonds) but is not expected to affect the overall results at the low levels of macroscopic deformation at which we limit our simulations.

## 2.4. *Numerical modeling*

As the bonds are submitted to large deformation (even for small macroscopic deformation), we used an arbitrary Lagrangian Eulerian (ALE) description [Huerta and Casadei, 1994, Wang and Gadala, 1997, Ghosh and Raju, 1996, Rodríguez-Ferran et al., 2002]. This description incorporates the advantages of the Lagrangian and Eulerian descriptions, and can avoid some of their drawbacks. Usually in solid mechanics, where we do not deal with large mass fluxes among different parts of the sample and the strains are not too large, the Lagrangian kinematics formulations are intensively used. However, because of severe distortion of elements in some practical problems, the determinant of the Jacobian matrix can become negative, which results in numerical errors. Eulerian
<figure><img src="image_7.png" /><figcaption>Figure 3. The initial mesh used for the quasi-static computations.</figcaption></figure>

methods coupled with an ALE description (which either uses a fixed mesh or adapts the mesh at each time step) can eliminate the problems associated with distorted meshes. With the ALE method adopted here, the computational grid can be moved arbitrarily, to optimize element shapes independently of material deformation.

For the numerical integration of the equilibrium of the hypo-elastic-plastic model described above, to find the Eulerian unknowns acting on the domain $D_{t+\Delta t}$, we used a classical implicit backward radial return mapping scheme [e.g., see Belytschko et al., 2000] with explicit treatment of the damage field on the finite element discretization of domain $D_t$. The discretization was done with [P2-continuous] finite elements for the velocity fields and [P1-discontinuous] for the other unknowns (e.g., stress, damage, plastic deformation).

Since the investigated processes are very different, the quasi-static and dynamic problems require two distinct meshes (one can see the difference by comparing Figures 3 and 8). For the quasi-static problem, as the beads are almost not deformed outside of their peripheral zone, it is not necessary to mesh them finely at the core, but only close to the bonds that need a fine mesh to correctly describe their large deformation. The initial mesh is then finer in the bonds and the periphery of the beads, with a given ratio between the size of the edges (see Figure 3). Moreover, to capture the shear bands, we used an adaptive mesh technique with respect to the plastic strain rate norm $|D^p|$. This means that regions where the plastic strain rate is larger will have a fine mesh, while elsewhere the mesh is coarse. The ratio between the sizes of the fine and coarse meshes was 1/4.

## 2.5. Results

### 2.5.1. Macroscopic level

Figure 4 shows the macroscopic normal load $\sigma_{\text{macro}}$ versus the imposed strain on the sample $\epsilon_{\text{macro}}$. The macroscopic load initially increases linearly with the displacement, which highlights an elastic behavior. The linear fit close to the origin (see Figure 4) gives an estimation of the initial (undamaged) macroscopic uniaxial strain modulus (P-wave modulus) $M_{\text{macro}}^0 \approx 12.6$ GPa, which is larger than the modulus of the bonds but much less than the bead modulus. After some point, the macroscopic strain-stress curve deviates from this linear trend and then reaches a plateau, which is a typical strain-stress curve for an elasto-plastic material without hardening.

However, the macroscopic damage $d_{\text{macro}}$ of the sample, defined as

$$d_{\text{macro}} = 1 - \frac{M_{\text{macro}}}{M_{\text{macro}}^0},$$

cannot be found from the strain-stress curve only. Indeed, to get the macroscopic uniaxial stress modulus we cannot use the secant strain-stress modulus and several unloading and reloading experiments during the main loading process are needed.

This kind of quasi-static unloading/reloading experiments with fluid pressure are possible at the laboratory scale or at a small geophysical scale corresponding to reservoirs for example (in the framework of the oil and gas industry). However, as far as we know these techniques have never been applied in Nature to a fault zone at the scale relevant
<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) Macroscopic stress $\sigma_{\text{macro}}$ (black; in × 33 MPa) and the averaged damage $d_{\text{ave}}$ (blue; in %) as functions of the strain $\varepsilon_{\text{macro}}$ (in %). A linear fit of the load (red) highlights its first elastic trend. The chosen states for the dynamic computations are indicated with vertical green lines. (b) Zoom in around the low strains.</figcaption></figure>

to earthquake studies. One can think about the tides effect but our study does not target it. Indeed, different techniques have been used to study the tide effect and actually highlighted the velocity dependence on tide deformation independently from the presence of a fault system [Sens-Schönfelder and Eulenfeld, 2019, Mao et al., 2019a, Reasenberg and Aki, 1974, Takano et al., 2014, Delorey et al., 2021]. Nevertheless, the continuous nature of tides implies that they are not a source of evolution damage, or over time scales and amplitudes that are incompatible with their measurements. In our perspective, nonlinearity and damage evolution are quite different. This is why we did not use the quasi-static unloading/reloading technique to measure the macroscopic damage but the wave speed probing (see next section), which is much more adapted to the targeted geophysical measurements.

We computed and plotted the microscopic averaged damage $d_{\text{ave}}(t)$, defined as the ratio between the integral of the damage, that vanishes outside the bonds, and the area of the sample $|\mathcal{D}_t|$ which is the area of the set of grains and bonds. Since plastic deformations preserve the area (volume in 3D) the only variation of the area $|\mathcal{D}_t|$ is due to elastic deformations. The deviation of the macroscopic load from a linear trend is clearly correlated to the accumulation of damage. Even when the load reaches a stationary value, $d_{\text{ave}}(t)$ keeps increasing (which can be

associated to anisotropic effects of the compression; see Section 4 for more details).

### 2.5.2. Microscopic level

Even for small macroscopic deformation, local plastic irreversible deformation occurs in the bonds. The first microscopic plastic and damage effects occur at $\varepsilon_{\text{macro}} = 0.096\%$ and $\varepsilon_{\text{macro}} = 0.132\%$, respectively (see Figure 4 for the damage). This is consistent with the chosen microscopic laws with successive plastic and damage thresholds (see Section 2.3 and the Appendix for more details). The damage process is very heterogeneous as it occurs along the force chains. Some of the bonds are not damaged, and some of them are completely damaged, although there are a lot of bonds in an intermediary state where the damage is localized in microscopic shear bands. We note that there is no macroscopic localization of the damage (i.e., no noticeable development of a macroscopic shear band).

Figure 5 shows the deformed configuration at the end of the edometric numerical experiment ($\varepsilon_{\text{macro}}^f = 1.58\%$), with a map of the deviatoric stress (top) that shows mainly vertical heterogeneous force chains, and a map of the mean stress (bottom). The damage distribution is plotted in Figure 6. Several examples of deformed and damaged bonds are shown in the enlarged boxes, and illustrate the heterogeneity of the local deformation. Even if at the macroscopic
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Top: Map of the (Von-Mises) equivalent stress (in MPa) that describes force chains in the cemented packing at the final stage of compression. Bottom: Map of the mean stress (in MPa).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Top: Map of the damage distribution at the final stage of compression. Bottom: Zoom in boxes with the cumulated plastic strain on a color scale.</figcaption></figure>
scale we deal with small strain-driven compression of the sample, at the microscopic scale the bonds are submitted to combinations of different large deformations: compression, traction, rotation, and shear. This is due to the geometrical complexity of the cemented granular material. In almost all cases the bond deformation is localized in a shear band that acts in each damaged bond. The orientation of these shear bands is also very heterogeneous.

Even at this low level of compaction we found that new contacts between the grains and the walls are established, but no new contact between the grains appear during the loading. That means that the P-wave propagation is not facilitated by these new inter-grains contacts as discussed in the context of “non-classical nonlinearity” or “clapping interface model” [Chaboche, 1992, Lemaitre and Desmorat, 2005, Pecorari and Solodov, 2006, Lyakhovsky et al., 2009].

### 3. Elastic wave probing

In this section, we present the simulations of the wave propagation, to highlight the wave speed sensitivity to the damage state of the sample. We want to show that it is possible to monitor the macroscopic damage through the wave speed record. For this, we model the propagation of small amplitude, purely reversible, non-destructive waves in the damaged material at different stages of the quasi-static loading. These stages are chosen to better understand the impact of large and small increments of damage and different regimes of the loading process (i.e., elastic at the beginning, then inelastic with increased damage).

#### 3.1. Problem setting

##### 3.1.1. Macroscopic level

At different stages $\epsilon_{\text{macro}}^*$, which correspond to the time $t^* \in [0, T^{\text{quasi-static}}]$ of the quasi-static loading process described in the previous section (indicated by green vertical lines in Figure 4), we analyze the wave propagation in the sample $\mathcal{D}_{\text{macro}}^* = \mathcal{D}_{t^*}^{\text{macro}}$. At the zero acoustic time (here $t$ represents the acoustic time, while the real time is $t^* + t$), a single sinusoidal longitudinal overstress pulse $\sigma_{\text{macro}}^*(t)$ (the total stress is $\sigma_{\text{macro}}^*(t^*) + \sigma_{\text{macro}}^*(t)$) is emitted at the top of the sample ($y = l$) in the $y$-direction and

propagates along the $y$-axis (see Figure 7). The transmitted P-wave is recorded on the bottom of the sample ($y = 0$). On the lateral walls, we impose frictionless contact boundary conditions. Two series of simulations are made for two frequencies of the pulse. The first one for 150 kHz, which is a low frequency for this sample as the corresponding wavelength is 20 mm in the glass, almost 50 times the mean radius of the beads. The second simulation is for 1.50 MHz, which we designate as a “medium” frequency as it corresponds to a wavelength of 2 mm in the glass, which is close to five times the mean radius of the beads. As expected, high-frequency simulations (at 7.50 MHz, corresponding to a wavelength of 0.4 mm, as the mean radius of the beads) show that the pulse does not (or hardly) propagate any more due to intense scattering, and these are not investigated here.

Note that in the seismological context, and for fault monitoring, the seismic velocity monitoring relies on relatively low frequency waves (less than a few Hz) due to the rapid attenuation of the waves in the complex, highly fractured rocks of the fault core. This means that the probe has a wavelength much larger than the typical grain size.

The macroscopic computations of the P-wave speed $V_P^*$ can be used to compute the macroscopic damage $d_{\text{macro}}^* = 1 - M_{\text{macro}}^* / M_{\text{macro}}^0$ (here $M_{\text{macro}}$ is the macroscopic uniaxial strain modulus). We know that $V_P^* = \sqrt{M_{\text{macro}}^* / \rho_{\text{macro}}^*}$ and by mass conservation, that $\rho_{\text{macro}}^0 = (1 - \epsilon_{\text{macro}}^*) \rho_{\text{macro}}^*$; then, we get the formula

$$d_{\text{macro}}^* = 1 - \left( \frac{V_P^*}{V_P^0} \right)^2 \frac{1}{1 - \epsilon_{\text{macro}}^*} \quad (1)$$

##### 3.1.2. Microscopic level

At different stages of the quasi-static loading process described in the previous section, and as indicated by green vertical lines in Figure 4, we stored the deformed meshes and the damage distributions. Let us denote the domain occupied by the cemented granular material by $\mathcal{D}^* = \mathcal{D}_{t^*}$ (with the partition into the phase grains $\mathcal{G}^* = \mathcal{G}_{t^*}$ and the phase bonds $\mathcal{B}^* = \mathcal{B}_{t^*}$), and the damage distribution computed in the quasi-static finite element simulations by $d^* = d_{t^*} : \mathcal{B}^* \to [0, 1]$.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Representation of the ultrasound probing numerical experiment.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** The initial mesh used for the dynamic computations.</figcaption></figure>

## 3.2. Mechanical and numerical modeling

As the waves have a small amplitude, there is no evolution of the damage and plastic strain during the dynamic process, which is said to be nonperturbative or nondestructive. This is why the small perturbation assumption can be considered to be valid and the material can be supposed to be linear elastic, with the mass density $\rho^*$ and the elastic coefficients $\lambda^*$, $\mu^*$, given by

$$ \rho^* = \begin{cases} \rho_g, & \text{in } \mathcal{G}^*, & \lambda^* = \begin{cases} \lambda_g, & \text{in } \mathcal{G}^*, \\ (1-d^*)\lambda_b, & \text{in } \mathcal{B}^* \end{cases} \\ \rho_b, & \text{in } \mathcal{B}^*, & \end{cases} \\ \mu^* = \begin{cases} \mu_g, & \text{in } \mathcal{G}^*, \\ (1-d^*)\mu_b, & \text{in } \mathcal{B}^*. \end{cases} $$

The dynamic problem is discretized in time by the classical implicit Newmark method, with $\beta = 1/4$ and $\gamma = 1/2$ [e.g., see Fung, 1997], which is unconditionally stable and thus allows much larger values of the time step than the critical Courant-Friedrichs-Lewy time step. For the space discretization, we used [P2-continuous] finite element discretization for the displacement field and [P1-discontinuous] for the stress field.

The quasi-static mesh used previously is not appropriate for wave propagation. For the dynamic wave propagation simulation, the domain $\mathcal{D}^*$ was remeshed to obtain a uniformly fine mesh (see Figure 8). Indeed, as the waves propagate everywhere in the sample, for the dynamic problem we need a
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Snapshots of the longitudinal stress $\sigma_{yy}$ field (as a color scale) at $t = 2 \ \mu s$ after the low-frequency excitation (top) and at $t = 2.64 \ \mu s$ after the medium-frequency excitation (top), in the undamaged sample (left) and after compression (right) corresponding to $\varepsilon_{\text{macro}}^* = 1.2\%$ and $d_{\text{ave}}(t^*) = 0.54\%$.</figcaption></figure>

regular mesh without any size contrast between the peripheral and the core edges of the beads. Since the open space results in a free boundary condition for the propagating waves, the beads size and the open spaces sizes are important in choosing the wavelength to probe the average material properties (see the discussion in the next section about the wave behavior for different frequencies).

### 3.3. Results

The low-frequency pulse (for which the wavelength in the glass is almost 50 times the mean radius of the beads) generates a well-established coherent wave, without remarkable multiply scattered waves. It is visible at $t = 2 \ \mu s$ after the excitation in the first line of Figure 9, which shows several snapshots of the longitudinal stress $\sigma_{yy}$ field in all of the sample. The left row shows the undamaged sample ($\varepsilon_{\text{macro}}^* = 0\%$ and $d_{\text{ave}}(t^*) = 0\%$), in which the coherent front has not yet reached the reception wall. The right row shows the same snapshot, but computed with a damaged sample ($\varepsilon_{\text{macro}}^* = 1.2\%$ and $d_{\text{ave}}(t^*) = 0.54\%$). This highlights a short delay of the coherent wave, which is especially visible in the beads close to the lower wall.

For the 1.50 MHz pulse (for which the wavelength in the glass is almost 5 times the mean radius of the beads), a similar delay can be observed between clear

coherent waves computed with the undamaged state and the same damaged state, as seen in the second line of Figure 9. However, in this case, the coherent waves are followed by multiply scattered waves, producing a tail in the signal similar to the seismic coda. Note the complexity of the field distribution at both the macroscopic scale and the bead scale, due to the high heterogeneity of the granular material.

The transmitted $\sigma_{yy}$-waves are recorded along the bottom wall (called the “reception wall”), as discretized in 1000 spatial points and shown for the 1.5 MHz pulse in Figure 10. They are gathered when they belong to the same bead of the reception side, making voids appear between the 26 beads of this side. One of these signals is plotted in Figure 11 (in green) for different states of damage of the sample. The sum of all these signals is visible as the blue curve inside the right box of Figure 10, which enhances the coherent arrivals. Indeed, unlike the coherent wave that is almost uniform along the y-axis, the codas depend on the recording position. So their sum, which is physically defined in the experiments by the large transducer in contact with all of the beads of the reception wall, tends to cancel, whereas the summed amplitudes of the coherent wave and its reflections are enhanced.

The signals plotted in Figure 11 are associated with different levels of damage as they correspond to the macroscopic strains $\varepsilon_{\text{macro}} = 0\%$, $\varepsilon_{\text{macro}} = 0.51\%$,
<figure><img src="image_14.png" /><figcaption>**Figure 10.** (a) Acoustic signals recorded along the reception wall (represented above) for a 1.5 MHz pulse. The green arrow indicates the position chosen to compare the codas in Figure 11. (b) Summed waves along the reception wall.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** The transmitted acoustic signals recorded for a 1.5 MHz pulse at the same position of the reception wall (indicated in green in Figure 10) for the undamaged sample (black) and after different steps of compression: $\varepsilon_{\text{macro}} = 0.51\%$ (red), $\varepsilon_{\text{macro}} = 0.54\%$ (green), and $\varepsilon_{\text{macro}} = 0.84\%$ (blue).</figcaption></figure>

$\varepsilon_{\text{macro}} = 0.54\%$ and $\varepsilon_{\text{macro}} = 0.84\%$. We observe that the second and third signals in Figure 11 are very similar, which is understandable, with the weak strain and the damage difference between them. The third signal in Figure 11 is slightly delayed compared to the second one, with a small amplitude loss, and this delay increases visibly with the acoustic time $t$ between 7 and 20 µs. After 20 µs, these signals are not coherent anymore. Similar, but more pronounced, observations can be made between these two signals and the first one computed with the undamaged sample. We note that the first coherent peak, which corresponds to the direct transmitted ballistic wave, has almost the same amplitude and time of flight for these three signals. This is not inconsistent, as at the same time

the thickness of the sample decreases and the wave velocity decreases. The contrast between the signals is much more obvious with the fourth signal computed, when the quasistatic loading begins to saturate after important damage to the sample. We observe a loss of amplitude of the transmitted ballistic wave and a time of flight increase. At later times, we observe important changes in the coda. These can be explained by the increasing scattering process of the waves due to two changes: the first one is the increasing impedance contrast between the glass and the damaged cement at the interface bead/bond linked to the degradation of the elastic parameters of the cement; the second one is the slight geometric evolution of the sample without true rearrangement,
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The wave velocity *V* in the sample during a quasistatic process versus the macroscopic strain ε<sub>macro</sub>. The velocity is calculated with the time of flight of the summed signals along the reception wall for a low-frequency pulse (150 kHz).</figcaption></figure>

as the contacts remain the same (i.e., no loss, no
creation).

The sum of the signals along the reception wall al-
lows estimation of the time of flight of the transmit-
ted pulse, and therefore of the velocity of the lon-
gitudinal coherent wave, as the ratio of the time of
flight and the thickness of the numerical sample at
the considered state of damage. We choose the time
of flight measured from the first peak. The veloc-
ities obtained are shown in Figure 12 for the low-
frequency pulse, to ensure that the coherent regime
is well established, and validate the approximation of
the effective medium theory [Digby, 1981]. The veloc-
ity can also be calculated from the mean time of flight
over the statistics of the signals recorded at 1000 posi-
tions along the reception wall, which gives very close
results, as can be seen in Figure 14b for a 1.5 MHz
pulse.

Figure 12 shows that for a low frequency pulse (150 kHz), the velocity first increases slightly in the elastic part of the deformation, when there is no damage (dV/V close to +1%); it then decreases dramatically (dV/V close to -10%). The first phase is associated to a slight increase of the wave velocity and can explained easily on the basis that at the beginning some additional disks enter into contact with the rigid piston, which creates new force chains that facilitate the wave propagation. The decrease in wave

velocity begins when the bonds are damaged; i.e.,
when the microscopic damage average d_ave starts
increasing significantly (see Figure 4). It is remark-
able that a very small amount of microscopic damage
average can imply such a dramatic loss. Moreover,
when the load reaches its stationary value, the veloc-
ity keeps decreasing, although by less, despite the al-
most constant microscopic damage average rate d_ave
(as also seen in Figure 4). This can be explained on
the basis that simulated waves propagate along the
deformation axis, which is also the privileged axis of
the first bond damaging, while in the last phase, the
microscopic damage concerns mainly the bonds ori-
entated mostly in the other directions (see Section 4
for more details). Moreover, it should not be forgotten
that in our simple mechanical model there is damage
saturation at d_max = 0.8, which implies that the veloc-
ity drop is also saturated in some sense.

In all of these considerations, it appears clear that
the microscopic damage average *d*<sub>ave</sub>, defined as the
ratio between the integral of the damage and the area
of the sample |∂<sub>*t*</sub>| and associated to the quasi-static
process, is not adapted to measure the damage of
the sample. Figure 13 shows the macroscopic dam-
age *d*<sub>macro</sub>, given by (1) and computed with the elas-
tic wave probing method. This overall damage pa-
rameter, which can be found without any computa-
tion of the quasistatic process, appears to be more
adapted to describing the complexity of the problem.
Indeed, for *ε*<sub>macro</sub> = 1% for instance, we found a value
of *d*<sub>macro</sub> that was almost 40 times greater than *d*<sub>ave</sub>,
and in the last phase associated to the load plateau,
the damage *d*<sub>macro</sub> also has a plateau, as expected.

To check the validity of the velocities calculated
with the time of flight in the case of a 1.5 MHz
pulse, we use the doublet technique [Poupinet et al.,
1984], which is also called the moving-window cross-
spectrum technique [Clarke et al., 2011] or coda wave
interferometry [Snieder et al., 2002]. We apply this to
all the consecutive couples *i* − 1, *i* corresponding to
the consecutive chosen quasistatic steps for the dy-
namical simulations. This ensures good temporal co-
herence between the signals *i* − 1 and *i*, as shown in
Figure 11, which is necessary to apply this technique.
This technique is especially interesting for subtle
changes of the material velocity, and it consists of
the evaluation of the relative velocity change δv<sub>ij</sub> =
(v<sub>j</sub> − v<sub>i</sub>)/v<sub>i</sub> of the medium at a given position be-
tween two temporal states *i* and *j*. We can also write
<figure><img src="image_17.png" /><figcaption>**Figure 13.** The macroscopic damage $d_{\text{macro}}$ computed with elastic wave probing (with a low-frequency pulse of 150 kHz) versus the macroscopic strain $\varepsilon_{\text{macro}}$.</figcaption></figure>

$\delta v_{ij} = \delta l_{ij} - \delta t_{ij}$ with $\delta l_{ij}$ linked to the strain (calculated from the quasistatic loading), $\delta t_{ij} = dt_{ij}/t$ as the relative change of the time of flight $dt_{ij}$ of the wave through the sample, and $t$ as the acoustic time, which is zero at the beginning of the excitation pulse. The procedure to evaluate $\delta t_{ij}$ has two steps: calculation of $dt_{ij}(t)$ for different times $t$ in the whole range of the acoustic time (so with the ballistic parts and the codas), as sufficient for the second step, which is linear regression of $dt_{ij}$ with $t$. This regression goes through the origin and its slope is directly $\delta t_{ij}$. Figure 14a shows the two contributions $\delta l_{ij}(t)$ and $\delta t_{ij}(t)$ as functions of the macroscopic deformation. The two contributions were calculated between consecutive pairs $i-1, i$ and were accumulated, to be compared to the initial step 0. The cumulative error from the last regression of the doublet technique is also plotted in Figure 14a. It increases with the macroscopic deformation, until almost 0.5% at the last step, as the nondeformed configuration is taken as the reference here.

Figure 14b shows the velocity change $\delta v_{ij}(t)$, which is the difference between these two contributions. The same error bars are also plotted. The velocity changes obtained with the times of flight and both of the summed signals and from their statistics (from the same dynamic simulations) are also plotted here. The results are very close.

4. Induced damage anisotropy

In this section, we investigate how the induced damage is related to the loading direction. We also investigate the loss of isotropy of the macroscopic wave propagation. As the effects are similar and more pronounced in the x-axis than in the y-axis direction, we consider the quasistatic edometric compression along the x-axis. For this simulation, the final strain is $\varepsilon_{\text{macro}} = 2.17\%$.

The “norm” of the deviatoric stress tensor is shown in Figure 15a, which corresponds to the last step of the computation. It shows that the force chains tend to be oriented along the axis of the loading. This assumes that the bonds oriented along this same axis are more solicited for the transmission of the stress, and as a consequence they would be more damaged than the others. This can be seen in Figure 15b, which shows the distribution of the damage parameter $d$ in the bonds at the last step. To rigorously demonstrate this observation, Figure 15c shows the value of the damage integrated over each bond independent of the others, as a function of its angle. Clear anisotropy appears and confirms the previous hypothesis: the damage is concentrated in the bonds with angles in $[\pi/6; 5\pi/6]$ modulo $\pi$.

Therefore we define the “x-bonds” and “y-bonds” as the two groups of bonds with angles between $[-\pi/6; \pi/6]$ and $[\pi/3; 2\pi/3]$ modulo $\pi$, respectively. Figure 15d shows the values of the damage $d_{\text{ave}}^x$ and $d_{\text{ave}}^y$ of the average damage $d_{\text{ave}}$ integrated over these two groups as functions of the macroscopic strain. The x-bonds concentrate always the greatest part of the total damage, and almost three times the damage of the y-bonds at the last step. The average damage initially increases quickly in the x-bonds, then keeps increasing, although more slowly, while it increases later in the y-bonds, with an intermediary and stable slope. First the x-bonds are mainly damaged and then the damage appears in the y-bonds, which are the last bonds possible to damage. This can explain why the normal load in the axis of the loading becomes stationary, whereas the total damage keeps increasing, as can be seen for compression along the y-axis in Figure 4.

These computations on the damage average suggest that the bounds orientated along the propagation axis are mainly implied in the speed wave drop. To confirm this assumption, we can test the
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Velocity changes computed with the doublet technique for a 1.5 MHz pulse. (a) Thickness (red) and delay (blue) changes as functions of the quasistatic step. Error bars on the delay come from the last regression of the doublet technique. (b) Velocity changes as a function of the quasistatic step and computed with the doublet technique (black), from the time of flight statistics (red) and the time of flight of the summed signals (blue).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Edometric compression along the x-axis. (a,b) Maps of the deviatoric stress and damage parameter *d*, respectively, at the last step of the quasi-static computation. (c) Damage integrated over the bonds independently as a function of their angle for different macroscopic strains. (d) Damage integrated over the “x-bonds” and “y-bonds” as a function of the macroscopic strain.</figcaption></figure>

anisotropy of the wave propagation by measuring the wave velocity along the x-axis and y-axis (called $V_x$ and $V_y$, respectively) in the configuration of this section with compression along the x-axis. We use low frequency pulses (150 kHz), as defined in Section 3, at six different steps of the same quasi-static compression simulation. With the associated wavelength of 20 mm, we consider that no multiple scattering occurs, which would considerably change the comparison between the two axes, as the wall is three times longer than the other one and would accumulate more scattering effects. The velocities are plotted in Figure 16 and are computed with the time of

flight of the first peak of the summed signals recorded along the corresponding wall and the thickness (constant for the propagation along the y-axis).

First, we check that $V_x$ and $V_y$ are almost the same at the initial state without compression (step 0), as there should not be any source of anisotropy at this step. We indeed find values that are close, as 2740 and 2690 m/s respectively, with a relative deviation of 1.8%. The small difference might be due to the geometry of the sample and the boundary conditions. Then we observe that both velocities are impacted by the loading and the induced damage, even $V_y$, and where $V_x$ shows a decrease of 7%, twice the decrease
<figure><img src="image_20.png" /><figcaption>**Figure 16.** $dV_x/V_x$ and $dV_y/V_y$ divided by their respective initial values as functions of the macroscopic deformation. The velocities are calculated with the time of flight of the summed signals along the reception wall for a low-frequency pulse (150 kHz).</figcaption></figure>

of $V_y$. This confirms the induced damage anisotropy for the wave propagation.

As they are calculated with the same method, we can also compare $V_x$ and $V_y$ to $V$, the velocity obtained from the compression and the wave propagation along the same y-axis represented in Figure 12a. $V$ and $V_x$ are very similar, with very close values and the same trend for $\varepsilon_{\text{macro}}$ between 0.3% and 1.5%, with maximal velocity losses $dV/V$ and $dV_x/V_x$ close to 10%. However, unlike $V_x$, $V_y$ and $V$ show a slight increase for short strains, which is due to the shorter size along the y-axis.

**5. Conclusion and discussion**

We have presented a series of numerical experiments to discuss the phenomenology responsible for the changes in seismic velocities during deformation of a cohesive granular medium. Our model is simplified to highlight the most important elements, and of course to guarantee the numerical stability of our results. We have considered a 2D model of a cohesive granular medium in an edometer compression test configuration, for which only monotonic loading experiments have been considered. The damage is associated with the nonlinear behavior of the bonds, in the form of a plasticity law. During compression, the

deformation quickly becomes very heterogeneous, with a concentration on the force chains where the bonds undergo all of the deformation modes.

In practice, the propagation of waves allows the changes in elastic properties to be highlighted by measuring the propagation speed. Our simulations allow us to observe the velocity drop in the deformation domain where the model is not very sensitive to the maximum damage characteristics that we had to introduce. We can deduce from the propagation velocity an effective macroscopic damage parameter, which turns out to be about 40 times larger than the average damage, which indicates the strong sensitivity of the waves to the deformation and the predominance of the weakening of the force chain elements. We have also shown the anisotropic character of the velocity reduction controlled by the direction of the imposed deformation.

At the microscopic level the damage model used for the bonds is isotropic and involves a Poisson ratio (or equivalently a ratio between P and S-wave velocities) which is not affected by damage. In contrast, at the macroscopic level the model loses the initial isotropy and we have noticed the differences in the P-wave velocities according to the direction of propagation. It will be interesting to see how the S-wave speeds are affected and to analyze how the ratio between P and S-wave speeds depends on the propagation direction. We note that the evolution of the ratio between P and S-wave velocities during the damage evolution is still an observational challenge for seismology. The reason is mainly the difficult localization and spatial extension of the changes that make precise quantitative measurement of the evolution of this ratio very difficult with present-day techniques.

To compare our results with seismic observations, a first order observable is the velocity itself. The earth crust has been subject to damage for a long time and damage has accumulated and persisted in the regions of strongest deformations, as active fault zones. The observed velocity in the shallow crust is diminishing in the vicinity of faults [Zigone et al., 2015], which show the presence of damage. The model qualitatively predicts these observations. Furthermore, measurable temporal changes in velocity can be observed in the Earth in response to various processes, which range from earthquake shaking to rain and Earth tides. These perturbations correspond to very small deformations, typically ranging
from $10^{-6}$ to $10^{-8}$. They would not induce significant global changes, but can affect locally the wave velocity. Our numerical experiments are limited to monotonically compressions and we do not have the precision to analyze the deformation levels encountered in seismology.

We have nevertheless shown high sensitivity of the wave velocity with the deformation during the damaging process for a cohesive granular medium as it has been reported before [Langlois and Jia, 2014]. We can try to put into perspective the simulations for our simplified model and the observations.

We can deduce from the simulations the slope of the velocity-deformation curve (e.g., Figure 12), which measures the sensitivity of the wave velocity. The curve shows an increase in sensitivity with damage. Natural observations do show significant variations in sensitivity as a function of rock conditions. This has been clearly shown on a regional scale in Japan [Brenguier et al., 2014], with active, highly fractured volcanic structures being the regions where the largest velocity drops are observed after large earthquakes.

The model we have studied allows us to highlight the importance of stress heterogeneity and localized nonlinearity at the microscopic scale for macroscopic behavior. Nevertheless, the model is very simple in terms of its geometry, and corresponds to a case where the jumps are few, leaving an important part to the elastic granular behavior. The parameters will have to be adapted so that the simulations correspond quantitatively to seismological observations. The complexity of geophysical environments, with multi-scale heterogeneities and the likely presence of fluids, makes the exercise impossible with the numerical means implemented here.

Based on Figure 12, the sensitivity of the relative velocity to deformation is of the order of 10 for our model in compressional experiments. Natural observations are for much smaller deformations than what is resolved in our simulations. It should be noted, however, that the observed sensitivity is much greater than that calculated in this study. For tectonically driven static expansions of the order of $10^{-6}$, a sensitivity of the order of 100 has been reported [Rivet et al., 2014]. This suggests very high spatial concentration of force chains and nonlinearity. In the extreme case of the response to Earth tides, the deformation is of the order of $10^8$ and the sensitivity

for very shallow materials would reach 10,000 [Mao et al., 2019b, Sens-Schönfelder and Eulenfeld, 2019]. The differences are significant, but the increase in sensitivity for the model can be produced by changing the geometry and the distribution of the elastic elements within the phenomenology described here, although at the cost of heavy computational effort. The most important aspect for comparison is the 2D limitation of the model. We have seen in 2D that the ratio between the macroscopic damage which governs the effective wave velocity and the mean damage is about 40, due to the concentration of the stresses on the force chains. In the transition to the 3D case, this concentration towards linear chains might result in a different (larger or equal) ratio. This leads to sensitivities of the order of some of those observed in Nature. The highest sensitivities to transient or periodic perturbations suggest a model for which the macroscopic response is controlled by very sparse and very localized chains of forces that are extremely sensitive to small macroscopic deformation perturbations. These conditions will be considered in further numerical studies.

## Conflicts of interest

Authors have no conflict of interest to declare.

## Acknowledgements

The authors acknowledge support from the European Research Council under the European Union Horizon 2020 research and innovation program (grant agreement no. 742335, F-IMAGE).

IRI acknowledges support of the Romanian Ministry of Research (project number PN-III-P4-PCE-2021-0921 within PNCDI III).

## Appendix. Mechanical model description

The movement (flow) in the Eulerian description is given by the velocity field, denoted $v(t, \cdot) : \mathcal{D}_t \to \mathbb{R}^d$. The rate of deformation and the spin rate are denoted by $D = D(v) = (\nabla v + \nabla^T v)/2$ and by $W = W(v) = (\nabla v - \nabla^T v)/2$, respectively. If we denote the Cauchy stress tensor by $\sigma(t, \cdot) : \mathcal{D}_t \to \mathbb{R}^{n \times n}$, then the momentum balance law for the quasi-static loading reads:

$$ \mathrm{div} \, \sigma(t) + \rho \mathbf{b}(t) = 0, \quad \text{in } \mathcal{D}_t, $$
where $\rho$ is the mass density and **b** are the body forces.

For the grains, we consider an isotropic hypo-elastic law [e.g., see Belytschko et al., 2000]

$$ \boldsymbol{\sigma}^{\nabla}(t) = \lambda_g \text{ trace } (\boldsymbol{D}\boldsymbol{I} + 2\mu_g \boldsymbol{D}), \text{ in } \mathcal{G}_t, $$

where $\lambda_g, \mu_g$ are the elastic coefficients of the grains, $\boldsymbol{\sigma}^{\nabla} = \dot{\boldsymbol{\sigma}} - \boldsymbol{W}\boldsymbol{\sigma} - \boldsymbol{\sigma}\boldsymbol{W}$ is the Jaumman time derivative, and $\dot{\boldsymbol{\sigma}} = \partial_t \boldsymbol{\sigma} + \boldsymbol{v} \cdot \nabla \boldsymbol{\sigma}$ is the material (total) time derivative.

For the mechanical modeling of the bond material, we chose a very simple elasto-plastic model, which couples isotropic ductile damage with Von-Mises plasticity. For this, we considered the additive decomposition of the rate deformation tensor into the elastic $\boldsymbol{D}^e$ and plastic rates $\boldsymbol{D}^p$ of deformation

$$ \boldsymbol{D} = \boldsymbol{D}^e + \boldsymbol{D}^p. $$

The weakening effect is characterized by the (phenomenological) damage parameter $d \in (0, 1)$. Following the strain equivalence principle proposed by Lemaitre [Lemaitre and Plumtree, 1979], the strain of damaged material ($d > 0$) subjected to the effective stress is the same as that of material without damage ($d = 0$). This means that the Lamé elastic coefficients have to be replaced by $(1-d)\lambda_b$ and $(1-d)\mu_b$, respectively, where $\lambda_b, \mu_b$ are the elastic coefficients of the undamaged bond material (see Table 1). For the elastic range, as for the grains, we considered the generalization of Hooke's law written in terms of the Jaumann rate of the Cauchy stress tensor

$$ \boldsymbol{\sigma}^{\nabla}(t) = (1-d)\lambda_b \text{ trace } (\boldsymbol{D}^e)\boldsymbol{I} + 2(1-d)\mu_b\boldsymbol{D}^e, \text{ in } \mathcal{B}_t. $$

The plastic rate of deformation is related to the Cauchy stress tensor through the flow rule associated to the classical Von-Mises yield criterion with no hardening. To be more precise, let $\mathcal{F}(\boldsymbol{\sigma}, d) = \sigma_{eq}/(1-d) - k$ be the yield function, with $k$ the yield limit and $\sigma_{eq} = \sqrt{2/3}|\boldsymbol{\sigma}^D|$ the Von-Mises stress ($\boldsymbol{\sigma}^D = \boldsymbol{\sigma} - (1/3)\text{trace}(\boldsymbol{\sigma})\boldsymbol{I}$ is the stress deviator). If we denote the accumulated plastic strain by $\boldsymbol{\epsilon}^p$ (given through the differential equation $\boldsymbol{\epsilon}^p = \sqrt{2/3}|\boldsymbol{D}^p|)$, then the flow rule reads

$$ \boldsymbol{D}^p = \frac{\boldsymbol{\epsilon}^p}{\sigma_{eq}}\boldsymbol{\sigma}^D, $$

while the loading-unloading conditions read

$$ \boldsymbol{\epsilon}^p \ge 0, \quad \mathcal{F}(\boldsymbol{\sigma}, d) \le 0, \quad \boldsymbol{\epsilon}^p \mathcal{F}(\boldsymbol{\sigma}, d) = 0. $$

We have considered here a very simple basic damage law [following Lemaitre and Chaboche, 1994],

where the damage is related only to the (accumulated) plastic strain $\boldsymbol{\epsilon}^p$.

$$ d = d^p(\boldsymbol{\epsilon}^p). $$

For simplicity the dependence of the damage on the accumulated plastic strain will be a piecewise linear function that involves an activation plastic strain $\boldsymbol{\epsilon}_{\text{activation}}^p$ (damage threshold) and a maximal level of damage $d_{\text{max}}$, corresponding to the plastic strain $\boldsymbol{\epsilon}_{\text{max}}^p$ (which might represent the rupture plastic strain), by

$$ d^p(\boldsymbol{\epsilon}^p) = \begin{cases} x \frac{d_{\text{max}}}{\boldsymbol{\epsilon}_{\text{max}}^p - \boldsymbol{\epsilon}_{\text{activation}}^p} [\boldsymbol{\epsilon}^p - \boldsymbol{\epsilon}_{\text{activation}}^p]_+, & \text{if } \boldsymbol{\epsilon}^p < \boldsymbol{\epsilon}_{\text{max}}^p, \\ d_{\text{max}}, & \text{if } \boldsymbol{\epsilon}^p \ge \boldsymbol{\epsilon}_{\text{max}}^p, \end{cases} $$

where $[x]_+ = (x + |x|)/2$ is the positive part function. Even if the constants $\boldsymbol{\epsilon}_{\text{activation}}^p, \boldsymbol{\epsilon}_{\text{max}}^p$ and $d_{\text{max}}$ depend generally on the tri-axiality ratio, in what follows we will assume that they are constants.

## References

Allam, A. A. et al. (2014). Seismic velocity structure in the hot springs and trifurcation areas of the San Jacinto fault zone, California, from double-difference tomography. *Geophys. J. Int.*, 198(2), 978–999. ISSN: 1365-246X, 0956-540X.

Barajas, A. et al. (2021). Separation of poroelastic and elastic processes of an aquifer from tectonic phenomena using geodetic, seismic, and meteorological data in the Pollino Region, Italy. *Geochem. Geophys. Geosyst.*, 22(11), article no. e2021GC009742. ISSN: 1525-2027, 1525-2027.

Belytschko, T., Liu, W. K., and Moran, B. (2000). *Non-linear Finite Elements for Continua and Structures*. John Wiley & Sons, Ltd, Chichester.

Brenguier, F., Campillo, M., Hadziioannou, C., et al. (2008). Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations. *Science*, 321(5895), 1478–1481. ISSN: 0036-8075, 1095-9203.

Brenguier, F., Campillo, M., Takeda, T., et al. (2014). Mapping pressurized volcanic fluids from induced crustal seismic velocity drops. *Science*, 345(6192), 80–82. ISSN: 0036-8075, 1095-9203.

Canel, V. et al. (2020). Monitoring of damage processes in cemented granular materials with acoustic emissions and seismic velocity reduction. Technical Report EGU2020-21891, Copernicus Meetings.
Cappa, F. et al. (2014). Off-fault long-term damage: A condition to account for generic, triangular earthquake slip profiles. *Geochem. Geophys. Geosyst.*, 15(4), 1476–1493. ISSN: 15252027.

Chaboche, J.-L. (1992). Damage induced anisotropy: On the difficulties associated with the active/passive unilateral condition. *Int. J. Damage Mech.*, 1(2), 148–171. ISSN: 1056-7895.

Chester, F. M. (1993). Internal structure and weakening mechanisms of the San Andreas fault. *J. Geophys. Res.*, 98(B1), 771–786.

Clarke, D. et al. (2011). Assessment of resolution and accuracy of the moving window cross spectral technique for monitoring crustal temporal variations using ambient seismic noise. *Geophys. J. Int.*, 186(2), 867–882. ISSN: 0956-540X.

Delorey, A. A. et al. (2021). Probing the damage zone at Parkfield. *Geophys. Res. Lett.*, 48(13), article no. e2021GL093518. ISSN: 1944-8007.

Digby, P. J. (1981). The effective elastic moduli of porous granular rocks. *J. Appl. Mech.*, 48(4), 803–808. ISSN: 00218936.

Dvorkin, J., Nur, A., and Yin, H. (1994). Effective properties of cemented granular materials. *Mech. Mater.*, 18(4), 351–366. ISSN: 0167-6636.

Evesque, P. (2000). Eléments de mécanique quasi-statique des milieux granulaires mouillés ou secs. *Poudres et grains*, NS(1), 1–155.

Fung, T. C. (1997). Unconditionally stable higher-order Newmark methods by sub-stepping procedure. *Comput. Methods Appl. Mech. Eng.*, 147(1), 61–84. ISSN: 0045-7825.

Ghosh, S. and Raju, S. (1996). R–S adapted arbitrary Lagrangian–Eulerian finite element method for metal-forming problems with strain localization. *Int. J. Numer. Methods Eng.*, 39(19), 3247–3272. ISSN: 1097-0207.

Hall, S. A. et al. (2008). Crack density tensor inversion for analysis of changes in rock frame architecture. *Geophys. J. Int.*, 173(2), 577–592. ISSN: 0956-540X.

Hamiel, Y. et al. (2009). Brittle deformation and damage-induced seismic wave anisotropy in rocks. *Geophys. J. Int.*, 178(2), 901–909. ISSN: 1365-246X.

Hemmerle, A., Schröter, M., and Goehring, L. (2016). A cohesive granular material with tunable elasticity. *Sci. Rep.*, 6(1), 1–11. ISSN: 2045-2322.

Huerta, A. and Casadei, F. (1994). New ALE applications in non-linear fast-transient solid dynamics. *Eng. Comput.*, 11(4), 317–345. ISSN: 0264-4401.

Johnson, P. A. and Rasolofosaon, P. N. J. (1996). Non-linear elasticity and stress-induced anisotropy in rock. *J. Geophys. Res. Solid Earth*, 101(B2), 3113–3124. ISSN: 2156-2202.

Johnson, P. A., Zinszner, B., and Rasolofosaon, P. N. J. (1996). Resonance and elastic nonlinear phenomena in rock. *J. Geophys. Res. Solid Earth*, 101(B5), 11553–11564. ISSN: 2156-2202.

Langlois, V. and Jia, X. (2014). Acoustic probing of elastic behavior and damage in weakly cemented granular media. *Phys. Rev. E*, 89(2), article no. 023206.

Lemaitre, J. and Chaboche, J.-L. (1994). *Mechanics of Solid Materials*. Cambridge University Press, Cambridge. ISBN: 978-0-521-47758-1.

Lemaitre, J. and Desmorat, R. (2005). *Engineering Damage Mechanics. Ductile, Creep, Fatigue and Brittle Failure*. Springer, Berlin, Heidelberg.

Lemaitre, J. and Plumtree, A. (1979). Application of damage concepts to predict creep-fatigue failures. *J. Eng. Mater. Technol.*, 101(3), 284–292. ISSN: 0094-4289.

Lyakhovsky, V., Ben-Zion, Y., and Agnon, A. (1997). Distributed damage, faulting, and friction. *J. Geophys. Res. Solid Earth*, 102(B12), 27635–27649. ISSN: 01480227.

Lyakhovsky, V. et al. (2009). Non-linear damage rheology and wave resonance in rocks. *Geophys. J. Int.*, 178(2), 910–920. ISSN: 0956540X, 1365246X.

Mao, S., Campillo, M., et al. (2019a). High temporal resolution monitoring of small variations in crustal strain by dense seismic arrays. *Geophys. Res. Lett.*, 46(1), 128–137. ISSN: 1944-8007.

Mao, S., Mordret, A., et al. (2019b). On the measurement of seismic travel-time changes in the time-frequency domain with wavelet cross-spectrum analysis. *Geophys. J. Int.*, 221(1), article no. ggz495. ISSN: 0956-540X, 1365-246X.

Margerin, L. et al. (2016). Sensitivity kernels for coda-wave interferometry and scattering tomography: Theory and numerical evaluation in two-dimensional anisotropically scattering media. *Geophys. J. Int.*, 204(1), 650–666. ISSN: 0956-540X, 1365-246X.

Meier, U., Shapiro, N. M., and Brenguier, F. (2010). Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise. *Geophys. J. Int.*, 181(2), 985–996. ISSN: 0956540X, 1365246X.
Nur, A. and Simmons, G. (1969). Stress-induced velocity anisotropy in rock: An experimental study. *J. Geophys. Res.* (1896–1977), 74(27), 6667–6674. ISSN: 2156-2202.

Obermann, A. et al. (2016). Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media. *Geophys. J. Int.*, 207(1), 59–66. ISSN: 0956-540X, 1365-246X.

Olivier, G. et al. (2015). Investigation of coseismic and postseismic processes using in situ measurements of seismic velocity variations in an underground mine. *Geophys. Res. Lett.*, 42(21), 9261–9269. ISSN: 1944-8007.

Pasqualini, D. et al. (2007). Nonequilibrium and nonlinear dynamics in berea and fontainebleau sandstones: Low-strain regime. *J. Geophys. Res. Solid Earth*, 112(B1), article no. B01204. ISSN: 2156-2202.

Pecorari, C. and Solodov, I. (2006). Nonclassical nonlinear dynamics of solid surfaces in partial contact for NDE applications. In Delsanto, P. P., editor, *Universality of Nonclassical Nonlinearity: Applications to Non-Destructive Evaluations and Ultrasonic*, pages 309–326. Springer, New York, NY. ISBN: 978-0-387-35851-2.

Poupinet, G., Ellsworth, W. L., and Frechet, J. (1984). Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California. *J. Geophys. Res. Solid Earth*, 89(B7), 5719–5731. ISSN: 2156-2202.

Radjai, F. and Dubois, F. (2011). *Discrete-Element Modeling of Granular Materials*. Wiley-Iste, London.

Reasenberg, P. and Aki, K. (1974). A precise, continuous measurement of seismic velocity for monitoring *in situ* stress. *J. Geophys. Res.* (1896–1977), 79(2), 399–406. ISSN: 2156-2202.

Rivet, D. et al. (2014). Seismic velocity changes, strain rate and non-volcanic tremors during the 2009–2010 slow slip event in Guerrero, Mexico. *Geophys. J. Int.*, 196(1), 447–460. ISSN: 0956-540X.

Rodríguez-Ferran, A., Pérez-Foguet, A., and Huerta, A. (2002). Arbitrary Lagrangian–Eulerian (ALE) formulation for hyperelastoplasticity. *Int. J. Numer. Methods Eng.*, 53(8), 1831–1851. ISSN: 1097-0207.

Roux, P. et al. (2016). A methodological approach towards high-resolution surface wave imaging of the

San Jacinto fault zone using ambient-noise recordings at a spatially dense array. *Geophys. J. Int.*, 206(2), 980–992. ISSN: 0956-540X, 1365-246X.

Sawicki, A. and Swidzinski, W. (1995). Cyclic compaction of soils, grains and powders. *Powder Technol.*, 85(2), 97–104. ISSN: 0032-5910.

Sayers, C. M. and Kachanov, M. (1991). A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics. *Intl J. Solids Struct.*, 27(6), 671–680. ISSN: 0020-7683.

Sayers, C. M. and Kachanov, M. (1995). Microcrack-induced elastic wave anisotropy of brittle rocks. *J. Geophys. Res. Solid Earth*, 100(B3), 4149–4156. ISSN: 2156-2202.

Schaff, D. P. and Beroza, G. C. (2004). Coseismic and postseismic velocity changes measured by repeating earthquakes: Coseismic and postseismic velocity changes. *J. Geophys. Res. Solid Earth*, 109(B10), article no. B10302. ISSN: 01480227.

Schubnel, A. and Guéguen, Y. (2003). Dispersion and anisotropy of elastic waves in cracked rocks. *J. Geophys. Res. Solid Earth*, 108(B2), article no. 2101. ISSN: 2156-2202.

Scuderi, M. M. et al. (2016). Precursory changes in seismic velocity for the spectrum of earthquake failure modes. *Nat. Geosci.*, 9(9), 695–700. ISSN: 1752-0908.

Sens-Schönfelder, C. and Eulenfeld, T. (2019). Probing the *in situ* elastic nonlinearity of rocks with earth tides and seismic noise. *Phys. Rev. Lett.*, 122(13), article no. 138501. ISSN: 0031-9007, 1079-7114.

Sens-Schönfelder, C. and Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. *Geophys. Res. Lett.*, 33(21), article no. L21302. ISSN: 1944-8007.

Snieder, R. et al. (2002). Coda wave interferometry for estimating nonlinear behavior in seismic velocity. *Science*, 295(5563), 2253–2255. ISSN: 0036-8075, 1095-9203.

Stanchits, S., Vinciguerra, S., and Dresen, G. (2006). Ultrasonic velocities, acoustic emission characteristics and crack damage of basalt and granite. *Pure Appl. Geophys.*, 163(5), 975–994. ISSN: 1420-9136.

Takano, T. et al. (2014). Seismic velocity changes caused by the earth tide: Ambient noise correlation
analyses of small-array data. *Geophys. Res. Lett.*, 41(17), 6131-6136. ISSN: 1944-8007.

TenCate, J. A., Smith, E., and Guyer, R. A. (2000). Universal slow dynamics in granular solids. *Phys. Rev. Lett.*, 85(5), 1020-1023. ISSN: 0031-9007, 1079-7114.

Touma, R. et al. (2022). Distribution of seismic scatterers in the San Jacinto fault zone, southeast of Anza, California, based on passive matrix imaging. *Earth Planet. Sci. Lett.*, 578, article no. 117304. ISSN: 0012-821X.

van Dinther, C., Margerin, L., and Campillo, M. (2021). Laterally varying scattering properties in the North Anatolian fault zone from ambient noise cross-correlations. *Geophys. J. Int.*, 225(1), 589-607. ISSN: 0956-540X.

Wang, J. and Gadala, M. S. (1997). Formulation and

survey of ALE method in nonlinear solid mechanics. *Finite Elem. Anal. Des.*, 24(4), 253-269. ISSN: 0168-874X.

Wang, Q.-Y. et al. (2017). Seasonal crustal seismic velocity changes throughout Japan. *J. Geophys. Res. Solid Earth*, 122(10), 7987-8002. ISSN: 2169-9356.

Zhang, T., Sens-Schönfelder, C., and Margerin, L. (2021). Sensitivity kernels for static and dynamic tomography of scattering and absorbing media with elastic waves: A probabilistic approach. *Geo-phys. J. Int.*, 225(3), 1824-1853. ISSN: 0956-540X, 1365-246X.

Zigone, D. et al. (2015). Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and love waves. *Pure Appl. Geophys.*, 172(5), 1007-1032. ISSN: 0033-4553, 1420-9136.