# Elastic Strain Effects on Wave Scattering: Implications for Coda Wave Interferometry (CWI)

Jerome Azzola, Jean Schmittbuhl, Dimitri Zigone, Olivier Lengliné, Frédéric Masson, Vincent Magnenet

► To cite this version:

Jerome Azzola, Jean Schmittbuhl, Dimitri Zigone, Olivier Lengliné, Frédéric Masson, et al. Elastic Strain Effects on Wave Scattering: Implications for Coda Wave Interferometry (CWI). Journal of Geophysical Research: Solid Earth, 2020, 125 (3), 10.1029/2019JB018974. http://dx.doi.org/10.1029/2019JB018977

HAL Id: hal-02550077

https://hal.science/hal-02550077

Submitted on 29 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

JGR Solid Earth

RESEARCH ARTICLE

10.1029/2019JB018974

**Key Points:**

* Impact of elastic deformation on CWI at laboratory scales is quantified following an experimental approach and an analog numerical model

* Relative contributions to travel time changes of local density changes from volumetric strain, of deformation of scatterers, and of acoustoelastic effects are quantified

* The numerical modeling of the strain influence on CWI highlights the dominant influence of acoustoelastics effects

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

J. Azzola,
jerome.azzola@unistra.fr

**Citation:**

Azzola, J., Schmittbuhl, J., Zigone, D., Lengliné, O., Masson, F., & Magnenet, V. (2020). Elastic strain effects on wave scattering: Implications for Coda Wave Interferometry (CWI). *Journal of Geophysical Research: Solid Earth, 125*, e2019JB018974. https://doi.org/10.1029/2019JB018974

Received 30 OCT 2019
Accepted 18 FEB 2020
Accepted article online 20 FEB 2020

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

# Elastic Strain Effects on Wave Scattering: Implications for Coda Wave Interferometry (CWI)

Jérôme Azzola¹ ID, Jean Schmittbuhl¹ ID, Dimitri Zigone¹ ID, Olivier Lengliné¹, Frédéric Masson¹ ID, and Vincent Magnenet² ID

¹Institut de Physique du Globe de Strasbourg/EOST, University of Strasbourg/CNRS, Strasbourg, France, ²ICUBE, University of Strasbourg/CNRS, Illkirch, France

**Abstract** Coda Wave Interferometry (CWI) is a highly sensitive monitoring technique built on the sensitivity of elastic coda waves to small changes in a diffusive medium. However, a clear connection between the physical processes involved in the evolution of the medium and the time changes observed by CWI has not been clearly described yet. Here, we quantify the impact of elastic deformation on CWI measurements at laboratory scales. We compare experimental results from wave scattering measurements during a uniaxial compression test to those of a numerical approach based on the combination of two codes (SPECFEM2D and *Code_Aster*), which allows us to model wave propagation in complex diffusive media during its elastic deformation. In both approaches, the reversible time delays measured between waveforms increase with the elastic deformation of the sample. From the numerical modeling, we gain insight to the relative contributions of different physical effects on the CWI measurement: local density changes from volumetric strain, the deformation of scatterers, and acoustoelastic effects. Our results suggest that acoustoelastics effects related to nonlinear elasticity are dominant.

## 1. Introduction

Imaging and monitoring techniques based on the arrival times of direct or scattered waves (using different types of waves like electromagnetic, acoustic, or elastic waves) are applied in many domains like medical imaging, nondestructive testing of materials, or seismic imaging with nominal scales ranging from mm to km (Kuperman et al., 2002; Mussett & Khan, 2000). The issue of monitoring complex media with noninvasive techniques has grown in interest in recent years and led to the development of innovative acoustic monitoring techniques (Anugonda et al., 2001; Wegler & Lühr, 2001). In such complex media, the energy of the direct arrivals is redistributed to multiply scattered waves with complex trajectories. When recorded, they typically correspond to the tail of the waveforms and are often referred as “coda” after the Italian word for tail. Multiply scattered waves sample the medium more extensively than direct or singly scattered waves because of the complexity of their paths (Sato et al., 2012). However, in highly heterogeneous media, scattered waves might be strongly attenuated which limit classical imaging and monitoring techniques (Aki & Richards, 2002). The phenomenon is common in natural systems such as volcanoes or geothermal reservoirs and also in anthropogenic structures made of concrete or rocks.

Coda waves are much more sensitive to changes occurring in the medium than direct arrivals. Measurements of phase or travel time variations in time, in multiple-scattered waves, can be exploited as a diagnostic of the weak perturbation in the medium. For instance, the comparison of the coda of seismic waves obtained at two different times from a fixed source enabled to monitor weak velocity changes in the Earth crust (Poupinet et al., 1984). This technique, called seismic doublet technique, has later been formalized under the denomination of Coda Wave Interferometry or CWI (e.g., Snieder, 2006). More recently, the use of the pervasive ambient seismic noise as a versatile seismic source has led to the development of Ambient Noise Interferometry, a technique that allows to reconstruct virtual active sources at every passive station, by cross-correlating continuous seismic noise records between each pair of stations (see Campillo and Roux (2015) for a review). This method not only revolutionized the use of seismic arrays to image the upper crustal structures of the Earth (Boué et al., 2014; Lin et al., 2013; Moschetti et al., 2007; Shapiro et al., 2005; Zigone et al., 2015, 2019) but also opened perspectives for monitoring the time evolution of structures with CWI. For example, small velocity changes have been monitored in natural structures like fault zones (Brenguier, Campillo, et al.,

<footnote>©2020 The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution-NonCommercial License,
which permits use, distribution and
reproduction in any medium, provided
the original work is properly cited and
is not used for commercial purposes.</footnote>
2008; Poupinet et al., 1984) or volcanoes where seismic velocity changes are highlighted a few days prior to eruptions (Ratdomopurbo & Poupinet, 1995; Sens-Schönfelder & Wegler, 2006). The monitoring has also been applied to geothermal reservoirs (Meunier et al., 2001).

Laboratory experiments on representative samples are necessary to understand the way CWI measurements are influenced by a local or global perturbation of the medium and to develop relevant monitoring techniques. The influence of various phenomena on wave velocity variations has been investigated. For instance, seismic velocity changes have been shown to be related to water saturation, air pressure, stress, or temperature changes (Grêt et al., 2006; Lockner et al., 1977; Planès & Larose, 2013; Snieder et al., 2002). For the case of concrete, using CWI with ultrasound, Planès & Larose, 2013) studied the impact on scattered wave propagation from global changes (thermal perturbations or acoustoelastic effects due to uniaxial stress) and from local changes (damage development which introduces new scatterers corresponding to local structural changes).

Laboratory experiments combined with numerical simulations are essential for the physical interpretation of CWI measurements (Azzola, Griffiths, et al., 2018). Indeed, their interpretation in terms of disruptive events is challenging since the measurements usually enclose at once all the contributions of the various physical sensitivities of the scattered waves. Deciphering a specific effect remains difficult from an experimental approach because of the structural complexity of scatterers in the material and the strong coupling between processes (e.g., microfracturing related to nonhomogeneous thermoelastic dilatation). Laboratory studies have already evidenced the effect of stress variations on velocity changes in strong scattering media such as concrete. Relative velocity variations are detected under various stress states whether it is for a few MPa with a relative resolution of $10^{-3}$ to $10^{-4}$ (Grêt et al., 2006) or under weaker load variations (50 kPa) (Larose & Hall, 2009). Otherwise, numerical studies report the effect of structural changes on the CWI measurements (Planès et al., 2014). But much less attention has been paid to the strain dependence of scattered wave velocities in heterogeneous media, and to our knowledge, very few studies address the effect of a reversible mechanical deformation on the scattered wavefield propagated within the medium.

In this work, we focus on the effect of elastic stress/strain deformation on CWI measurements. We based our study on the fine comparison of an extended experimental approach and an advanced numerical modeling. Both include the same parts: seismic wave scattering and mechanical deformation. They aim at developing a comparative forward modeling of the effect on the scattering elastic waves (i.e., the coda), of stress changes (from 5 up to 50 MPa), in the elastic regime, during a uniaxial compression test of the same sample. From the reproduction of similar effects, by both approaches, we intend to capture in a heuristic approach, the effect of nonlinear elasticity on wave scattering. On the one hand, one expects the experiments to include all the contributions of the processes involved and their complex coupling. On the other hand, the flexibility of the numerical simulation enables to tune parameters or effects and to decipher the different physical phenomena at the origin of the measurement, like the relative contributions of local density changes in the volumetric strain, the impact of the scatter deformation, or the acoustoelastic effects related to nonlinear elasticity. The combination of an experimental approach with an analog numerical model (Azzola, Griffiths, et al., 2018) opens up perspectives to assess the relative contributions of deformation and velocity changes during CWI measurements.

Section 2 of the present manuscript introduces the problem we address: the effect of the reversible deformation, in the elastic domain, on the scattered wavefield. We first detail the experimental method based on a simple synthetic sample. A decimetric block is made of homogeneous Aluminum 2017A (Au4G) (i.e., Duraluminium) but includes heterogeneities as a set of cylindrical holes designed to be responsible for a quasi-2-D multiply scattered elastic wavefield. The block is mechanically loaded using a uniaxial imposed displacement during the propagation of scattered elastic waves from a set of point sources. We then detail the complementary numerical method that we developed for a quantitative comparison, where elastic waves are propagated in a similar numerical sample with the same geometry, mechanical properties and propagating elastic wave velocities as the experimental block, in order to obtain the evolution of synthetic scattered waveforms during the same imposed deformation. The numerical method combines two approaches: a finite element method for the quasi-static deformation and a spectral element method for the wave propagation. Both numerical methods are adjusted to conserve the grid mesh during the whole simulation. In section 3, we compare first the wave scattering properties of the analog and laboratory samples, without any loading,
<figure><img src="image_4.png" /></figure>

**Figure 1.** Sketch describing the principle of the setup. The mechanical elastic deformation of the sample—a holed Au4G (Aluminium 2017A) block— is obtained from the step-by-step increase of the displacement $\delta$ applied along a 62 mm long “loading line,” that is, the contact surface between the sample and the piston of a servo-controlled uniaxial press (represented by a red line). A Ricker wavelet is sent at the middle of the left face of the sample at each step of the displacement. Waveforms are recorded at the opposing face with the receiver. Background colors show the trace of strain tensor computed by *Code_Aster* when a displacement of $\delta = 80$ µm is applied.

and second their mechanical responses during uniaxial load including the local measurement of the deformation field. We also include reversibility tests to check the elastic behavior of the sample. Section 4 presents the principle of CWI and the laboratory and numerical CWI measurements during deformation. This section also details, using the numerical approach, the effects of the sample boundaries deformation, of the local density changes related to the local volumetric strain, and of nonlinear elastic effects, that is, acoustoelastic effects. Section 5 discusses two main limitations of our approach: the sensitivity of high-order elastic coefficients and the 2-D approximation proposed in our simulations. Finally, conclusions are proposed in section 6.

## 2. Materials and Methods

### 2.1. Objectives of the Study

Figure 1 presents the principle of our approach, designed to study the sensitivity of elastic wave scattering to the reversible mechanical deformation of the medium. For this, we performed a large set of laboratory and numerical experiments on a decimetric elastic sample that includes a prescribed set of cylindrical wave scatterers. The loading consists of a displacement $\delta$ applied at the top of the sample that is increased step by step in the elastic regime. At each step of the loading, an impulsive acoustic wave is sent through the diffusive sample using a point source. The transmitted waves are recorded on the opposite face of the sample by a linear array of piezoelectric receivers. Waveform records include direct wave arrivals and also reflected and coda wave arrivals. The waveforms recorded for different loading states—referred to in the following as “perturbed” waveforms—are compared to the one recorded when the sample is unloaded—referred to in the following as “reference” waveform—using an interferometric technique: the CWI (Snieder, 2006).

The specificity here is to provide extended direct observations on wave scattering during the elastic deformation of the same sample through two different but complementary approaches: a laboratory experimental approach that is detailed in section 2.2 and a forward numerical approach to model the laboratory sample and experimental observations, as detailed in section 2.3. Laboratory measurements encompass on the one hand the contributions of all the distinctive processes involved in the sample's evolution

while numerical simulations allow to test the approximations necessary to reproduce the experimental observations and to explore separately the different physical processes at play.

### 2.2. Experimental Approach

#### 2.2.1. Laboratory Sample Description

Figure 2 presents the experimental setup used in the laboratory to measure changes in wave scattering induced by the mechanical loading of the sample. It shows the sample made of Aluminum 2017A (Au4G), commonly labeled Duraluminium or Dural, an aluminum, and copper alloy. The material has been chosen for its homogeneity at the scales considered for the study and its suitability for a precise design of the geometrical heterogeneities introduced inside. The size of the sample is 112 mm in height, 95 mm in width, and 32 mm thick.

To induce wave scattering, we introduced scatterers in the Au4G sample by drilling a set of 70 parallel cylindrical holes with a diameter of 6 mm through the whole thickness of the sample (see Figures 1 and 2). By changing the number and the size of the cylindrical inclusions, one expects to modify the scattering properties of the propagation medium. The design of the scatterers in the sample has been optimized numerically (see Azzola, Griffiths, et al., 2018) by adjusting the number of holes and their size in order to produce a multiply scattered wave field (see section 3). Considering the same geometry of the propagation
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Image of the experimental apparatus. The piston of a servo-controlled 10 T press loads a holed block of Au4G (70 holes of radius 3 mm). Two linear acoustic arrays of 32 sensors are held in direct contact with the boundaries of the sample. At each step of the displacement imposed by the piston, the element of the transducer on the left which is located at the middle of the array sends a pulse wave. This wave is recorded by each of the 32 elements transducer located on the opposite face of the sample. An image correlation technique is used to estimate the displacement field by comparing pictures of the sample taken at each step of the loading procedure.</figcaption></figure>

medium, the authors have tested different hole configurations and we use here the selected design of the sample. The length of the holes (32 mm, i.e., equal to the thickness of the sample) is significantly longer than the wavelength (5.8 mm) making the medium quasi-2-D in terms of wave scattering (Derode et al., 2001; Tourin et al., 2000). Au4G is a homogeneous and isotropic material, and its mechanical properties are well known. Young's modulus, Poisson ratio, and density provided by the manufacturer are, respectively, $E = 74$ GPa, $\nu = 0.33$, and $\rho = 2790$ kg m$^{-3}$. The compressive elastic limit of elasticity is $Re = 250$ MPa.

In the laboratory, we carried out several tests in order to ascertain the physical parameters supplied by the manufacturer. A typical uniaxial compressive test on a cylindrical sample of Au4G of radius 31 mm (similar to the loading piston radius) has been performed to check the value of Young's modulus. The macroscopic vertical strain $\epsilon_{YY}$—calculated from the relative variation of height of the sample—varies linearly with the vertical stress $\sigma_Y$—calculated by division of the force $F$ applied by the bearing surface of the piston $S$. The best linear regression fitted to the data points gives a slope of 74 GPa in accordance with Young's modulus provided by the manufacturer (see Figure S1 in supporting information). We also assessed the $P$ wave velocity in Duraluminium by measuring the first arrivals in an unstressed sample without scatterers and made of the same material (see Text S1). The $P$ wave velocity directly measured in an unstressed and scatterers-free Au4G sample, that is, $v_P = 6,209 \pm 47$ m s$^{-1}$, is in good agreement with the one calculated from quasi-static values of elastic moduli and density.

### 2.2.2. Laboratory Setup Description

Figure 2 also presents the mechanical loading of the sample. It is held between a vertical piston loaded by a servo-controlled 10 T press and a plate of Au4G attached to the force sensor below. The force applied by the piston is measured from a force sensor (REP-transducer of type TC4 with a maximum load capacity of 30 T) and is recorded using a NI-6133 data acquisition card. The maximal stress recorded in the experiments is 95 MPa.

Two linear acoustic arrays of 32 piezoelectric transducers (Imasonic) are placed on the opposing faces of the sample as shown in Figure 2. All the piezoelectric sensors have a resonance frequency of 500 kHz. To ensure an efficient acoustic coupling and to reduce the acoustic impedance contrast between the sensors and the sample, we used a high-viscosity Shear Wave Couplant gel from Sofranel. Both transducers are held in direct contact with the boundaries of the sample (Figure 2). Because of the fluid property of the coupling material, the sample deforms independently of the sensors. After each step increase of the displacement, we

send an electrical pulse at one of the transducers located at midheight. It has the temporal shape of a standard Ricker wavelet, that is, the second derivative of a Gaussian, whose standard definition is given for a dominant frequency $f_0$ in equation (1).

$$R(t) = [1 - 2[\pi f_0 t]^2] e^{-(\pi f_0 t)^2}, \quad (1)$$

In our approach, the central frequency is $f_0 = 400$ kHz (wavelength of 5.8 mm). The emission simultaneously triggers the acquisition at the 32 piezoelectric receivers of the opposite linear array which faces the source. The emission and acquisition are repeated 100 times at each step increase of the displacement to average the response and increase the signal to noise ratio. The acquisition time for the laboratory experiments is 0.4 ms with a sampling rate of 20 M samples per second using a 64 channel Lecoeur Electronic acquisition system.
### 2.2.3. Measurement of Displacement

The 2-D displacement field along one of the lateral sample surfaces (i.e., facing the camera) is measured using an image correlation technique. High-resolution pictures of the sample are taken at each step of the loading process with a Nikon D800 camera (Figure 2). They are compared to a reference one, acquired when the sample was unstressed. A time-resolved digital particle image velocimetry (DPIV) tool (PIV-lab, Thielicke & Stamhuis, 2014) is employed to compare by cross correlation the sample images recorded for different loading conditions.

The principle of DPIV is to divide a pair of images into multiple zones of interest and to find the translation vector that allows an optimal overlapping of both image tiles. This optimization is performed by calculating the correlation coefficient between the color levels of each zone. The displacement field between the two images is deduced by calculating a translation vector for each pair of image tiles. We base our approach on the use of 64 pixel square interrogation areas that are spaced by 10 pixels. The correlation matrix is computed in the frequency, which enables to achieve better computational costs than direct approaches based on the computation of the correlation matrix in the spatial domain. Subpixel precision is reached in the measurement of the displacement field using a Gauss algorithm for the interpolation.

For a homogeneous material such as Dural, the detection of the translation vector is potentially perturbed by the low texture of the image tiles, so the maximum correlation can be poorly resolved. Inhomogeneous lighting also affects the correlation signal (Raffel et al., 2018). The images are thus preprocessed to improve the accuracy of the vector detection. The aim of the preprocessing is to artificially generate texture in the image by amplifying the color gradients and to reduce the influence of oversaturated areas. Even if artificial, this texture is attributed to the material and not to lighting or other elements not related to the deformation of the sample. This preprocessing enables the maximum correlation to be better defined and thus increases the accuracy of the DPIV tool (Thielicke, 2014).

The resolution of the displacement field is conditioned by the resolution of the pictures and by the parameters of the DPIV routine. The minimum displacement measurable is less than the pixel length, as subpixel precision is reached through interpolation. The resolution of the pictures is here approximately 50 pixels per millimeter. The density of the displacement field is determined by the spacing between each square interrogation area, which is here of 10 pixels, that is, about 0.2 mm.

The variation in height of the sample is finally computed from the difference in mean vertical displacements between top and bottom boundaries of the sample. To account for the local fluctuations of the displacement field computed from the cross-correlation procedure, we estimate the mean displacement along the top and bottom boundaries by calculating the averaged vertical displacement along these lines.

## 2.3. Numerical Approach

The aim of the simulation is to replicate the laboratory experiments presented in section 2.2 to study the impact of the mechanical loading of the sample on wave scattering with a complementary technique. We describe here the numerical sample and the principle of the numerical scheme used to replicate the laboratory experiment.

### 2.3.1. Numerical Sample Description

The numerical sample aims to accurately reproduce both the macroscopic mechanical behavior and the wave scattering properties of the laboratory sample. Our numerical simulations are performed in two dimensions (2-D) both for the sample deformation and the wave propagation using the same sample dimensions (height and width) as the laboratory sample. We discuss in section 5.2 the influence of considering the full 3-D geometry of the sample. All physical parameters of the numerical sample are isotropic and similar to the laboratory sample. We also reproduce in the numerical sample the set of holes drilled in the Au4G sample with the same dimensions and same locations. The medium is discretized using GMSH—a GNU General Public License 3-D finite element mesh generator. We use a 2-D finite element mesh grid composed of 235,611 quadrilateral four-noded elements with a characteristic length $l_c = 0.5$ mm, which is the characteristic distance between two nodes of the mesh.

### 2.3.2. Numerical Scheme Description

The numerical scheme is the same as proposed by Azzola, Schmittbuhl, et al. (2018). It relies on the combination of a finite element deformation code (Code ASTER) and a spectral element wave propagation code
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Power spectral density (PSD) of (a) the unfiltered “reference” waveform obtained at the receiver from the simulation (blue line) and measured in the laboratory experiments (red line) with an unstressed sample, (b) after filtering in the (300–600) kHz frequency band.</figcaption></figure>

(SPECFEM2D). The mesh grid is iteratively deformed using the finite element software *Code_ASTER*, which is an open source general FEA software. Each step of the loading consists in the resolution of a static mechanical problem. The model assumes plane deformation. The bottom of the sample is fixed without any horizontal or vertical displacement. The right and left boundaries are stress free. We impose a given displacement, noted $\delta$, along a 62 mm long “loading line” centered at the top of the sample which transcribes in 2-D, the contact area between the piston and the laboratory sample. The displacement field at each step of the deformation process is used to deform subsequently the mesh grid.

The deformed mesh grid is then an input for SPECFEM2D (a 2-D spectral element code, Komatitsch & Vilotte, 1998) that is used to simulate the wave propagation within the medium. A Ricker wavelet of dominant frequency $f_0 = 400$ kHz (see equation (1)) is emitted at the source as a force point source acting perpendicularly to the sample boundary. The waveforms are sampled at the opposite side along 32 locations that mimic the laboratory 32 sensor array. The waveforms are computed with a sampling rate of 10 G samples per second, 500 times higher than the experimental sampling rate to limit numerical dispersion. Due to the difference in sampling rates, the original waveforms were resampled to ensure that the parameters of the two approaches were consistent.

In both the simulations and the experiments, the Ricker-shaped source function has a central frequency $f_0 = 400$ kHz. However, in the laboratory experiments, the actual Ricker source corresponds to the electric signal sent to the sensor and not the shape of the mechanical signal sent through the sample. The source emitted in the medium is actually a signal resulting from the convolution of the source function (i.e., the Ricker) by the transfer function of the sensor, which is peaked at 500 kHz. In the simulation, the transfer function of the emitting source is considered as perfect: The source function is applied without any convolution with a transfer function. To check that the source used in the simulations and the pulse emitted in the laboratory illuminate the medium in the same frequency range, we compared their respective power spectral density. Figure 3a shows the power spectral density computed for the “reference” seismogram recorded at the receiver, experimentally (red line) and numerically (blue line), without any filtering of the signals. Both spectra are very similar after filtering in this frequency range (300–600) kHz with a central frequency around 400 kHz for both sources (see Figure 3b). The similarity has been checked by computing the normalized correlation coefficient between both spectra displayed in Figure 3b. It reaches a value of 0.96.

# 3. Wave Scattering Properties of the Samples and Elastic Loading

## 3.1. Wave Field Scattering

The comparison of laboratory and simulated CWI measurements requires that the numerical model accurately reproduces the laboratory sample properties. Before comparing the mechanical behavior of the sample in the simulation and in the laboratory, we characterize their scattering properties, in order to verify that the behavior of the laboratory sample is retrieved in the simulation and also to check that both samples are in the strong scattering regime needed for the CWI analysis. The characterization is based on two measurements, each of them bringing arguments for the validation of the hypothesis related to the frame of the
study. We measure the temporal variation of the seismic energy for a wave propagation in the unstressed sample and estimate from these measurements the mean-free paths from a least squares inversion procedure using a diffusion model (Dainty & Toksoz, 1977; Dainty & Toksöz, 1981; Olivier et al., 2015).

### 3.1.1. Energy Density and Diffusion Model

The waveforms are recorded at each sensor of the receiver line in the unstressed sample until 0.4 ms, in the laboratory and in the numerical approach. From the measurement, we compute the temporal variation of the energy density $W(t)$ as a function of time, for a fixed distance $D$ between the source and the receiver:

$$W(t) = f(t)^2 + H[f[t]]^2, \quad (2)$$

where $f$ (in m s⁻¹) refers to the waveform measured at the receiver. $H$ is the Hilbert transform. This measurement characterizes the temporal variation of the seismic energy measured for fixed positions of the source and of the receiver.

Then, we estimate for comparison the mean-free path $ℓ$ in the numerical and the laboratory samples. It is the mean distance between two scattering events and a key parameter in the characterization of the scattering effect. We used a least squares minimization procedure to compare the temporal variation of the energy density to the one predicted by a diffusion model. The diffusion model expresses the temporal variation of the seismic energy of multiply scattered body waves, measured at a distance $D$ from the source as (Dainty & Toksoz, 1977; Dainty & Toksöz, 1981; Margerin, 2005; Trégourès & van Tiggelen, 2002):

$$Wt = E_0 \left( \frac{4 \pi v_s t \ell}{3} \right)^{-\frac{1}{2}} e^{-\left( \frac{v_s t}{\ell_a} - \frac{3 D^2}{4 v_s t \ell} \right)}, \quad (3)$$

where $E_0$ is the source energy, $v_s$ is the S wave velocity, $ℓ_a$ is the intrinsic absorption characteristic length, and $ℓ$ the mean-free path. We use a normalized logarithmic expression of the temporal variation of the seismic energy, in the same way as Olivier et al. (2015). By choosing an adequate reference in time $t_1$, so that $W_1(t_1) = 1$ J m⁻³ and by taking the natural logarithm of a referenced function $U(t)$ defined below, Wegler and Lühr (2001) show that equation (3) can be expressed as a sum of three dimensionless terms. Each of them consists in the product of a function of time ($t^0$, $t^1$, and $t^{-1}$) with a constant ($a_1$, $a_2$, and $a_3$, respectively) that depends on parameters of the model:

$$\ln U(t) = a_1 + a_2 t + a_3 \frac{1}{t}, \quad (4)$$

with $U(t) = (W(t)/W_1)(t/t_1)^{3/2}$, and $a_1 = (v_s t_1 / \ell_a + 3 D^2 / 4 v_s t_1 \ell) / \ell_a$; $a_2 = -v_s / \ell_a$; $a_3 = -3 D^2 / 4 v_s \ell$.

From the temporal variation measured in the laboratory and in the numerical sample, we compute this logarithmic and normalized expression of the energy density by multiplying equation (2) by the factor $t^{3/2}$ and by taking its natural logarithm, such that one reconstructs the left side of equation (4). The dimensionless coefficients $a_1$, $a_2$, and $a_3$ are measured in the same way as Olivier et al. (2015) by applying a least squares minimization to compare the observed measurements to the ones predicted by the diffusion model for body waves described by equation (3). The mean-free path is calculated from $a_3$ using $v_s = 3,158$ m s⁻¹ and the distance $D$ between the source and the receiver, which depends on the position of the sensor of the receiver line at which the methodology is applied.

### 3.1.2. Seismic Energy

The average temporal variation of the logarithmic and normalized expression of the function $W(t)$, $<\ln U(t)>$, is presented in Figure 4 for the experiment and the numerical simulation as continuous lines (respectively black and red curves). It is the spatial average of the functions computed from the response of each of the 32 receivers deployed along the side of the sample. The cross-correlation coefficient between both functions computed for a time window between 0 and 0.4 ms is 0.94. The high cross-correlation coefficient proves the strong agreement between the behavior of the numerical sample and the one of the laboratory sample. Both curves show a strong increase in the measured seismic energy up to a maximum (around 0.2 ms). It is followed by a progressive stabilization for the numerical simulation and a slow decrease for the laboratory experiments. The difference in behavior is explained by the fact that the intrinsic absorption of the
<figure><img src="image_7.png" /></figure>

**Figure 4.** Averages of the functions $\ln U(t)$ calculated from the responses of the 32 sensors distributed on the side of the sample, computed from laboratory measurements (plain black curve) and simulated waveforms (plain red curve), in an unstressed medium. $\ln U(t)$ is the temporal variation of the logarithmic expressions of the normalized energy density $U(t)$. The best fits to laboratory and numerical measurements determined by a least squares method using the proposed diffusion model (see text for details) are plotted as dashed curves (respectively dashed black and dashed red lines). To illustrate the spread in laboratory measurements, we have included a gray shading around the mean trend (black curve) measuring the standard deviation at each time sample. Experimental and numerical measurements show a different trend beyond 1–1.5 s.

material that exists in the experimental sample is not introduced in the simulation. In the following CWI measurements, we limit our coda window to 0.1 ms in order to stay in a regime where intrinsic attenuation is negligible.

### 3.1.3. Mean-Free Paths and Strong Scattering Regime

The best fit to the laboratory (dashed black) and numerical (dashed red) data, determined by the least squares method and using the diffusion model (equation (4)), is represented in Figure 4. Thirty-two distinctive estimates of the mean-free path are inverted using the records of each of the 32 experimental receivers, by applying the same computation to each receiver after adjusting the value of the distance to the source $D$. Variations among each set of estimates characterize the sensitivity of the measurement to the receiver position. We computed the average and the standard deviation of the set of estimates: $\ell = 11.5 \text{ mm} \pm 2 \text{ mm}$ for the experiments and $\ell = 11.0 \text{ mm} \pm 0.5 \text{ mm}$ for the numerical simulations. We observe that the mean-free paths inverted from laboratory and simulated records are similar within the error bars. From the similarity of the mean-free paths and the high cross correlation of the functions deduced from laboratory and numerical records represented in Figure 4, we conclude that the scattering properties of the numerical sample reproduce well those of the laboratory sample.

In a multiple-scattering regime, the wavefield emitted encounters the heterogeneities several times before being recorded by the receiver. It requires that the wavelength $\lambda$ of the incident wave is of same length or shorter than the size of the heterogeneities. Here, the diameter of the holes is 6 mm and the wavelength calculated from the wave velocities and frequencies used in the simulations is 5.8 mm. The condition neces-

sary to assess strong scattering is that the mean-free path $\ell$, the wavelength $\lambda$, the defect size $d$, and the distance from source to receivers $D$ satisfy the inequality: $\lambda < d \le \ell < D$ (Planès & Larose, 2013). In our simulations, we neglect attenuation effects and we do not consider in the inequality the intrinsic absorption characteristic length $\ell_a$ included in the inequality proposed by Planès and Larose (2013). The mean-free path inverted using simulated record ($\ell = 11.0 \text{ mm}$) or laboratory records ($\ell = 11.5 \text{ mm}$) satisfies the inequality. This result along with the progressive stabilization of the seismic energy highlighted from Figure 4 proves that the samples satisfy both numerically and experimentally, the multiple-scattering condition, a requirement necessary to apply CWI.

## 3.2. Quasi-Static Linear Elastic Deformation

Using *Code\_Aster*, we compute at each step of the loading the displacement of each node of the elements constitutive of the mesh grid. In the simulation, the displacement field has a spatial resolution at the scale of the length of each finite element, which is of the order of the characteristic length of the mesh grid—$l_c = 0.5 \text{ mm}$ (see Figure 5, right). The displacement field everywhere within a given element is discretized from the displacement computed by *Code\_ASTER* at each Gauss point and from the shape function of the elements. The components of the strain tensor at each Gauss point are calculated by computing the gradients of the displacements evaluated at these points. The strain components are finally obtained at each node by extrapolation, using the shape functions of the elements. The trace of the strain tensor extrapolated at each node of the mesh is displayed as background of Figure 1 for an applied displacement $\delta = 60 \mu\text{m}$. In Figure 5, the vertical component of the displacement $u_Y(X, Y)$ is plotted for an imposed displacement of $\delta = 60 \mu\text{m}$. It shows a weak influence of the holes and a limited effect of the boundary conditions. We compare this numerical computation to the experiments (Figure 5). The best resolution in the experiments is achieved by using cross correlation of images of the deformed sample with the image of the unstressed sample (see section 2.2.3). Left panel in Figure 5 shows the amplitude of the displacement field $u_Y(X, Y)$ obtained from the cross-correlation procedure using the same color panel to compare visually both images. Despite the local fluctuations in the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: experimental displacement field $u_Y(X, Y)$ calculated from image correlation between a picture of the unstressed sample and one captured after imposing a displacement of $\delta = 60$ µm to the top of the sample (corresponding to an applied force of $F = 75$ kN). Right: numerical displacement field $u_Y(X, Y)$ computed by Code_Aster when applying the same loading and using the same color scale. The holes introduced in the sample are represented by black circles with the same size and spatial locations. The red bar indicates the extent of the contact zone between the sample and the piston.</figcaption></figure>

experimental displacement field and a blurring effect around the holes (they appear larger than their real size), Figure 5 shows strong similarities in terms of amplitudes and spatial distribution between the simulated displacement field and the one calculated from laboratory experiments. The vertical displacement field computed from the laboratory measurements shows significant speckle pattern because of the method used for its assessment. Any spatial derivative of it will magnify this noise. As a result, the strain field computed from experimental images is too noisy to allow a direct comparison between the measurements obtained in the laboratory and in the simulation. We thus restrict the comparison between both approaches to the analysis of the displacement field.

Figure 6 better highlights the amplitude similarities by comparing the magnitude of the vertical displacement $u_Y(X = 0, Y)$ extracted along the central line passing through the samples. Measurements obtained from Code_Aster (red dots) are compared to the estimates obtained from laboratory measurements (black dots). The vertical displacements are averaged over a width equal to the characteristic length of the elements, $l_c = 0.5$ mm. At this scale, the experimental and simulated displacement fields agree on average, which verifies the concordance between both approaches in terms of local strain.

At the sample scale, we also compare the mechanical behavior of the samples. Colored crosses in Figure 7 show the experimental behavior of the sample by indicating the force $F$ (kN) as a function of the variation of height $\delta$ (µm) measured during three distinctive experiments. We emphasize that the same sample is subjected three times to the same multistage load. The repetitive experiments intend to test the reversibility and the repeatability of the behavior highlighted from experimental measurements. The latter are fitted with a linear regression with a determination coefficient $r^2$ of 0.96 and with a slope of $1.31 \times 10^9$ N m$^{-1}$. The determination coefficient result suggests that the observed measurements are well replicated by a linear relationship between force and displacement. The mechanical behavior simulated using Code_Aster is shown as a black line (Figure 7). The force applied to the sample is calculated every 10 µm (black dots in Figure 7). Since only an elastic constitutive law is introduced in the simulation, the force varies linearly with the displacement applied to the sample. The slope is $1.29 \times 10^9$ N m$^{-1}$, very close to the experimental behavior with only a 1.3% difference in the measured slopes. The very similar linear behavior between laboratory and numerical measurements shows the good agreement between both approaches and evidences a macroscopic linear elastic deformation of both samples.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The vertical displacements $u_Y (X = 0, Y)$ measured by Code_Aster (red dots) along a vertical line centered at the middle of the sample ($X = 0$) are compared to the laboratory measurements obtained from the image correlation procedure. Laboratory and simulated measurements are obtained for an imposed displacement of $\delta = 60$ µm, which corresponds to a loading force of $F = 75$ kN.</figcaption></figure>

The superposition of the measurements for repetitive loading highlights that no residual deformation exists. It demonstrates the reversibility of the deformation of the laboratory sample and shows that there is no inelastic deformation, like damage or plasticity. The multiscale comparison of the displacement measurements, that is, the similar spatial distribution of the local displacement field (Figures 5 and 6) and the common linear macroscopic behavior (Figure 7), is essential to verify that the numerical approach reproduces well the mechanical behavior of the experimental sample in the limit of the resolution of laboratory measurements. This complements the similar wave scattering properties, as shown in the previous section.

We finally checked that the medium behaves elastically in our numerical simulations, at the scale of the finite elements and in the range of displacement explored despite possible stress concentrations in the vicinity of the holes. The verification is done by computing Von Mises's stresses along the meshes. To

be in the elastic regime, Von Mises's stresses have to be smaller than the elastic limit in compression, *Re*. If not, the material is expected to enter the plasticity domain. Von Mises's criterion has been shown to be well suited for metal like Au4G (Ford & Alexander, 1963). In our simulations, the criterion is verified in the bulk of the sample, for the range of displacement explored in the laboratory. When the displacement is lower than $\delta = 80$ µm, the ratio of Von Mises's stress and the elastic limit *Re* is below 0.75. The criterion is nevertheless exceeded on two points of the top boundary, more precisely for the elements in contact with both ends of the loading line (see Figure S2 in the supporting information).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Force (kN) applied by the piston on the sample, as a function of the variation of height (µm) measured by image correlation (colored crosses) during three distinctive loadings in the laboratory. The laboratory measurements are compared to the simulated measurements from Code_Aster (black dots) which are fitted by a linear regression (plain black line).</figcaption></figure>

## 4. CWI Monitoring

The similarity in mechanical and acoustic behavior highlighted previously between our laboratory experiment and numerical simulations legitimizes the comparison of laboratory and simulated CWI measurements. CWI is applied to the simulated and experimental waveforms, to assess time delays during the first 0.1 ms. As presented in section 3, the waveform is weakly influenced by the intrinsic absorption of the material during that period (see Figure 4 and related text in section 3.1). In the following paragraphs, we first present the laboratory CWI measurements and then detail the results of the numerical model.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Waveforms recorded during the laboratory experiments filtered in the (300–600) kHz range. The “reference” waveform (black line) is superimposed to the “perturbed” waveform (red line) (a) recorded when applying a displacement of 50 µm (which corresponds to an applied force F = 60 kN). Zooms in two 6.5 µs long windows at the beginning (b) and at the end (c) of the coda.</figcaption></figure>

## 4.1. Laboratory Measurements

Figure 8 shows an example laboratory record on which the CWI technique is applied. The waveform is recorded when applying a 50 µm displacement to the sample and compared to the “reference” waveform for no loading. Both waveforms are band-pass filtered from 300 to 600 kHz. We first see that no direct wave arrivals are clearly visible. This result supports that the coda is in the multiple-scattering regime, in the frequency range chosen. In Figures 8b and 8c, we focus on two distinctive time periods within the signal by zooming into two 6.5 µs long widows at the beginning and at the end of the coda. The zoom for early times (Figure 8b) shows that the reference and perturbed waveforms are similar and superimpose well. The relative wave paths are expected to be quasi-direct with little influence from the deformation. At later times (Figure 8c), travel time delays clearly appear between the “reference” and “perturbed” waveforms. We interpret this difference in behavior between the early part of the waveform and the last part, that is, the coda part, as the larger sensitivity of multi-ply scattered waves to slight changes in the medium compared to direct arrivals (Snieder, 2006). With our uniaxial loading, the “perturbed” waveform appears compressed in time. Applying CWI techniques intends to quantify this time compression and to estimate the relative variation of delays with time, *dt/t*, that quantifies the perturbation in the waveforms.

Sens-Schönfelder and Wegler (2006) propose to estimate the relative time delays, *dt/t*, by searching the factor by which the time axis of one of the traces has to be elongated (apparent decrease in velocity) or compressed (apparent increase in velocity) to obtain the best correlation with the other trace, that is, the stretching method. The “stretched” version $S'(t)$ of the reference signal $S(t)$ is then obtained as $S'(t) = S[(1 + a)t]$ where a is the stretching coefficient. The coefficient a is obtained using the following procedure. All waveforms are normalized, and the laboratory waveforms are interpolated using a piecewise cubic spline interpolation in order to adjust the sampling to the finer sampling rate used in the simulation. We use here a stretching of the whole 0.1 ms long waveform. Before application of the stretching technique, a passband filter is applied to the simulated and laboratory waveforms from 300 to 600 kHz (see section 2.3.2 for details about the frequency band selection). In the laboratory experiments and in the simulation, waveforms are recorded on 32 distinctive receivers equally distributed along the sample side. The method is applied on “reference” and “perturbed” waveforms recorded at each receiver in order to analyze the spatial variability of our measurements. We thus discuss the strain sensitivity of the relative time delay *dt/t* averaged on the 32 estimates obtained at each position of the receiver.

Figure 9 shows the relative time shift measured from the “stretching” technique in the laboratory experiments (colored crosses) as a function of the displacement applied at the top of the sample. Figure 9 depicts a clear dependency of both variables. It can be fitted by a linear trend with a slope of $-2.6 \times 10^{-5} \mu\text{m}^{-1}$ (black line in Figure 9) for which we calculated a sum of squared residuals of $1.7 \times 10^{-6} \mu\text{m}$.
<figure><img src="image_12.png" /></figure>

**Figure 9.** Comparison between the relative time shifts *dt/t* of three distinctive laboratory experiments (colored crosses) and the simulated ones as a function of the displacement applied to the sample. In the numerical approach, we model the influences on *dt/t* of the change of shape of the medium and of wave velocity variations due to density variations and due to acoustoelastic effects using third-order elastic coefficient (TOEC) of Dural (red dotted line). Relative time shifts are measured by the stretching technique and are the mean values calculated from the 32 distinctive estimates obtained along the receiver array.

The CWI measurements represented in Figure 9 are those obtained during three successive and independent experiments. During each experiment, the sample is subjected to an identical loading path. The reproducibility of the laboratory CWI measurements proves that the mechanisms affecting the CWI delay measurements are reversible, which is in accordance with the reproducibility of the macroscopic loading curve (Figure 7), characteristic of an elastic deformation of the sample.

## 4.2. Numerical Simulations

In the laboratory, the scattered wavefield carries the signature of the different phenomena induced by the elastic deformation of the propagation medium and responsible for travel time variations. The numerical modeling enables to go beyond the limits of the laboratory experiments and to test the relative contributions of the different effects that are expected to influence the wave scattering, which opens perspectives to improve our understanding of the measurements. Snieder et al. (2002) interpret relative time delay changes as small variations of intrinsic seismic velocities, due to factors such as water saturation or air pressure variations, change in temperature changes, or change in the stress state (Grêt et al., 2006; Lockner et al., 1977; Planès & Larose, 2013; Snieder et al., 2002). According to this interpretation, time delays are supposed to be related to an intrinsic velocity change, such that $dt/t = -dv/v$. Effects of deformation on the wave scattering are accordingly neglected, being assumed to be a few orders of magnitude lower than the effects of intrinsic velocity changes. Conversely, the identification of nonnegligible reversible thermomechanical deformation contributions in CWI measurements has been the object of recent studies (Azzola, Griffiths, et al., 2018; Azzola, Schmittbuhl, et al., 2018). The numerical model proposed in the present work follows these studies and aims at assessing from the comparison of two approaches (analog and laboratory approaches), the relative contribution of the deformation of the sample and of the wave velocity variations in the CWI measurements which is typically a difficult task from a single approach. The three following subsections describe the three physical phenomena included progressively into the simulation and present the related measurements.

### 4.2.1. Numerical Modeling of the Effect of Sample Boundary Deformation

The uniaxial loading of the sample deforms the medium and its scatterers. We model here its influence on the scattered wavefield. To do so, we limit the update of the model, at each step of the loading, to the deformation of the mesh grid using *Code\_Aster* without changing any property of the model. We then propagate the wavefield in the deformed mesh grid using SPECFEM2D. Keeping the density, Young's modulus and Poisson ratio of the sample constant, we can test the sole influence of the change in shape of the medium on the CWI results without imposing any wave velocity variations. In doing so, we focus on the impact of the elastic deformation of the boundaries of the sample and of its scatterers, which is a contribution usually neglected when following a classical interpretation of the CWI measurements in terms of wave velocity variations. As shown in Figure 10 (red solid line—Model 1), the simulation predicts in this case, as for the experiments, relative time delays *dt/t* that linearly evolve with the stress $\sigma_Y = F/S$ where *F* is the force applied to the sample and *S* the contact surface between the sample and the piston. However, their amplitude in Model 1 is significantly smaller than the laboratory measurements. The linear regression fitted to the simulated CWI measurements shows a slope of $-3.5 \times 10^{-6} \text{ MPa}^{-1}$. The contribution of the changes in shape of medium and its included scatterers reproduces a similar linear behavior but only explains 11% of the CWI signal amplitude measured in the experiment.

### 4.2.2. Numerical Modeling of the Influence of Density Variations

We introduce here wave velocity variations of the sample that only result from the local changes of the sample density $\ell$ following equations (6) and (7) and keeping the other elastic moduli ($E$, $\nu$) constant. Using these equations that are specific to the elastic deformation regime, we show that density effects contribute to CWI measurements by weakly reducing the time delays (see Figure 10—Model 2). Assuming small deformation, volumetric strain $\epsilon_{vol}$ equals the trace of the strain tensor Tr($\epsilon$) and the relative variation of volume $dVol$.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison between the linear regression fitted to the time shifts $dt/t$ of the laboratory experiments (black line) and the ones obtained from simulations based on three distinctive models as a function of the stress $\sigma_Y$. Model 1 corresponds to the relative contribution of the change of sample shape including the boundaries of the holes (i.e., the scatterers) (red plain line). Model 2 is an extended version of Model 1 with also the wave velocity variations due to density variations (red dashed line). The red dotted line (Model 3) depicts the measurements when modeling all the influences on $dt/t$, that is, the impact of the change of sample shape, the change in the wave velocity variations due to density variations, and those due to acoustoelastic effects using third-order elastic coefficient (TOEC) of Duraluminium. Relative time shifts are measured by the stretching technique and are the mean values calculated from the 32 distinctive estimates obtained along the receiver array.</figcaption></figure>

$Vol_0$. Assuming a constant mass, the density is related to its initial value $\rho_0$ and the volumetric strain as follows:

$$ \rho = \frac{\rho_0}{\varepsilon_{vol} + 1}, \qquad (5) $$

$$ v_P = \sqrt{\frac{\nu E + E (1 - 2 \nu)}{\rho (1 + \nu) (1 - 2 \nu)}}, \qquad (6) $$

$$ v_S = \sqrt{\frac{E}{2 \rho (1 + \nu)}}. \qquad (7) $$

The density is calculated in SPECFEM2D for each spectral element according to equation (5) using the trace of the strain tensor computed by Code_Aster at a given step of the loading procedure. Wave velocities are derived using equations (6) and (7), keeping the elastic moduli $E$ and $\nu$ constant (i.e., assuming that the sample is loaded without damage). Because wave velocities are related to the square root of the inverse of density (see equations (6) and (7)), one expects that relative time delays decrease owing to the effect of density variations and of the medium compression. Indeed, the contribution of wave velocity variations due to density changes from volumetric strain is of opposite sign to those of the other modeled phenomena. As a result, adding the effect of wave velocity changes due to local density changes to the contribution of the change in shape of the sample leads to a decrease in the relative delays $dt/t$ (red dashed line in Figure 10). The CWI measurements ($dt/t$) measured using Model 2, which includes density effects and geometrical deformation of the sample, explain only 7% of the CWI signal measured in the experiments.

### 4.2.3. Numerical Modeling of the Influence of the Acoustoelastic Effect

Finally, we explore the impact of acoustoelastic effect on elastic wave velocities, related to nonlinear elasticity (Murnaghan, 1951). Acoustoelastic effects refer to the changes in the elastic wave velocities related to stress-induced effects. It has been theoretically explained for a Green elastic or hyperelastic medium, that is, an elastic medium where stresses depend only on the current state of deformation (Smith, 1963). In this case, the stored elastic energy per unit of undeformed volume $W_{elas}$ is a polynomial function of the strain tensor $\varepsilon$:

$$ W_{elas}(\varepsilon) = \frac{1}{2} C^{(2)}_{ijkl} \varepsilon_{ij} \varepsilon_{kl} + \frac{1}{6} C^{(3)}_{ijklmn} \varepsilon_{ij} \varepsilon_{kl} \varepsilon_{mn} + \dots, \qquad (8) $$

where the coefficients in $C^{(2)}_{ijkl}$ are the second-order elastic constants and in $C^{(3)}_{ijklmn}$ the third-order elastic constants (TOEC) (Smith, 1963; Toupin & Bernstein, 1961). Classically, for small deformation, only the second-order elastic constants are considered. They describe the first term in the expansion of the elastic energy and refer to the first-order linear elasticity in the expansion of the stresses as a function of the strains. When a higher expansion of the stored elastic energy is considered, TOEC that correspond to a quadratic expansion of the stresses in terms of strain are relevant. Hughes and Kelly (1953) showed in this case that, even for an isotropic elastic material, the velocities of the compression and shear waves can differ by a small amount which is proportional to the load. Subsequently, the acoustoelastic theory predicts an anisotropic effect of the uniaxial elastic deformation of the propagation medium on the elastic wave velocities. In particular, the anisotropic behavior highlighted is due to nonlinearities in the strain-displacement relationship. We introduced in our model the constitutive relations of acoustoelasticity that have been confirmed for metallic materials (Hughes & Kelly, 1953). They are derived from Murnaghan's theory of finite deformations (Murnaghan, 1951) using TOEC. These coefficients (TOEC) can be determined experimentally by measuring the change in acoustic phase velocity as a function of stress (Si-chaib et al., 2011; Smith, 1963). Here, we used
Murnaghan's TOEC L, M, and N. These parameters have been measured for Duraluminium by Chatellier (1987) and by Si-chaib et al. (2011), which employ these constants in simulations. The three TOEC are respectively L = -51 ± 5 GPa, M = -325 ± 2 GPa, and N = -351 ± 6 GPa. Second-order elastic constants Λ and μ (or Lamé constants) are calculated using Young's modulus E and the Poisson ratio ν of the Au4G sample as follows:

$$ \Lambda = \frac{E \nu}{(1 + \nu)(1 - 2 \nu)}, \qquad (9) $$

$$ \mu = \frac{E}{2(1 + \nu)}, \qquad (10) $$

We obtained respectively Λ = 54 GPa and μ = 28 GPa. In order to model the acoustoelastic effects in SPECFEM2D, we defined, for each spectral element of the mesh, pseudo second-order elastic coefficients C*_{ijkl}. These coefficients are calculated from the second-order constants and TOEC and using the principal components ε₁, ε₂, and ε₃ of the strain tensor computed at each step of the loading procedure from Code_Aster, according to Hughes and Kelly (1953):

$$ \sigma_{ij} = \sigma^{0}_{ij} + \sum_{kl} C^{*}_{ijkl} \varepsilon_{kl}, \qquad (11) $$

where σ_{ij} is the current stress tensor, σ^{0}_{ij} is the initial stress tensor, and

$$ C^{*}_{ijkl} = [\Lambda + 2[L - \Lambda - M]\varepsilon_{vol} + 2[\Lambda + M]\varepsilon_{l} + \varepsilon_{k} - 2\mu\varepsilon_{j}]\delta_{i,j}\delta_{k,l} \\ + [\mu + (\Lambda + M - \mu)\varepsilon_{vol} + 2\mu(\varepsilon_{i} + \varepsilon_{j} + \varepsilon_{l})(\delta_{i,k}\delta_{j,l} + \delta_{i,l}\delta_{j,k})] + \frac{1}{2}N\sum_{v}(\delta_{j,v,k}^{i,v,l} + \delta_{j,v,l}^{i,v,k})\varepsilon_{vol}, \qquad (12) $$

with ε_{vol} the trace of the strain tensor and where Λ and μ are ordinary Lamé constants and L, M, and N the third-order Murnaghan constants. In equation (12), the δs are the general Kronecker delta functions. Acoustoelastic effects on scattered waves have been predicted in the case of the interaction of two monochromatic plane waves by Jones and Kobett (1963)), but the influence of the sample stressing on the CWI measurements has, to our knowledge, not been explored yet. The relative time delays dt/t measured for different loading states exhibit a significant impact of the acoustoelastic effects (Figures 9 and 10). The curve when including these effects (red dotted line—Model 3) shows a clear increase of the relative delays with the loading and an interesting superposition to the experimental results. The CWI measurements obtained with Model 3 overestimate by 8% the laboratory observations. This superposition, together with the strong discrepancy highlighted in Figure 10 between Models 1 and 3, shows that acoustoelastic effects have a major impact on the CWI measurements. The strength of the related signal suggests that monitoring nonlinear elastic effects could be carried out based on the analysis of strongly scattered waves.

# 5. Discussion

## 5.1. Uncertainties Related to TOEC Measurements

Considering the preponderant contribution to the CWI measurements of the acoustoelastic effects, we discuss here the sensitivity of the results of Model 3 (Figure 10) toward the parameters involved in the computation of the pseudo second-order elastic coefficients C*_{ijkl}, that is, the Lamé second-order elastic parameters and the Murnaghan third-order elastic parameters. We first tested the impact of increasing or decreasing all three Murnaghan parameters concomitantly by 10% in the simulations. This value of 10% is a large proxy of the uncertainty on the measurements of the third-order elastic parameters. Indeed, the third-order elastic parameters proposed for Au4G by Chatellier (1987) are measured with an uncertainty of 10% on L, of 0.6% GPa on M, and of 1.7% GPa on N. Following the change in parameters, we observed an 8% variation of the measured relative time delays. Second, we impose a variation of 10% of the second-order parameters, which are calculated using the density, Young's modulus, and the Poisson ratio. This induces an 11% variation of the measured relative time delays dt/t.

However, the concomitant change of all three TOECs does not depict the influence of each parameter on the elastic wave velocities and thus on the CWI measurements. To gain more insights in the influence of these
elastic constants on the CWI measurements, we cross tested and applied a 10% decrease—or increase—in the value of a single, of two, or of all three coefficients concomitantly. In each case, we quantified the variation from the relative travel time changes observed with Model 3 and using the initial set of parameters. The results of the sensitivity tests are tabulated in the supporting information (Text S4), and we describe here the two major observations.

First, the increase—respectively decrease—of each third-order coefficient causes an increase—respectively decrease—in the relative time delays measured in the simulation. Then, when testing the impact of a multivariate change in the TOECs, we show that the effects add up linearly. For example, we measure a 3.2% increase in relative time delays when we apply a variation of 10% only to the value of the parameter *M*. This variation is respectively 4.1% when applying the change to the value of the coefficient *N*. When we test the impact of a concomitant increase in *M* and *N*, we consequently measure a 7.3% variation in the relative time delays, that is, the sum of the variations measured when increasing a single parameter.

This sensitivity test shows that the laboratory CWI measurements are in the range of the simulated results when considering an uncertainty of 10% on the three TOECs, that is, *L*, *M*, and *N*.

In order to catch closely the physics expressed in the laboratory experiments, the simulations necessitate a finer tuning of these mechanical parameters. It would need the direct measurements of the second and third elastic parameters of our holed Au4G sample, which was beyond the focus of this study.

## 5.2. Approximation of the 2-D Modeling

An important simplification of our numerical modeling is the use of a 2-D representation of the 3-D experimental sample. The advantages of the 2-D numerical models are principally linked to computational performances. They are also related to the significant simplification of the discretization of numerical samples, in particular when they involve a complex assembly of scatterers. We justify in the following section the choice of the 2-D approximation and treat the arguments in favor of the bidimensional modelization of both the elastic deformation of the sample and of the wave propagation.

The choice of the plane deformation modeling of the experiments in 2-D is motivated by the invariance of the geometry of the sample in the *Z* direction, with cylindrical holes drilled through the 32 mm thick block. In a typical plane strain problem, the loading is such that the deformation is essentially exerted in the cutting plane (*X*, *Y*), which is a relevant assumption in our case of a step-by-step loading imposed from the top of the sample. To verify the pertinence of this approximation, we modeled the 3-D mechanical response of the holed Au4G sample with a 32 mm thickness. Using *Code_Aster* for an applied displacement (δ = 75 µm), we compared the map of the trace of the local strain tensor, computed in the middle plane (*Z* = 0) of the 3-D sample to that of the 2-D sample (see Figure 11). It shows a similar pattern of the local volumetric strain with a more pronounced volumetric deformation in the 2-D case. We quantified locally the difference highlighted in Figure 11 by averaging the trace of the strain tensor along a vertical line at the middle (*X* = 0) of the sample presented in Figure 11. The mean values are 1.15 × 10⁻⁴ in the 2-D case and 6.34 × 10⁻⁵ in the 3-D case. The difference in the magnitude of the trace of the tensor components, which is of 45%, shows the contribution of the orthogonal component of the strain tensor ε<sub>ZZ</sub>, which is different from zero in the three-dimensional case. However, the local measurement obtained along the vertical line of Figure 11 does not depict the distribution of the magnitude of the strain components throughout the entire sample. We averaged the amplitude of the tensor components computed by *Code_Aster* over the entire 3-D medium. Our measurements show that the orthogonal component ε<sub>ZZ</sub> represents 20% in magnitude of the predominant component ε<sub>YY</sub>. If the magnitude of the orthogonal component ε<sub>ZZ</sub> of the strain tensor is not negligible given the magnitude of the strain components exerted in the cutting plane (*X*, *Y*) of the 3-D medium, the ratio between ε<sub>ZZ</sub> and ε<sub>YY</sub> shows that the latest is predominant, which legitimizes our simplification of the three-dimensional mechanical problem into a 2-D problem in plane strain.

In the bidimensional modelization of the wave propagation, we neglect the effects on the wave scattering in the *Z* direction. This approximation is mainly due to geometric considerations. The cylindrical holes are drilled through the entire thickness of the laboratory sample, and they are significantly longer than the wavelength. The geometry of the sample is designed such that the wave scattering is expected to occur principally in the cutting plane (*X*, *Y*) (Derode et al., 2001; Tourin et al., 2000).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison between the trace of the strain tensor computed in *Code_Aster* within the 2-D sample and within the middle plane (Z = 0) in the 3-D simulation of the elastic deformation of the holed Au4G sample. Black line represents the median line in the 2-D and 3-D models along which the trace of the strain tensor is averaged to quantify the difference in behaviors.</figcaption></figure>

**5.3. Application of the Numerical Model**

Under the assumption of the stability of the sources of the noise, the CWI of stacked Ambient Noise Cross-Correlation Functions allows typically to follow the fine-scale evolution of the medium, mostly visible in the coda part (e.g., Brenguier, Shapiro, et al., 2008; Brenguier, Campillo, et al., 2008; Campillo & Roux, 2015; Sens-Schönfelder et al., 2014; Sens-Schönfelder & Wegler, 2011). However, the interpretation of the temporal changes highlighted in different frequency bands remains uncertain, as the mechanisms at their origin are not clearly identified (Sens-Schönfelder & Wegler, 2006; Yates et al., 2019).

CWI and in particular Ambient Noise Interferometry are highly sensitive monitoring techniques with a very important potential of applications but are presently lacking precise physical interpretations of the measured signals. Our approach opens perspectives for new interpretations of these CWI measurements. Indeed, we identify a clear contribution of the volumetric elastic deformation. We decipher the dominating effect on wave scattering during the reversible elastic deformation of the sample: acoustoelastic effects.

The results provided by our analysis have implications for the monitoring of natural geological structures, such as a fractured reservoir or porous media. The upscaling of the proposed numerical model would provide new outcomes for the monitoring of geological reservoirs or large-scale crustal faults. Our results might therefore provide new perspectives to study the application of Ambient Noise Interferometry to monitor the reservoir evolution in a variety of operating environments, such as injection and extraction of fluid (Hillers et al., 2015).

**6. Conclusion**

We monitor the wave scattering changes within an elastic sample that includes a set of scatterers, during a uniaxial compression, using a CWI method. We compare results from two approaches. One is based on a laboratory experiment where a reversible elastic deformation is applied to a homogeneous Au4G sample that includes a set of cylindrical holes designed to favor multiple scattering. The other is a forward numerical modeling of the same sample with the same geometry, set of scatterers, mechanical behavior, and under the same loading conditions. The numerical model allows to focus on the relative signatures (a) of the
sample deformation without intrinsic wave velocity variations, (b) of the local density changes due to local volumetric strain, and (c) of acoustoelastic effects due to nonlinear elasticity.

Our approach shows the dominant influence of nonlinear acoustoelastic effects on the scattered wavefield and subsequently on the CWI measurements. When modeling CWI measurements without these acoustoelastic effects, that is, with only the impact of the deformation of the sample boundaries and of scatterers and/or of the density change related to the volumetric strain, only 7% to 11% of the relative time delays measured in the laboratory is numerically reproduced. We note however that besides the dominance of the acoustoelastic effects, the effect of shape changes of the medium, which is a strain effect without any velocity changes, is not negligible since it explains about 11% of the measured relative time delays. This result highlights that the strain effects on the CWI measurements are not negligible even if they remain small compared to the acoustoelastic effects. The identification of different contributions to the travel time delays quantified by CWI is the basis for future theoretical developments.

By comparing the signature on wave scattering of different physical phenomena involved in the reversible elastic deformation of the sample, our approach opens new perspectives for the interpretation of CWI measurements that could have implications for the monitoring infrastructures or natural geological objects (e.g., fractured reservoir and porous medium).

## Availability of Data and Materials

Data displayed or used to generate figures and plots are accessible online (DOI: 10.6084/m9.fishshare.7485158).

## Competing Interests

The authors declare no competing financial interest.

## Funding

This work has been published under the framework of the LABEX ANR-11-LABX-0050-G-EAU-THERMIE-PROFONDE and benefits from a funding from the state managed by the French National Research Agency (ANR) as part of the "Investments for the Future" program. It has also been funded by the EU projects DESTRESS (EU H2020 research and innovation program, Grant Agreement 691728).

## Acknowledgments

The authors would like to thank Alain Steyer for sample preparation and his help with the construction of the experimental setup. For the use of SPECFEM2D, a free software is downloadable online (https://github.com/geodynamics/specfem2d.git); we thank the Computational Infrastructure for Geodynamics (http://geodynamics.org) which is funded by the National Science Foundation under Awards EAR-0949446 and EAR-1550901. For the use of Code_Aster, an open-access software is downloadable online (https://www.code-aster.org/); we thank Electricité de France (EDF). The manuscript benefited from constructive comments by Kasper van Wijk, an anonymous reviewer, Associate Editor Nori Nakata, and Editor Yehuda Ben-Zion.

## References

Aki, K., & Richards, P. G. (2002), Quantitative seismology.

Anugonda, P., Wiehn, J. S., & Turner, J. A. (2001). Diffusion of ultrasound in concrete. *Ultrasonics*, 39(6), 429-435. https://doi.org/10.1016/S0041-624X(01)00077-4

Azzola, J., Griffiths, L., Schmittbuhl, J., Zigone, D., Magnenet, V., Masson, F., et al. (2018). Coda wave interferometry during the heating of deep geothermal reservoir rocks. *Geothermal Energy*, 6(1), 21. https://doi.org/10.1186/s40517-018-0107-2

Azzola, J., Schmittbuhl, J., Zigone, D., Magnenet, V., & Masson, F. (2018). Direct modeling of the mechanical strain influence on coda wave interferometry. *Journal of Geophysical Research: Solid Earth*, 123, 3160-3177. https://doi.org/10.1002/2017JB015162

Boué, P., Roux, P., Campillo, M., & Briand, X. (2014). Phase velocity tomography of surface waves using ambient noise cross correlation and array processing. *Journal of Geophysical Research: Solid Earth*, 119, 519-529. https://doi.org/10.1002/2013JB010446

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N. M., Nadeau, R. M., & Larose, E. (2008). Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations. *Science*, 321(5895), 1478-1481. https://doi.org/10.1126/science.1160943

Brenguier, F., Shapiro, N. M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., & Nercessian, A. (2008). Towards forecasting volcanic eruptions using seismic noise. *Nature Geoscience*, 1(2), 126-130. https://doi.org/10.1038/ngeo104

Campillo, M., & Roux, P. (2015). Crust and lithospheric structure—Seismic imaging and monitoring with ambient noise correlations. In *Treatise on Geophysics*, (pp. 391-417). Amsterdam: Elsevier. https://doi.org/10.1016/B978-0-444-53802-4.00024-5

Chatellier, J.-Y. (1987), Simulation acousto-élastoplastique du formage des métaux : Application à la détermination des contraintes résiduelles (PhD Thesis).

Dainty, A. M., & Toksoz, M. N. (1977). Elastic wave propagation in a highly scattering medium—A diffusion approach. *Journal of Geophysics Zeitschrift Geophysik*, 43, 375-388.

Dainty, A. M., & Toksoz, M. N. (1981). Seismic codas on the Earth and the Moon: A comparison. *Physics of the Earth and Planetary Interiors*, 26(4), 250-260. https://doi.org/10.1016/0031-9201(81)90029-7

Derode, A., Tourin, A., & Fink, M. (2001). Random multiple scattering of ultrasound. I. Coherent and ballistic waves. *Physical Review E*, 64(3). https://doi.org/10.1103/PhysRevE.64.036605

Ford, H., & Alexander, J. M. (1963), Advanced mechanics of materials.

Grêt, A., Snieder, R., & Scales, J. (2006). Time-lapse monitoring of rock properties with coda wave interferometry: Time-lapse monitoring of rock properties. *Journal of Geophysical Research*, 111, B03305. https://doi.org/10.1029/2004JB003354
Hillers, G., Husen, S., Obermann, A., Planès, T., Larose, E., & Campillo, M. (2015). Noise-based monitoring and imaging of aseismic transient deformation induced by the 2006 Basel reservoir stimulation. *Geophysics, 80*(4), KS51-KS68. https://doi.org/10.1190/geo2014-0455.1

Hughes, D. S., & Kelly, J. L. (1953). Second-order elastic deformation of solids. *Physical Review, 92*(5), 1145–1149. https://doi.org/10.1103/PhysRev.92.1145

Jones, G. L., & Kobett, D. R. (1963). Interaction of elastic waves in an isotropic solid. *The Journal of the Acoustical Society of America, 35*(1), 5–10. https://doi.org/10.1121/1.1918405

Komatitsch, D., & Vilotte, J.-P. (1998). The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. *Bulletin of the Seismological Society of America, 88*(2), 368–392.

Kuperman, W. A., Montagner, J.-P., & Tourin, A. (2002). *Imaging of complex media with acoustic and seismic waves*, (Vol. 82). Berlin: Springer Science & Business Media.

Larose, E., & Hall, S. (2009). Monitoring stress related velocity variation in concrete with a 2 × 10⁻⁵ relative resolution using diffuse ultrasound. *The Journal of the Acoustical Society of America, 125*(4), 1853–1856. https://doi.org/10.1121/1.3079771

Lin, F.-C., Li, D., Clayton, R. W., & Hollis, D. (2013). High-resolution 3D shallow crustal structure in Long Beach, California: Application of ambient noise tomography on a dense seismic array. *Geophysics, 78*(4), Q45-Q56. https://doi.org/10.1190/geco2012-0453.1

Lockner, D. A., Walsh, J. B., & Byerlee, J. D. (1977). Changes in seismic velocity and attenuation during deformation of granite. *Journal of Geophysical Research, 82*(33), 5374–5378. https://doi.org/10.1029/JB082i033p05374

Margerin, L. (2005). In A. Levander, & G. Nolet (Eds.), *Introduction to radiative transfer of seismic waves, Geophysical Monograph Series*, (Vol. 157, pp. 229–252). Washington, D. C: American Geophysical Union. https://doi.org/10.1029/157GM14

Meunier, J., Huguet, F., & Meynier, P. (2001). Reservoir monitoring using permanent sources and vertical receiver antennae: The Céré-la-Ronde case study. *The Leading Edge, 20*(6), 622–629. https://doi.org/10.1190/1.1439008

Moschetti, M. P., Ritzwoller, M. H., & Shapiro, N. M. (2007). Surface wave tomography of the western United States from ambient seismic noise: Rayleigh wave group velocity maps: U.S. Surface Wave Tomography. *Geochemistry, Geophysics, Geosystems, 8*, Q08010. https://doi.org/10.1029/2007GC001655

Murnaghan, F. D. (1951). *Finite deformation of an elastic solid*. New York: John Wiley & Sons.

Mussett, A. E., & Khan, M. A. (2000). *Looking into the Earth: An introduction to geological geophysics*. Cambridge: Cambridge University Press.

Olivier, G., Brenguier, F., Campillo, M., Lynch, R., & Roux, P. (2015). Body-wave reconstruction from ambient seismic noise correlations in an underground mine. *Geophysics, 80*(3), KS11-KS25. https://doi.org/10.1190/geo2014-0299.1

Planès, T., & Larose, E. (2013). A review of ultrasonic Coda Wave Interferometry in concrete. *Cement and Concrete Research, 53*, 248–255. https://doi.org/10.1016/j.cemconres.2013.07.009

Planès, T., Larose, E., Margerin, L., Rossetto, V., & Sens-Schönfelder, C. (2014). Decorrelation and phase-shift of coda waves induced by local changes: Multiple scattering approach and numerical validation. *Waves in Random and Complex Media, 24*(2), 99–125. https://doi.org/10.1080/17455030.2014.880821

Poupinet, G., Ellsworth, W. L., & Frechet, J. (1984). Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California. *Journal of Geophysical Research, 89*(B7), 5719–5731. https://doi.org/10.1029/JB089iB07p05719

Raffel, M., Willert, C. E., Scarano, F., Kähler, C. J., Wereley, S. T., & Kompenhans, J. (2018). *Particle image velocimetry: A practical guide*. Berlin: Springer.

Ratdomopurbo, A., & Poupinet, G. (1995). Monitoring a temporal change of seismic velocity in a volcano: Application to the 1992 eruption of Mt. Merapi (Indonesia). *Geophysical Research Letters, 22*(7), 775–778. https://doi.org/10.1029/95GL00302

Sato, H., Fehler, M. C., & Maeda, T. (2012). *Seismic wave propagation and scattering in the heterogeneous Earth*, (Vol. 496). Berlin: Springer.

Sens-Schönfelder, C., Pomponi, E., & Peltier, A. (2014). Dynamics of Piton de la Fournaise volcano observed by passive image interferometry with multiple references. *Journal of Volcanology and Geothermal Research, 276*, 32–45. https://doi.org/10.1016/j.jvolgeores.2014.02.012

Sens-Schönfelder, C., & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. *Geophysical Research Letters, 33*, L21302. https://doi.org/10.1029/2006GL027797

Sens-Schönfelder, C., & Wegler, U. (2011). Passive image interferometry for monitoring crustal changes with ambient seismic noise. *Comptes Rendus Geoscience, 343*(8-9), 639–651. https://doi.org/10.1016/j.crte.2011.02.005

Shapiro, N. M., Shapiro, N., Campillo, M., Stehly, L., & Ritzwoller, M. H. (2005). High-resolution surface-wave tomography from ambient seismic noise. *Science, 307*(5715), 1615–1618. https://doi.org/10.1126/science.1108339

Si-chaib, M. O., Djelouah, H., & Nour, A. (2011). Study of shear ultrasonic waves propagating in materials under static uniaxial loadings: Simulation, shear waves, acoustoelasticity, mechanical stresses, materials. *Strain, 47*, 134–143. https://doi.org/10.1111/j.1475-1305.2010.00749.x

Smith, R. T. (1963). Stress-induced anisotropy in solids—The acousto-elastic effect. *Ultrasonics, 1*(3), 135–147. https://doi.org/10.1016/0041-624X(63)90003-9

Snieder, R. (2006). The theory of Coda Wave Interferometry. *Pure and Applied Geophysics, 163*(2-3), 455–473. https://doi.org/10.1007/s00024-005-0026-6

Snieder, R., Grét, A., Douma, H., & Scales, J. (2002). Coda Wave Interferometry for estimating nonlinear behavior in seismic velocity. *Science, 295*(5563), 2253–2255. https://doi.org/10.1126/science.1070015

Thielicke, W. (2014). The flapping flight of birds: Analysis and application. [S.n.].

Thielicke, W., & Stamhuis, E. J. (2014). PIVlab—Towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB. *Journal of Open Research Software, 2*. https://doi.org/10.5334/jors.bl

Toupin, R. A., & Bernstein, B. (1961). Sound waves in deformed perfectly elastic materials. Acoustoelastic effect. *The Journal of the Acoustical Society of America, 33*(2), 216–225. https://doi.org/10.1121/1.1908623

Tourin, A., Derode, A., Peyre, A., & Fink, M. (2000). Transport parameters for an ultrasonic pulsed wave propagating in a multiple scattering medium. *The Journal of the Acoustical Society of America, 108*(2), 503–512. https://doi.org/10.1121/1.429580

Trégourès, N. P., & van Tiggelen, B. A. (2002). Quasi-two-dimensional transfer of elastic waves. *Physical Review E, 66*(3). https://doi.org/10.1103/PhysRevE.66.036601

Wegler, U., & Lühr, B. G. (2001). Scattering behaviour at Merapi volcano (Java) revealed from an active seismic experiment. *Geophysical Journal International, 145*(3), 579–592.

Yates, A. S., Savage, M. K., Jolly, A. D., Caudron, C., & Hamling, I. J. (2019). Volcanic, coseismic, and seasonal changes detected at White Island (Whakaari) Volcano, New Zealand, using seismic ambient noise. *Geophysical Research Letters, 46*, 99–108. https://doi.org/10.1029/2018GL080580
Zigone, D., Ben-Zion, Y., Campillo, M., & Roux, P. (2015). Seismic tomography of the southern california plate boundary region from noise-based Rayleigh and Love waves. *Pure and Applied Geophysics*, *172*(5), 1007–1032. https://doi.org/10.1007/s00024-014-0872-1

Zigone, D., Ben-Zion, Y., Lehujeur, M., Campillo, M., Hillers, G., & Vernon, F. L. (2019). Imaging subsurface structures in the San Jacinto fault zone with high-frequency noise recorded by dense linear arrays. *Geophysical Journal International*, *217*(2), 879–893. https://doi.org/10.1093/gji/ggz069