Vincent Canel
email: canelvince@gmail.com

Michel Campillo
email: michel.campillo@univ-grenoble-alpes.fr

Xiaoping Jia
email: xiaoping.jia@espci.fr

Ioan R Ionescu
email: ioan.r.ionescu@gmail.com

New Developments in Passive Seismic Imaging and Monitoring / Nouveaux développements dans le domaine de l'imagerie et de la surveillance sismique passive Damage in cohesive granular materials: simulations and geophysical implications

Keywords: Granular medium, Damage, Elastic wave, Simulation, Earthquake

The aim of this paper is to test a simple damage model of a cohesive granular medium to study the relationship between the damage and velocity of elastic waves. Our numerical experiments of edometric compression show that the microscopic deformation quickly becomes very heterogeneous, while our simulations of elastic waves propagation show that a small amount of damage induces a dramatic decrease in the elastic velocity. This shows that cohesive discrete media are very sensitive to strain field heterogeneity, and that the wave velocities in these media can measure subtle transient deformation processes, such as earthquake initiation phases.

Introduction

The damage state of rocks in a fault core is a controlling element of seismic rupture and of the characteristics of fault evolution [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF][START_REF] Cappa | Off-fault long-term damage: A condition to account for generic, triangular earthquake slip profiles[END_REF]. The interaction between slip and damage processes is poorly understood, especially as the nature and actual level of rock damage at depth are also poorly known despite some recent studies [START_REF] Delorey | Probing the damage zone at Parkfield[END_REF]. On a geological time scale (i.e., over several seismic cycles), fault zones evolve, with the development of damage structures in their vicinity, partly due to dynamic stresses during rupture [START_REF] Chester | Internal structure and weakening mechanisms of the San Andreas fault[END_REF]. The existence of this damaged zone is traditionally shown at depth by seismic tomography, which indicates a lower velocity zone that extends several kilometers in depth [START_REF] Zigone | Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and love waves[END_REF][START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto fault zone using ambient-noise recordings at a spatially dense array[END_REF][START_REF] Allam | Seismic velocity structure in the hot springs and trifurcation areas of the San Jacinto fault zone, California, from doubledifference tomography[END_REF]. More recently, it has been shown that the existence of a several kilometers thick zone of strong scattering around the North Anatolian Fault is required to explain the regional distribution of multiply scattered coda wave energies [START_REF] Van Dinther | Laterally varying scattering properties in the North Anatolian fault zone from ambient noise cross-correlations[END_REF]. Passive scatterer imaging was made possible by a novel aberration correction method [START_REF] Touma | Distribution of seismic scatterers in the San Jacinto fault zone, southeast of Anza, California, based on passive matrix imaging[END_REF] that confirmed the presence of intense fracturing at depth around the San Jacinto fault in California. Indeed, measurements of the temporal evolution of seismic wave velocities give information on the mechanical state of the medium [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[END_REF][START_REF] Schaff | Coseismic and postseismic velocity changes measured by repeating earthquakes: Coseismic and postseismic velocity changes[END_REF][START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]. Noise-based passive monitoring allows to envision long-term continuous monitoring of the seismic velocity [START_REF] Brenguier | Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations[END_REF], 2014[START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF][START_REF] Rivet | Seismic velocity changes, strain rate and non-volcanic tremors during the 2009-2010 slow slip event in Guerrero, Mexico[END_REF]. 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 [START_REF] Wang | Seasonal crustal seismic velocity changes throughout Japan[END_REF] and very short time resolution [Mao et al., 2019b, Sens-Schönfelder and[START_REF] Sens-Schönfelder | Probing the in situ elastic nonlinearity of rocks with earth tides and seismic noise[END_REF].

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 andWegler, 2006, Barajas et al., 2021], temperature [START_REF] Meier | Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise[END_REF], or tides [START_REF] Sens-Schönfelder | Probing the in situ elastic nonlinearity of rocks with earth tides and seismic noise[END_REF]. As the causes of external forcing are fairly well known, correction strategies have been proposed [START_REF] Wang | Seasonal crustal seismic velocity changes throughout Japan[END_REF]. Assessing these changes in the vicinity of rupture zones will improve through the development of sensitivity kernels for spatial inversion methods [START_REF] Margerin | Sensitivity kernels for coda-wave interferometry and scattering tomography: Theory and numerical evaluation in two-dimensional anisotropically scattering media[END_REF][START_REF] Obermann | Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media[END_REF][START_REF] Van Dinther | Laterally varying scattering properties in the North Anatolian fault zone from ambient noise cross-correlations[END_REF][START_REF] Barajas | Separation of poroelastic and elastic processes of an aquifer from tectonic phenomena using geodetic, seismic, and meteorological data in the Region, Italy[END_REF][START_REF] Zhang | Sensitivity kernels for static and dynamic tomography of scattering and absorbing media with elastic waves: A probabilistic approach[END_REF].

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 [START_REF] Brenguier | Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations[END_REF] and transient tectonic deformations [START_REF] Rivet | Seismic velocity changes, strain rate and non-volcanic tremors during the 2009-2010 slow slip event in Guerrero, Mexico[END_REF] has been shown in the form of a velocity drop followed by a slow recovery similar to laboratory observations [START_REF] Tencate | Universal slow dynamics in granular solids[END_REF]. These changes are also observed at the time scale of days in a mining environment [START_REF] Olivier | Investigation of coseismic and postseismic processes using in situ measurements of seismic velocity variations in an underground mine[END_REF] and in laboratory experiments of stick-slip and slow slip [START_REF] Scuderi | Precursory changes in seismic velocity for the spectrum of earthquake failure modes[END_REF], 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 [START_REF] Nur | Stress-induced velocity anisotropy in rock: An experimental study[END_REF] that could be linked to the elastic nonlinear behavior of rocks [Johnson et al., 1996, Johnson and Rasolofosaon, 1996[START_REF] Pasqualini | Nonequilibrium and nonlinear dynamics in berea and fontainebleau sandstones: Low-strain regime[END_REF]. Direct comparisons between rock degradation (linked to a damage quantified by the distribution of microcracks) and wave velocity measurements could be highlighted experimentally [START_REF] Hamiel | Brittle deformation and damage-induced seismic wave anisotropy in rocks[END_REF]. 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 [START_REF] Sayers | A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics[END_REF], 1995[START_REF] Schubnel | Dispersion and anisotropy of elastic waves in cracked rocks[END_REF][START_REF] Stanchits | Ultrasonic velocities, acoustic emission characteristics and crack damage of basalt and granite[END_REF][START_REF] Hall | Crack density tensor inversion for analysis of changes in rock frame architecture[END_REF]. For the following numerical simulations, we consider granular materials rather than real rocks as they have been studied as synthetic rocks in previous experiments [START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]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 [START_REF] Tencate | Universal slow dynamics in granular solids[END_REF] observed in the earth after a strong drop [START_REF] Brenguier | Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations[END_REF].

Here we use a simple cemented granular media model [START_REF] Dvorkin | Effective properties of cemented granular materials[END_REF][START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF][START_REF] Hemmerle | A cohesive granular material with tunable elasticity[END_REF] 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 [START_REF] Evesque | Eléments de mécanique quasistatique des milieux granulaires mouillés ou secs[END_REF][START_REF] Sawicki | Cyclic compaction of soils, grains and powders[END_REF][START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF], 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 [START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]. 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 [START_REF] Radjai | Discrete-Element Modeling of Granular Materials[END_REF], 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.

Quasistatic loading

Macroscopic setting

Let us first describe the numerical edometric experiment. The "sample" (or cell) of the cemented granular material initially occupies the 2D domain D macro 0 = (0, L) × (0, l 0 ) (with L the constant length and l 0 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 and then the macroscopic displacement ∆l = |l -l 0 | evolve linearly with time. The microscopic deformation process is heterogeneously activated by the macroscopic deformation ε macro = ∆l /l 0 , which increases slowly from 0 to the final strain ε f macro = 1.58%. During the loading process, to find the macroscopic stress σ macro = F macro /L, we compute the resultant force F macro that acts on the upper plate. In this configuration, the macroscopic Hooke's law at the macroscopic level is written as σ macro = M macro ε macro , where M macro is the macroscopic uniaxial strain modulus or P-wave modulus (M = λ + 2µ = E (1 -ν)/((1 + ν)(1 -2ν)) with λ and µ the macroscopic first and second Lamé coefficients and E and ν the macroscopic Young and Poisson moduli respectively).

Microscopic setting

The cemented granular material, which occupies at the moment t ∈ [0, T quasi-static ] the domain D t ⊂ D macro t is composed of grains, denoted by G t , and bonds, denoted by B t , i.e. Dt = Ḡt ∪ Bt and G t ∩B t =

. For the initial geometric settings of grains G 0 we used 382 disks (to model glass beads) of radii 300 µm to 600 µm. The beads are connected by the cement B 0 composed of 956 bonds and gathered in a cell D macro 0 of size 3 cm × 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 ε f macro . 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 probability distribution of the orientation angles of the bonds (i.e., angle between the (Ox) 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 π-periodic: each bond is indeed indexed on the basis of both its associated bead, with an angle θ, and the other associated bead, with an angle θ + π.

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 symetrical 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 The index 0 refers to the state without damage.

there is no correlation between length and orientation angle.

Mechanical modeling

For the mechanical modeling of the bond material, which is detailed in Appendix, we chose the model of [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF] coupling isotropic ductile plastic damage and with an elasto-plastic law in the framework of an Eulerian (large deformations) description [see for instance [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF].

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 λ and µ written in terms of the Jaumann 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 ∈ [0, 1] and its evolution law is related to the cumulated plastic strain ϵ 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 [START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]Jia, 2014, Canel et al., 2020]. The elastic coefficients of the undamaged material are given in Table 1 and the yield limit is κ 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 ϵ p,activation = 2%, we consider the beginning of damage process, while the maximal level of damage is d max = 0.8 corresponding to ϵ 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 [START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]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.

Numerical modeling

As the bonds are submitted to large deformation (even for small macroscopic deformation), we used an arbitrary Lagrangian Eulerian (ALE) description [START_REF] Huerta | New ALE applications in non-linear fast-transient solid dynamics[END_REF][START_REF] Wang | Formulation and survey of ALE method in nonlinear solid mechanics[END_REF][START_REF] Ghosh | R-S adapted arbitrary Lagrangian-Eulerian finite element method for metal-forming problems with strain localization[END_REF][START_REF] Rodríguez-Ferran | Arbitrary Lagrangian-Eulerian (ALE) formulation for hyperelastoplasticity[END_REF]. 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 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 +∆t , we used a classical implicite backward radial return mapping scheme [e.g., see [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] 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-discontinous] 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 and8). 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.

Results

Macroscopic level

Figure 4 shows the macroscopic normal load σ macro versus the imposed strain on the sample ε 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 0 macro ≈ 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 strainstress curve for an elasto-plastic material without hardening.

However, the macroscopic damage d macro of the sample, defined as

d macro = 1 - M macro M 0 macro
, 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 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[START_REF] Reasenberg | A precise, continuous measurement of seismic velocity for monitoring in situ stress[END_REF][START_REF] Takano | Seismic velocity changes caused by the earth tide: Ambient noise correlation analyses of small-array data[END_REF][START_REF] Delorey | Probing the damage zone at Parkfield[END_REF]. 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 ave (t ), defined as the ratio between the integral of the damage, that vanishes outside the bonds, and the area of the sample |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 |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 ave (t ) keeps increasing (which can be associated to anisotropic effects of the compression; see Section 4 for more details).

Microscopic level

Even for small macroscopic deformation, local plastic irreversible deformation occurs in the bonds. The first microscopic plastic and damage effects occur at ε macro = 0.096% and ε 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 (ε f macro = 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 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 heterogenous.

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 intergrains contacts as discussed in the context of "nonclassical nonlinearity" or "clapping interface model" [START_REF] Chaboche | Damage induced anisotropy: On the difficulties associated with the active/passive unilateral condition[END_REF][START_REF] Lemaitre | Engineering Damage Mechanics. Ductile, Creep, Fatigue and Brittle Failure[END_REF][START_REF] Pecorari | Nonclassical nonlinear dynamics of solid surfaces in partial contact for NDE applications[END_REF][START_REF] Lyakhovsky | Non-linear damage rheology and wave resonance in rocks[END_REF].

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

Problem setting

Macroscopic level

At different stages ε * macro , which correspond to the time t * ∈ [0, T 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 D * macro = D macro t *

. At the zero acoustic time (here t represents the acoustic time, while the real time is t * + t ), a single sinusoidal longitudinal overstress pulse σ * macro (t ) (the total stress is σ macro (t * ) + σ * 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 mean 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 * macro = 1 -M * macro /M 0 macro (here M macro is the macroscopic uniaxial strain modulus). We know that V * P = M * macro /ρ * macro and by mass conservation, that ρ 0 macro = (1 -ε * macro )ρ * macro ; then, we get the formula

d * macro = 1 - V * P V 0 P 2 1 1 -ε * macro .
(1)

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

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 ρ * and the elastic coefficients λ * , µ * , given by

ρ * = ρ g , in G , ρ b , in B * , λ * = λ g , in G * , (1 -d * )λ b , in B * , µ * = µ g , in G * , (1 -d * )µ b , in B * .
The dynamic problem is discretized in time by the classical implicit Newmark method, with β = 1/4 and γ = 1/2 [e.g., see [START_REF] Fung | Unconditionally stable higherorder Newmark methods by sub-stepping procedure[END_REF], 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 [P2continuous] finite element discretization for the displacement field and [P1-discontinous] for the stress field.

The quasi-static mesh used previously is not appropriate for wave propagation. For the dynamic wave propagation simulation, the domain 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

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 µs after the excitation in the first line of Figure 9, which shows several snapshots of the longitudinal stress σ y y field in all of the sample. The left row shows the undamaged sample (ε * macro = 0% and d 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 (ε * macro = 1.2% and d 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 σ y y -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 ε macro = 0%, ε macro = 0.51%, 10) for the undamaged sample (black) and after different steps of compression: ε macro = 0.51% (red), ε macro = 0.54% (green), and ε macro = 0.84% (blue). ε macro = 0.54% and ε 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, as the contacts remain the same (i.e., no loss, no creation).

The sum of the signals along the reception wall allows estimation of the time of flight of the transmitted pulse, and therefore of the velocity of the longitudinal 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 velocities obtained are shown in Figure 12 for the lowfrequency pulse, to ensure that the coherent regime is well established, and validate the approximation of the effective medium theory [START_REF] Digby | The effective elastic moduli of porous granular rocks[END_REF]. The velocity can also be calculated from the mean time of flight over the statistics of the signals recorded at 1000 positions 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 begining 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 remarkable that a very small amount of microscopic damage average can imply such a dramatic loss. Moreover, when the load reaches its stationary value, the velocity keeps decreasing, although by less, despite the almost constant microscopic damage average rate ḋ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 orientated 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 velocity drop is also saturated in some sense.

In all of these considerations, it appears clear that the microscopic damage average d ave , defined as the ratio between the integral of the damage and the area of the sample |D t | and associated to the quasi-static process, is not adapted to measure the damage of the sample. Figure 13 shows the macroscopic damage d macro , given by (1) and computed with the elastic wave probing method. This overall damage parameter, which can be found without any computation of the quasistatic process, appears to be more adapted to describing the complexity of the problem. Indeed, for ε macro = 1% for instance, we found a value of d macro that was almost 40 times greater than d ave , and in the last phase associated to the load plateau, the damage d macro 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 [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[END_REF], which is also called the moving-window crossspectrum technique [START_REF] Clarke | Assessment of resolution and accuracy of the moving window cross spectral technique for monitoring crustal temporal variations using ambient seismic noise[END_REF] or coda wave interferometry [START_REF] Snieder | Coda wave interferometry for estimating nonlinear behavior in seismic velocity[END_REF]. We apply this to all the consecutive couples i -1, i corresponding to the consecutive chosen quasistatic steps for the dynamical simulations. This ensures good temporal coherence 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 i j = (v jv i )/v i of the medium at a given position between two temporal states i and j . We can also write Figure 13. The macroscopic damage d macro computed with elastic wave probing (with a low-frequency pulse of 150 kHz) versus the macroscopic strain ε macro . δv i j = δl i j -δt i j with δl i j linked to the strain (calculated from the quasistatic loading), δt i j = dt i j /t as the relative change of the time of flight dt i j 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 δt i j has two steps: calculation of dt i j (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 i j with t . This regression goes through the origin and its slope is directly δt i j . Figure 14a shows the two contributions δl i j (t ) and δt i j (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 δv i j (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.

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 ε 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 [π/6; 5π/6] modulo π.

Therefore we define the "x-bonds" and "y-bonds" as the two groups of bonds with angles between [-π/6; π/6] and [π/3; 2π/3] modulo π, respectively.

Figure 15d shows the values of the damage d x ave and d y ave of the average damage d 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 yaxis 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 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 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 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 ε 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.

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 [START_REF] Zigone | Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and love waves[END_REF], 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 monotonic 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 [START_REF] Langlois | Acoustic probing of elastic behavior and damage in weakly cemented granular media[END_REF]. 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 [START_REF] Brenguier | Mapping pressurized volcanic fluids from induced crustal seismic velocity drops[END_REF], 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 [START_REF] Rivet | Seismic velocity changes, strain rate and non-volcanic tremors during the 2009-2010 slow slip event in Guerrero, Mexico[END_REF]. 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[START_REF] Sens-Schönfelder | Probing the in situ elastic nonlinearity of rocks with earth tides and seismic noise[END_REF]]. 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.

where ρ is the mass density and b are the body forces.

For the grains, we consider an isotropic hypoelastic law [e.g., see [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] σ ∇ (t ) = λ g trace (D)I + 2µ g D, in G t , where λ g , µ g are the elastic coefficients of the grains, σ ∇ = σ -W σ -σW is the Jaumman time derivative, and σ = ∂ t σ + v • ∇σ 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 D e and plastic rates D p of deformation

D = D e + D p .
The weakening effect is characterized by the (phenomenological) damage parameter d ∈ (0, 1). Following the strain equivalence principle proposed by Lemaitre [START_REF] Lemaitre | Application of damage concepts to predict creep-fatigue failures[END_REF], 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 )λ b and (1-d )µ b , respectively, where λ b , µ 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 We have considered here a very simple basic damage law [following [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF],

where the damage is related only to the (accumulated) plastic strain ε p . d = d p (ε 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 ε p activation (damage threshold) and a maximal level of damage d max , corresponding to the plastic strain ε p max (which might represent the rupture plastic strain), by

d p (ε p ) =      x d

Figure 1 .

 1 Figure 1. Representation of the quasi-static edometric numerical experiment.

Figure 2 .

 2 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 θ of a bond.

Figure 3 .

 3 Figure 3. The initial mesh used for the quasi-static computations.

Figure 4 .

 4 Figure 4. (a) Macroscopic stress σ macro (black; in × 33 MPa) and the averaged damage d ave (blue; in %) as functions of the strain ε 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.

Figure 5 .

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

Figure 6 .

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

  , we stored the deformed meshes and the damage distributions. Let us denote the domain occupied by the cemented granular material by D * = D t * (with the partition into the phase grains G * = G t * and the phase bonds B * = B t * ), and the damage distribution computed in the quasi-static finite element simulations by d * = d t * : B * → [0, 1].

Figure 7 .

 7 Figure 7. Representation of the ultrasound probing numerical experiment.

Figure 8 .

 8 Figure 8. The initial mesh used for the dynamic computations.

Figure 9 .

 9 Figure 9. Snapshots of the longitudinal stress σ y y field (as a color scale) at t = 2 µs after the lowfrequency excitation (top) and at t = 2.64 µs after the medium-frequency excitation (top), in the undamaged sample (left) and after compression (right) corresponding to ε * macro = 1.2% and d ave (t * ) = 0.54%.

Figure 10 .

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

Figure 11 .

 11 Figure11. The transmitted acoustic signals recorded for a 1.5 MHz pulse at the same position of the reception wall (indicated in green in Figure10) for the undamaged sample (black) and after different steps of compression: ε macro = 0.51% (red), ε macro = 0.54% (green), and ε macro = 0.84% (blue).

Figure 12 .

 12 Figure 12. The wave velocity V in the sample during a quasistatic process versus the macroscopic strain ε macro . The velocity is calculated with the time of flight of the summed signals along the reception wall for a low-frequency pulse (150 kHz).

Figure 14 .

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

Figure 15 .

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

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

Figure 16 .

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

σ

  ∇ (t ) = (1d )λ b trace (D e )I + 2(1d )µ b D e , in 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 F (σ, d ) = σ eq /(1d )k be the yield function, with k the yield limit and σ eq = 2/3|σ D | the Von-Mises stress (σ D = σ -(1/3)trace(σ)I is the stress deviator). If we denote the accumulated plastic strain by ε p (given through the differential equation εp = 2/3|D p |), then the flow rule readsD p = εp σ eq σ D ,while the loading-unloading conditions read εp ≥ 0, F (σ, d ) ≤ 0, εp F (σ, d ) = 0.

Table 1 .

 1 Parameters for simulations

		ρ (kg/m 3 ) E 0 (GPa) ν 0 (GPa) M 0 (GPa) c p0 (m/s) c s0 (m/s)
	Beads	2600	60	0.24	70.72	5200	3100
	Bonds	1200	3.5	0.33	5.18	2100	1000

  + = (x + |x|)/2 is the positive part function. Even if the constants ε max depend generally on the tri-axiality ratio, in what follows we will assume that they are constants.

	max p max -ε p activation d max , ε	[ε p -ε	p activation ] + , if ε p < ε if ε p ≥ ε	p max , p max ,
	where [x] p activation , ε

p max and d

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

Conflicts of interest

Authors have no conflict of interest to declare.

Appendix. Mechanical model description

The movement (flow) in the Eulerian description is given by the velocity field, denoted v(t , •) :

The rate of deformation and the spin rate are denoted by D = D(v) = (∇v + ∇ T v)/2 and by W = W (v) = (∇v -∇ T v)/2, respectively. If we denote the Cauchy stress tensor by σ(t , •) : D t → R n×n , then the momentum balance law for the quasi-static loading reads: div σ(t ) + ρb(t ) = 0, in D t ,