# Numerical and experimental study of ultrasonic seismic waves propagation and attenuation on high quality factor samples

Marine Deheuvels, Florian Faucher, Daniel Brito

► To cite this version:

Marine Deheuvels, Florian Faucher, Daniel Brito. Numerical and experimental study of ultrasonic seismic waves propagation and attenuation on high quality factor samples. Geophysical Prospecting, 2023, 10.1111/1365-2478.13465. hal-04336798

HAL Id: hal-04336798

https://hal.science/hal-04336798v1

Submitted on 13 Dec 2023

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 CC BY 4.0 - Attribution - International License
# Numerical and Experimental Study of Ultrasonic Seismic Waves Propagation and Attenuation on High Quality Factor Samples

Marine Deheuvels<fnref n="1" />,<fnref n="2" /> | Florian Faucher<fnref n="2" /> | Daniel Brito<fnref n="1" />

### Correspondence

Marine Deheuvels
Email: marine.deheuvels@inria.fr

### Funding information

E2S UPPA PhD grant

We propose an approach for measuring seismic attenuation at ultrasonic frequencies on laboratory-scale samples. We use a Gaussian filter to select a bandwidth of frequencies to identify the attenuation in a small window and, by moving the window across the frequency content of the data, we determine the frequency-dependent attenuation function. We assess the validity of the method with 3D numerical simulations of seismic wave propagation across different sample geometries, using free surface boundary conditions. We perform the simulations using viscoelastic media under various seismic attenuation models. Our numerical results indicate that we can successfully recover the representative viscoelastic attenuation parameters of the media, regardless of the sample geometry, by processing the seismic signal recorded either within the volume or at the boundaries. Due to the equipartition phenomenon, the energy of S-waves is consistently higher in seismic records than that of P-waves. Therefore, we systematically recover the attenuating properties of S-waves in the medium. We also conduct experiments of seismic wave propagation on samples of aluminum and Fontainebleau sandstone to validate our approach with real data. The quality factor of the S-wave $Q_s$ in the aluminum medium increases from 300 to 7000 between 60 kHz and 1.2 MHz. The Fontainebleau sandstone,

<footnote n="1">¹Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, TotalEnergies, LFCR, France</footnote>

<footnote n="2">²Inria MAKUTU - Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS, France</footnote>
which is more attenuating, exhibits a Qs that increases from 200 at 60 kHz to 1000 at 1.2 MHz. With our approach, we are not only able to recover the attenuation property, but also identify the frequency-dependent attenuation model of the samples. Our method allows for seismic attenuation recovery at ultrasonic frequencies in low-attenuating media.

### KEYWORDS

Attenuation, Ultrasounds, Viscoelasticity, Rock Physics, Signal processing, Time-harmonic waves

## Conflict of interest

All authors declare that they have no conflict of interest.

## Data availability

The data that support the findings of this study are available on request from the corresponding author.

# 1 | INTRODUCTION

The amplitude of a propagating seismic wave decays due to geometrical spreading, scattering, and intrinsic attenuation. Attenuation is of particular interest in seismic datasets acquired in porous media (see Dasgupta and Clark (1998); Quan and Harris (1997) for attenuation determination at a field scale) as it relates to physical properties of rocks such as porosity, density of fractures, and saturation (Pang et al., 2019). The attenuation is frequency-dependent (Müller et al., 2010; Gurevich and Pevzner, 2015; Carcione, 2007) and is usually associated with the dimensionless quality factor (Q-factor) describing the ability of a medium to propagate a wave (Knopoff and MacDonald, 1958).

Viscoelastic materials are described using two components: an elastic one (represented by the elastic modulus) that does not attenuate, and a damping one (represented by the viscous modulus) that causes attenuation. Multiple viscoelastic models exist to cover the different attenuation behaviour of materials (Biot, 1954; Golden and Graham, 1988; Carcione, 2007). For porous rocks, Biot (1956a,b) developed poroelastic approaches that describe wave propagation in an elastic porous matrix filled with fluids, and in this case the intrinsic attenuation of the medium is defined as the dissipation of the wave due to friction at the fluid/solid interface. Although poroelastic models are more accurate in modeling wave propagation in rocks, they can be approximated by viscoelastic models (Bardet, 1992; Morochnik and Bardet, 1996), which are easier to consider numerically, even though they may lead to inaccuracy near the boundaries (Geertsma and Smit, 1961). This is our choice for this work.

The frequency-dependent attenuation or Q-factor of a natural rock can be studied at the laboratory scale through different methods (see Jackson (1993) and Subramaniyan et al. (2014) for reviews) where various external conditions are applied to generate propagating waves on a wide range of frequencies, from a few Hz to MHz. At relatively low frequencies (from 10 Hz to 100 kHz), the resonant bar method - also called free oscillations method - excites eigenmode oscillations (Norris and Young, 1970) whose characteristics depend on the attenuation of the sample.
This experimental method allowed for instance to highlight changes in attenuation with humidity content in air (Pandit and King, 1979) or fluid saturation (Winkler and Nur, 1979, 1982) in Berea sandstone. In Fontainebleau Sandstone, Bourbie and Zinszner (1985) measured the attenuation depending on the frequency and water saturation up to 10 kHz; they found for example a Q-factor of 100 in dry conditions around 10 kHz. Another experimental method at low frequencies is the method of forced oscillations: it consists in imposing sinusoidal mechanical stress on a sample and measuring the associated mechanical strain (Spencer Jr., 1981). The shift between stress and strain relates to Q-factor of the samples. This technique is used up to 100 kHz for high attenuation values (Q-factor < 300) (see Tisato and Madonna (2012); Tisato et al. (2014); Pimienta et al. (2015b, 2017) for measurements on sandstones).

For studies at ultrasound frequencies (around 1 MHz), the spectral ratio technique is largely used for high at-
tenuation samples (Q-factor < 200 in general): it consists in computing a relative spectral amplitude of the signal's
study between a reference sample and the one of interest (Winkler and Plona, 1982) or between an initial and a final
state of the sample (Johnston et al., 1979; Toksöz et al., 1979). The spectral ratio technique allows the study of the
attenuation in a frequency range fixed by the frequency bandwidth of the seismic pulsed source. It is used in all kind
of materials such as natural rocks (Winkler and Plona, 1982; Johnston et al., 1979; Toksöz et al., 1979; Adam et al.,
2009), synthetic media (Bourbié and Nur, 1984), loose (non-consolidated) materials (Molyneux and Schmitt, 2000),
soils (Leong et al., 2004) or highly heterogeneous materials (Molero et al., 2010). Another approach is the reflection
method that analyzes the amplitudes of reflected waves within samples at ultrasonic frequencies as well to determine
the attenuation; this method is used for example by Best et al. (2007) and Agersborg et al. (2008) in sandstones or
carbonates at 750 kHz pulse-frequency in moderately low Q-factor materials (< 100). One difficulty in studying the
attenuation is also that the attenuation model of the material with respect to frequency is unknown and should be
recovered in addition to the physical viscous properties. In the aforementioned references, the attenuation model
is usually assumed to be known apriori; however with the method we propose, we show that we can identify the
attenuation model for the investigated frequency bandwidth.

In our study, we propose an approach to recover the frequency-dependent attenuation of low attenuation mate-
rials (attenuation β < 10 000 Np s⁻¹), such as salt (Manthei et al., 2006), volcanic rocks (Rao et al., 2002) or dry pure
sedimentary rocks (Bourbie and Zinszner, 1985) at a laboratory scale (around decimeter-sized samples). The additional
interest of recovering a representative attenuation model and parameters with respect to frequency may be for imag-
ing purposes, such as quantitative inversion in viscoelastic media, where one needs initial guess (Faucher and Scherzer,
2023). We describe in Section 2 the theoretical framework of the study in the frequency domain where a mechanical
wave is propagating and attenuated. In Section 3, we propose a method to recover the frequency-dependent viscoelas-
tic parameters using the multiple wave reflections from the free-surface boundaries at every face of the sample. After
some time, the beam spreading is no longer visible in the signals due to multiple reflections. Also, these multiple
reflections cause equipartition phenomenon, where P- and S- waves energies are partitioned in favor of the S-waves.
We then use a Gaussian filter to recover the attenuation with frequency, independently of the source wavelet because
we employ the relative amplitude decay in the seismogram. In Section 4, the method is tested on 3-D viscoelastic
numerical experiments, that is, considering synthetic data. In Section 5 we carry out laboratory experiments in an
aluminum cube sample and in a dry Fontainebleau sandstone cylinder, where we recover the attenuation properties
from the measurements of the global wavefield at ultrasound frequencies, between 60 kHz and 1.2 MHz. Finally,
we discuss the experimental data in the frame of the theoretical and numerical viscoelastic modeling in Section 6 and
give our conclusions in Section 7.
## 2 | VISCOELASTIC WAVE PROPAGATION IN ISOTROPIC MEDIA

In this section we provide the mechanical wave equations that are considered for viscoelastic media. We consider the time-harmonic formulation of the wave problems, that is, we work in the frequency domain. This approach is more convenient to handle the frequency-dependent attenuation behaviour of materials (Carcione, 2007).

### 2.1 | Time-harmonic wave propagation

With the time-harmonic formulation, the time-domain vector displacement field $\hat{u}(x, t)$ (given in [m] with $x$ the spatial position and $t$ the time in [s]) and stress tensor $\hat{\sigma}(x, t)$ (in [Pa]) are considered in the form of,

$$ \hat{u}(x, t) = u(x, \omega) \cdot e^{i\omega t}, \qquad (1a) $$

$$ \hat{\sigma}(x, t) = \sigma(x, \omega) \cdot e^{i\omega t}, \qquad (1b) $$

with the harmonic pulsation $\omega = 2\pi f$ where $f$ is the ordinary frequency given in [Hz]. The equation of motion for elastic waves is,

$$ -\rho(x)\omega^2 u(x, \omega) = \nabla \cdot \sigma(x, \omega) + f(x, \omega), \qquad (2) $$

where $f$ is an internal volume source. The constitutive law of the material can be written as

$$ \sigma(x, \omega) = M(x, \omega) : \varepsilon(x, \omega) = \frac{1}{2} M(x, \omega) : (\nabla u(x, \omega) + (\nabla u(x, \omega))^T), \qquad (3) $$

where $\varepsilon$ is the strain and $^T$ denotes the transposed. Here, the medium is described by the viscoelastic stiffness tensor $M$ (in [Pa]) and density $\rho$ (in [kg m$^{-3}$]). To handle attenuation in the frequency domain, the viscoelastic tensor $M$ is complex-valued and frequency-dependent.

### 2.2 | Viscoelastic rheologic models

There exist several models of attenuation, each leading to a different definition of $M$ (Carcione, 2007; Lakes, 2009). To illustrate attenuation, viscoelastic rheologic models are represented with springs and dashpots in one dimension. The spring represents the elastic component $M_e$ in [Pa] while the dashpot represents the viscous component $M_\eta$ in [Pa s]. The combination of these components in series or in parallel - or both - allows us to describe several viscoelastic behaviors as a function of frequency. In our work, we consider three different attenuation models: the Maxwell, the Kelvin-Voigt and the Zener viscoelastic models. They are depicted in Table 1 and we refer to Carcione (2007) for more details. In the case of the Zener model (6), one works with relaxation times $\tau_\sigma$ and $\tau_\varepsilon$ given in [s]. Note that at a fixed frequency, one can consider the different models equivalently (Imperiale et al., 2020).

### 2.3 | Viscoelastic quantities for isotropic media

In the context of an elastic isotropic medium, the viscoelastic tensor $M$ is defined from the two Lamé parameters (which are complex-valued and frequency-dependent) $\lambda$ and $\mu$. Omitting the space and frequency dependency for
<table><caption>**TABLE 1** One-dimensional representation and constitutive laws of the Maxwell, Kelvin-Voigt and Zener viscoelastic models.</caption><thead><tr><th></th><th>Maxwell model</th><th>Constitutive law</th></tr></thead><tbody><tr><td>1D representation</td><td>M<sub>e</sub> M<sub>η</sub> σ ε</td><td>(M<sub>η</sub><sup>-1</sup> + iω M<sub>e</sub><sup>-1</sup>) σ = iωε. (4)</td></tr><tr><td></td><td></td><td></td></tr><tr><td></td><td>Kelvin-Voigt model</td><td>Constitutive law</td></tr><tr><td>1D representation</td><td>M<sub>e</sub> M<sub>η</sub> σ ε</td><td>σ = M<sub>e</sub>ε + iωM<sub>η</sub>ε. (5)</td></tr><tr><td></td><td></td><td></td></tr><tr><td></td><td>Zener model</td><td>Constitutive law</td></tr><tr><td>1D representation</td><td>M<sub>e<sub>1</sub></sub> M<sub>e<sub>2</sub></sub> σ ε</td><td>(1 + iωτ<sub>σ</sub>) σ = M<sub>e</sub>(1 + iωτ<sub>ε</sub>) ε, (6)</td></tr><tr><td></td><td></td><td>with (M<sub>e<sub>1</sub></sub> + M<sub>e<sub>2</sub></sub>) M<sub>θ</sub> = M<sub>e<sub>1</sub></sub> M<sub>e<sub>2</sub></sub>, (M<sub>e<sub>1</sub></sub> + M<sub>e<sub>2</sub></sub>) τ<sub>σ</sub> = M<sub>η</sub> and M<sub>e<sub>1</sub></sub> τ<sub>ε</sub> = M<sub>η</sub>.</td></tr></tbody></table>

clarity, we have,

$$ M = \begin{pmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{pmatrix}. \quad (7) $$

For the different models presented in Table 1 under the assumption of isotropy, the elastic tensor $M_e$ and viscous tensor $M_\eta$ retain the same pattern of non-zeros as given by (7). We denote by $\lambda_e$ and $\mu_e$ the entries for $M_e$, and by $\eta_\lambda$ and $\eta_\mu$ the entries of $M_\eta$.

The wave speed of the body waves ($v_p$ for the longitudinal P-wave and $v_s$ for the shear S-wave) are defined from the Lamé parameters and therefore complex-valued, such that,

$$ v_p(\mathbf{x}, \omega) = \sqrt{\frac{\lambda(\mathbf{x}, \omega) + 2\mu(\mathbf{x}, \omega)}{\rho(\mathbf{x})}}, \quad \text{and} \quad v_s(\mathbf{x}, \omega) = \sqrt{\frac{\mu(\mathbf{x}, \omega)}{\rho(\mathbf{x})}}. \quad (8) $$

The level of attenuation is quantified by the imaginary part of the wavenumbers $k_p$ and $k_s$, respectively associated
with the P- and S-wavespeeds. We write (Carcione, 2007)

$$k_p(\mathbf{x}, \omega) = \frac{\omega}{v_p(\mathbf{x}, \omega)} = k_{\text{real}_p}(\mathbf{x}, \omega) + i\alpha_p(\mathbf{x}, \omega), \quad (9a)$$

$$k_s(\mathbf{x}, \omega) = \frac{\omega}{v_s(\mathbf{x}, \omega)} = k_{\text{real}_s}(\mathbf{x}, \omega) + i\alpha_s(\mathbf{x}, \omega), \quad (9b)$$

with the spatial attenuation parameters $\alpha_{\{p,s\}}$ given in [Np m⁻¹] (or [m⁻¹]), which causes the amplitude of a plane wave to exponentially decay in space (see equation (14a)). The temporal attenuation parameters $\beta_{\{p,s\}}$ in [Np s⁻¹] (or [s⁻¹]) are defined such that,

$$\beta_p(\mathbf{x}, \omega) = \alpha_p(\mathbf{x}, \omega)\text{Re}(v_p(\mathbf{x}, \omega)), \quad \text{and} \quad \beta_s(\mathbf{x}, \omega) = \alpha_s(\mathbf{x}, \omega)\text{Re}(v_s(\mathbf{x}, \omega)). \quad (10)$$

The quality factors (or Q-factors) associated with the P- and S- waves, respectively $Q_p$ and $Q_s$, can be rewritten as a function of the complex wavenumbers:

$$Q_p(\mathbf{x}, \omega) = \frac{\text{Re}(v_p(\mathbf{x}, \omega)^2)}{\text{Im}(v_p(\mathbf{x}, \omega)^2)}, \quad \text{and} \quad Q_s(\mathbf{x}, \omega) = \frac{\text{Re}(v_s(\mathbf{x}, \omega)^2)}{\text{Im}(v_s(\mathbf{x}, \omega)^2)}. \quad (11)$$

The quality factor can be rewritten as function of the complex wavenumbers:

$$Q_p(\mathbf{x}, \omega) = \frac{k_{\text{real}_p}(\mathbf{x}, \omega)^2 - \alpha_p(\mathbf{x}, \omega)^2}{2 k_{\text{real}_p}(\mathbf{x}, \omega) \alpha_p(\mathbf{x}, \omega)} = \frac{k_{\text{real}_p}(\mathbf{x}, \omega)}{2 \alpha_p(\mathbf{x}, \omega)} - \frac{\alpha_p(\mathbf{x}, \omega)}{2 k_{\text{real}_p}(\mathbf{x}, \omega)}, \quad (12)$$

and similarly for $Q_s$. In the context of weak attenuation, it is assumed that $k_p \gg \alpha_p$ such that we can neglect the term $\alpha_p/(2k_{\text{real}_p})$ resulting in the following approximations:

$$Q_p(\mathbf{x}, \omega) \approx \frac{k_{\text{real}_p}(\mathbf{x}, \omega)}{2 \alpha_p(\mathbf{x}, \omega)} = \frac{\omega}{2 \beta_p(\mathbf{x}, \omega)} \quad \text{and} \quad Q_s(\mathbf{x}, \omega) \approx \frac{k_{\text{real}_s}(\mathbf{x}, \omega)}{2 \alpha_s(\mathbf{x}, \omega)} = \frac{\omega}{2 \beta_s(\mathbf{x}, \omega)}. \quad (13)$$

Quality factors $Q_p$, $Q_s$ and attenuation factors $\alpha_p$, $\alpha_s$ and $\beta_p$, $\beta_s$ are used to describe the intrinsic attenuation of a wave, and the use of viscoelasticity is a mean to model this attenuation. In the following, we analyze and compare the behavior of waves for viscoelastic numerical simulations using the Maxwell and the Kelvin-Voigt rheologic models for different viscous parameters.

## 3 | METHODOLOGY

The approach we develop uses the multiple reflections of seismic waves to recover the attenuation property. We consider laboratory-sized samples of relatively weak attenuation and ultrasonic frequency waves. The waves reverberate multiple times within the sample over a short time period. The distinctive aspect of our method is that it uses the measurement of all types of seismic waves (P, S, and surface waves) in a single global signal envelope at any seismic receiver location. In our method, we take advantage of the reflections from the free-surface boundary conditions such that the wave energy is conserved within the medium, and the beam spreading effect is no longer visible in the data. In this case, the interference between the different reflected / converted waves creates an equilibrium of energy that is called equipartition, a phenomenon also used for coda waves (Ryzhik et al., 1996; Margerin et al., 2000). After
several reflections (Snieder, 2002), the energy of the wavefield is dominated by that of the S-waves. Therefore, the observed signals are overwhelmed by the S-waves and, consequently, the attenuation property we can recover will be representative of the S-waves attenuation. Furthermore, our method uses a relative amplitude decay of the signal, which makes it independent of the source wavelet, which makes it convenient for reproducibility.

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Flowchart of the signal processing for the recovery of the frequency-dependent attenuation. The numerical illustrations correspond to the configuration described in Table 2 case d, using a central frequency of 70 kHz. a) Normalized time-domain seismogram. b) Normalized frequency domain Fourier transform: [—] real part, [___] imaginary part, [---] amplitude spectrum, [---] amplitude spectrum of the Ricker seismic source. c) Gaussian filtering. d) Normalized inverse Fourier transform of c). e) Hilbert envelope of d) and $\beta(\omega_c)$ recovery.</figcaption></figure>

We present our approach for recovering the representative attenuation coefficients as a function of frequency in the flowchart of Figure 1, where we illustrate with a specific example. The data shown in Figure 1 corresponds to a 100 kHz Ricker seismic source emitted at time $t = 0$ in a medium following the viscoelastic model of case d (Table 2). Figure 1a and b show the seismic signal recorded by receiver R3 for the $u_z$ component, respectively in time
and frequency domains. In time, the seismogram displays an exponential characteristic decay of amplitude (equation (14b), modified from Carcione (2007)) of the following form:

$$ \hat{u}(\mathbf{x}, t) = \mathbf{u}_0 \cdot e^{-\alpha \mathbf{x}} \cdot e^{i(\omega t - k_{\text{real}} \mathbf{x})} \quad (14a) $$

$$ \Leftrightarrow \hat{u}(\mathbf{x}, t) = \mathbf{u}_0 \cdot e^{-\beta t} \cdot e^{i(\omega t - k_{\text{real}} \mathbf{x})}, \quad (14b) $$

with $u_0$ the initial wave amplitude, $e^{-\beta t}$ its amplitude decay with attenuation parameter $\beta$, and $e^{-i(\omega t - k_{\text{real}} \mathbf{x})}$ its oscillation at frequency $\omega$. As the seismogram is composed of a mix of P- S- and surface waves, we quantify a representative parameter $\beta$ instead of specific $\beta_p$ and $\beta_s$, respectively related to P- and S-waves.

The method proceeds as follows. (i) We need to select a small bandwidth of frequency that will evolve from low to high-frequency content. The use of a squared filter on the frequency data would lead to artifacts in time domain, this is why a Gaussian filter is used on the frequency domain data (Figure 1b) to select a relatively small bandwidth. In our case, we select a Gaussian filter of ± 10 kHz, span around the central chosen frequency $\omega_c$ (Figure 1c). (ii) After filtering, we apply an inverse Fourier transform in order to build the amplitude decay as a function of time for the small frequency bandwidth selected (Figure 1d). (iii) The next step consists in using the Hilbert envelope of the signal (Taner et al., 1979): the attenuation coefficient value $\beta$ is given by the decrease of the envelope as a function of time, and is computed through the natural logarithm of the signal (equation (14b), Figure 1e). (iv) Finally, the Q-factor is deduced from $\beta$ using equations (10) and (13). Then, we can move to the next frequency of interest and repeat these steps in order to cover the entire frequency content of the recorded data.

Selecting the appropriate width of the Gaussian frequency filter is a crucial aspect of determining the Q-factor and should be customized based on the analyzed dataset: insufficient frequencies within the chosen data window can cause instability in the inverse Fourier transform and prevent the method from accurately determining the attenuation value (Figure 2a). On the other hand, if the width of the Gaussian filter is excessively large, the decrease in amplitude observed in Figure 2c will not follow a linear trend due to the mixing of too many frequencies. Indeed, in the case of a large bandwidth, we are less sensitive to the variations of attenuation with respect to frequency. Considering that attenuation is frequency-dependent, the waves of high-attenuation frequencies decrease faster in the beginning of the seismogram, while the waves of low-attenuation frequencies remain energetic until the end of the recorded data. Currently, the optimal filter width is determined empirically in order to obtain a precise recovery of the attenuation model (Figure 2b). Further analysis would be required to analyze quantitatively how the size of the window depends on, e.g., the wavelength and signal-to-noise ratio.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Comparison of Hilbert envelopes after 3 gaussian filtering of distinct window widths. The viscoelastic simulations parameters correspond to the Table 2 case d, using the central frequency of 70 kHz. (a) ±0.5 kHz, (b) ±10 kHz, (c) ±50 kHz.</figcaption></figure>
# 4 | 3D VISCOELASTIC NUMERICAL SIMULATIONS

In this section, we detail our method to recover the representative viscoelastic properties based on wave measurements. Subsequently, we carry out numerical simulations to verify the performance of our approach. The numerical calculations are performed in samples whose dimensions are identical to the dimensions of the samples used in the laboratory experiments of Section 5.

## 4.1 | Numerical set-up configurations

We perform isotropic viscoelastic 3-D numerical simulations of wave propagation for the Maxwell and Kelvin-Voigt models of attenuation given in Table 1. The simulations are carried out through two types of geometries: a cube of size 7 cm × 7 cm × 7 cm and a cylinder of height 9 cm and diameter 5 cm (see Figure 3a and b respectively). We consider free surface boundary conditions such that,

$$ \sigma \cdot n = g, \quad \text{on the boundary,} \qquad (15) $$

where n denotes the normal direction. Here, we impose the source g on the boundary, and use a point-source positioned on the upper surface of the sample.

<figure><img src="image_4.png" /><figcaption>**FIGURE 3** Numerical set-up used for the numerical simulations. (a) 7x7x7 cm³ cube with the origin O located at the lower-left corner of the cube. The source is located at [3.5 6 7] cm on the top surface of the cube. The receiver R1 is located in the bulk of the cube at the position [2 2.5 4] cm, whereas R2 at [7 3.5 3.5] cm and R3 at [3.5 3.5 0] cm are on the surface. (b) 9 cm height and 5 cm diameter cylinder with origin O located at the center of the cylinder. The seismic source is located at [1 0.0 4.5] cm on the surface, the receiver R1 is at [1 1 1] cm inside the cylinder, R2 at [0 2.5 1] cm and R3 at [0 1 -4.5] cm on the surface.</figcaption></figure>

We carry out simulations in the frequency domain and use open-source software Haken (Faucher, 2021). The simulations are performed for frequencies between 20 Hz to 300 kHz with a step of 20 Hz. We then apply an inverse Fourier transform to obtain the time-domain seismogram for a duration of 50 ms, where we consider a Ricker wavelet source of 100 kHz peak-frequency. The synthetic signals are extracted for three positions of receiver (see Figure 3), located either inside the volume (R1) or at the boundary of the sample (R2 on a lateral face and R3 on the opposite face with respect to the source).

Five different simulations (cases a to e) are analyzed using fixed elastic parameters $λ_e = 55.95$ GPa, $μ_e = 25.95$ GPa, $ρ = 2700$ kg m⁻³ and varying viscous parameters ($η_λ, η_ν$) as detailed in Table 2: the choice of viscous parameters
may differ significantly between cases, while we maintain a comparable Q-factor using the specificities of each of the constitutive models used (Table 1). As an illustration of the variability of the attenuation induced by the Maxwell - Kelvin-Voigt - Zener viscoelastic models, Figure 4 (top) and (bottom) shows respectively the behavior of $\alpha_s$ and of the S-quality factor $Q_s$ as a function of frequencies. We see that the viscoelastic Maxwell and Kelvin-Voigt models lead respectively to a linear increase and decrease of the quality factor with frequencies, while the Zener viscoelastic model follows the behavior of the Maxwell and Kelvin-Voigt models, respectively at low and high frequencies. Moreover, we notice in Figure 4 (bottom) that around the frequency of 200 kHz, the Zener model can be used in the case of a locally constant quality factor (Imperiale et al., 2020) that can be extended to generalized viscoelastic models for a broader bandwidth of constant quality factor (Blanc et al., 2016).

<table><caption>**TABLE 2** Viscoelastic models and viscous parameters used for the numerical simulations. The elastic parameters $\lambda_e = 55.95$ GPa and $\mu_e = 25.95$ GPa are kept constant.</caption><thead><tr><th>Case</th><th>Geometry</th><th>Model</th><th>η<sub>λ</sub></th><th>η<sub>μ</sub></th><th>Q<sub>p</sub> (100kHz)</th><th>Q<sub>s</sub> (100kHz)</th></tr></thead><tbody><tr><td>a</td><td>Cube</td><td>Maxwell</td><td>100 MPa s</td><td>2.5 MPa s</td><td>612</td><td>60</td></tr><tr><td>b</td><td>Cube</td><td>Maxwell</td><td>10 MPa s</td><td>20 MPa s</td><td>291</td><td>484</td></tr><tr><td>c</td><td>Cube</td><td>Kelvin-Voigt</td><td>100 Pa s</td><td>10 Pa s</td><td>1430</td><td>4130</td></tr><tr><td>d</td><td>Cube</td><td>Kelvin-Voigt</td><td>20 Pa s</td><td>40 Pa s</td><td>1717</td><td>1033</td></tr><tr><td>e</td><td>Cylinder</td><td>Maxwell</td><td>50 MPa s</td><td>5 MPa s</td><td>350</td><td>121</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** S-waves attenuation coefficient (top) and quality-factor (bottom) behavior as a function of frequency for [---] Maxwell viscoelastic model with $\eta_\mu = 20.65$ MPa s, [---] Kelvin-Voigt viscoelastic model with $\eta_\mu = 20.65$ MPa s and [---] Zener viscoelastic model with $\tau_e = 796.57$ ns and $\tau_\sigma = 794.98$ ns. The parameters are chosen such that all attenuation models coincide at frequency 200 kHz.</figcaption></figure>
## 4.2 | Analysis of the numerical experiments

We present the results of the numerical simulations described in Table 2, which serve to validate our methodology and test its robustness. Figure 5 shows the propagating seismic wavefield at frequency 400 kHz, represented in a cross-section and at the surface. We provide the three components of the displacement, $u_x$, $u_y$, and $u_z$ in Figure 5 (left), (middle) and (right) respectively. In this specific numerical simulation, the receiver R3 is located exactly in a plane of symmetry of the domain, aligned with the seismic source. This configuration results in destructive interference for the wavefield component $u_x$ at this receiver since waves arrive from opposite directions simultaneously and cancel each other out, which means that there is no seismic energy in the $u_x$ component for receiver R3 in Figure 5. This configuration is highly improbable in practice, but with numerical simulations, singular point measurements such as this one can occur. Thus, we discard the $u_x$ component of receiver R3 from our analysis and consider the other components in our signal processing.

<figure><img src="image_6.png" /><figcaption>**FIGURE 5** Visualization of the real part of a seismic numerical wavefield generated in a cube at 400 kHz on a cross-section, lateral and upper boundaries. Amplitude of the displacement for $u_x$ (left), $u_y$ (middle) and $u_z$ (right) components. Red asterisk: source location. Blue dot: receiver R3 location.</figcaption></figure>

In the first simulation presented in Figure 6 case a, waves are propagating in a medium where the attenuation of the S-waves $\beta_S$ is 10 times higher than the attenuation of the P-waves $\beta_P$. Despite the P-waves being less attenuated, we recover numerically with our method - by far - a representative attenuation very close to the attenuation of the S-waves (see Figure 6 case a). For the sake of clarity in Figure 6, we have distinguished only the signal recorded by R1 in the bulk from the receivers R2 and R3 at the surface, with no distinctions between components x,y or z for all the receivers. The important information from Figure 6a is that we recover approximately the same attenuation value regardless of the location of the receivers and regardless of the component of the displacement field.

We have observed that the attenuation values recovered using our method consistently match the attenuation values of S-waves in the numerical simulations. This can be explained by the partitioning of energy between P- and S-waves, referred to as the *equipartition* phenomenon. It corresponds to the energy of seismic wavefield reaching a state of equilibrium after sufficient time of propagation. It has been proven analytically by Ryzhik et al. (1996) and Sánchez-Sesma et al. (2008) that after the equipartition time of propagation, the S-waves are more energetic than the P-wave. This property has also been observed numerically by Margerin et al. (2000); Sánchez-Sesma et al. (2018), and experimentally by Shapiro et al. (2000); Hennino et al. (2001); Margerin et al. (2009). In addition, Snieder (2002) further derived an analytical formula to calculate the equipartition time $\tau_{PS}$ needed to reach the state of balance
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Attenuation β (case a) and Q-factor recovery from a numerical signal generated in a 7 cm × 7 cm × 7 cm cube (cases a to d) using Maxwell viscoelastic attenuation (cases a and b) or Kelvin-Voigt viscoelastic attenuation (cases c and d) models in the media. Case e: Q recovery in a 9 cm height and 5 cm diameter cylinder with Maxwell attenuation. See Figure 3 for the geometries and Table 2 for the viscoelastic parameters. [---] Analytical βₚ or Qₚ, [---] Analytical βₛ or Qₛ, [•] recovered from the signal at the surface of the medium (receiver 2: uₓ, uᵧ and u_z components; receiver 3: uᵧ and u_z components), [•] recovered from the signal inside the medium (receiver 1: uₓ, uᵧ and u_z components).</figcaption></figure>

between the energy of the P- and S- waves in the case without attenuation, it corresponds to:

$$ \tau_{PS} = \frac{a(v_s^3 + 2v_p^3)}{2v_p v_s^3}. \qquad (16) $$

Using values for vₚ = 6320 m s⁻¹, vₛ = 3100 m s⁻¹ and a characteristic length of a = 7 cm, we expect to reach the
state of equilibrium between P- and S- waves around $τ_{PS} = 100$ µs, a time small compared to the total duration of our numerical simulations, confirming that we reach an equipartition state very fast in the numerical experiments. Indeed, the experiments we have carried out (see Table 2 and Figure 6 cases a to e) show that we systematically retrieve the attenuation of the S-waves, for simulations of Kelvin-Voigt and Maxwell viscoelastic models (that have different Q-factor behavior with respect to frequency) in a cube or a cylinder, regardless of which receiver location or which wavefield component is used for the determination of the attenuation.

Furthermore, we note that the maximum penetration depth of Rayleigh waves is of one wavelength (Kosevich and Syrkin, 1985; Soczkiewicz, 1997; Vinh and Malischewsky, 2007). In our configuration, it means that for frequencies higher than 150 kHz, the receiver R1 located at 2 cm depth captures only P- and S- body seismic waves while receivers R2 and R3 (which are on the boundaries) capture both body and surface waves. Here, we recover the same representative attenuation at all receivers. This can be explained as surface waves are formed by interferences of P- and S-waves (Takeuchi and Saito, 1972), hence after equipartition the properties of the S-waves are dominant in the observed surface waves, as it is the case for the interior wavefield.

## 4.3 | Conclusions of numerical simulations

Using numerical experiments, we have shown that for each receiver position, each component of the displacement measured, and for different geometries of samples, we can recover the frequency-dependent attenuation parameter of S-waves. This is explained by the equipartition phenomenon, that is, a state of balance of energy between propagating waves after multiple reflections and conversions within the medium. This results in the S-waves being more energetic than the P-waves, hence dominating the signals. We demonstrate that the method for the recovery of the attenuation parameter is efficient and works well with numerical simulations. In the next section, we investigate the case of experimental measurements with different rock samples.

# 5 | LABORATORY-SCALE EXPERIMENTAL MEASUREMENTS

We consider different samples of relatively low seismic attenuation and use ultrasonic wave measurements to recover their attenuation properties. The specificity of laboratory experiments is that the measurements are obtained at the surface of the sample and not in the bulk. In addition, having in mind that we want to reach the state of equipartition of the wave's energy (discussed in the numerical section), we need to record the signals for a sufficiently long time.

## 5.1 | Experimental set-up

The sample is laid on a platform as seen in Figure 7a. A seismic waveform is generated by a waveform generator Tabor 8024 and is emitted by a Panametrics piezoelectric transmitter (PZT) through the sample, the original signal being amplified by an amplifier Falco Systems WMA-300. The subsequent seismic displacement in the normal direction of the sample's surface is recorded by a laser Doppler interferometer Polytec VFX-I-120, and processed by an oscilloscope Keysight DSO-S 054A (Shen et al., 2022) on the opposite and on the lateral face of each sample with respect to the source (Figure 7b and Figure 10b). We further use two different piezoelectric sources to compare a large frequency bandwidth from 60 kHz to 1200 kHz. One source has a nominal frequency of 250 kHz (Panametrics V150) and the other one of 1 MHz (Panametrics V102). These two piezoelectric transducers are of same dimension (2.54 cm in diameter) in order to prevent biases due to uncertainties in the mechanical contact zone between the source and the
sample. Also, we used the same Gaussian filter width as in our numerical data processing, that corresponds to ± 10 kHz.

## 5.2 | Attenuation and quality factor recovery on natural samples

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

**FIGURE 7** (a) Picture of the experimental set-up with the aluminum 7 cm × 7 cm × 7 cm cube, a piezoelectric transducer (source) is located on top of the sample and the laser Doppler interferometer (receiver) records the wavefield on the lateral surface. (b) Location of the two receivers points.

### 5.2.1 | Aluminum 5083C

Aluminum 5083C is considered as a homogeneous material and low attenuating compared to natural rocks; in some cases, aluminum is even used as a non-attenuating reference material for relative measurements of attenuation (see for instance Molyneux and Schmitt (2000)). In the following, we aim to precisely recover the frequency-dependent attenuation and Q-factor of a 7 cm × 7 cm × 7 cm aluminum cube (Figure 7a) in the range of ultrasonic frequencies, and to identify the viscoelastic model followed by the sample. We measure the density $ρ = 2650 \text{ kg m}^{-3}$, and the P- and S- waves velocities respectively of $6320 \text{ m s}^{-1}$ and $3100 \text{ m s}^{-1}$ via transmission measurements between the source and the receiver.

Following the acquisition setup, the seismic signals are recorded during 50 ms with a stack of 3000 traces for the two different piezoelectric transducer sources as shown in Figure 8. The signal-to-noise ratio of the 250 kHz piezoelectric is better compared to the 1 MHz; for both experiments, the duration of the exploitable seismic signal is much longer than the estimated equipartition time (near $100 \text{ µs}$ from equation (16)). As shown in Figure 8, the frequency content of the records are from 60 to 400 kHz for the transducer of nominal frequency of 250 kHz (Figure 8a), and from 300 to 1200 kHz for the transducer of nominal frequency of 1 MHz (Figure 8b). The overlap in frequency between both experiments is used to evaluate the robustness of the method previously presented to compute the seismic attenuation.

We show in Figure 9 (top) the results of the attenuation quantification from the experiments performed in the aluminum cube, and the equivalent Q-factor in Figure 9 (bottom). The attenuation $β$ has roughly a constant value around $560 \text{ Np s}^{-1}$ throughout all the scanned frequencies, from 60 to 1200 kHz; the consequence (see equation (13)) is that the Q-factor linearly increases from 300 to 7000. Note that the data coming from the two sources and
<figure><img src="image_9.png" /><figcaption>**FIGURE 8** Experimental recorded seismogram and amplitude spectrum (normalised) of the Fourier transform for the receivers positioned on the lateral face of the aluminum cube with (a) 250 kHz PZT source and (b) 1 MHz PZT source.</figcaption></figure>

different receiver's locations recombine nicely to give quantitative and reliable values of the Q-factor with frequencies in aluminum. In order to understand and interpret the increase of the quality factor with frequencies, we remind that the definition of the Q-factor is roughly the inverse of the attenuation per wavelength (or per cycle) of a traveling seismic wave (Aki and Richards, 2002). The consequence is that for the same value of amplitude decay of a seismic wave as a function of time ($\beta$) for various frequencies, the Q-factor will increase as the wavelength decreases (see equation (13)).

Following the results shown in Figure 9 (constant $\beta$ and increasing Q-factor), we deduce that the attenuation model of the aluminum cube follows a Maxwell viscoelastic model. Using Table 1, we can further relate the Q-factor to the viscoelastic physical parameter associated with the S-waves assuming that the measured attenuation follows the S-waves attenuation. Deducing the elastic parameter $\mu_e = 25.5$ GPa from Equation (8), we compute the viscous parameter $\eta_\mu = 22.8$ MPa s for the Aluminum 5083C within the frequency bandwidth 60 kHz to 1.2 MHz using equations (4) and (11). We detail in Appendix A.1 the calculations of the viscous parameter.

### 5.2.2 | Fontainebleau sandstone

We now consider a sample of dry Fontainebleau sandstone which is well-cemented and composed of nearly 100% of quartz, such that it has a relatively low attenuation (Pimienta et al., 2015a). The sandstone sample is a cylinder of 9 cm in height and 5 cm in diameter pictured in Figure 10a, whose density is $\rho = 2500$ kg m$^{-3}$, the P-waves velocity is 5415 m s$^{-1}$ (measured) and the S-waves velocity is 3400 m s$^{-1}$ (Pimienta et al., 2018). Given the sample's height of 9 cm, the equipartition state is expected to be reached at 75 µs (equation (16)). The seismic signals recorded in the Fontainebleau Sandstone cylinder with the two piezoelectric transmitters V150 and V102 during 20 ms with a stack of 3000 traces are shown in Figure 11. We observe that the attenuation is stronger at high frequencies (PZT V102)
<figure><img src="image_10.png" /><figcaption>**FIGURE 9** Attenuation $\beta$ (top) and Q-factor (bottom) recovery on the aluminum cube from measurements: [+] on the lateral face for the 250 kHz PZT source, [+] on the opposite face for the 250 kHz PZT source, [O] on the lateral face for the 1 MHz PZT source, [O] on the opposite face for the 1 MHz PZT source. [---] corresponds to the best fit of a Maxwell model corresponding to parameters $\mu_e$ = 25.5 GPa and $\eta_\mu$ = 22.8 MPa s.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 10** (a) Picture of the Fontainebleau sandstone cylinder. (b) Location of the two receivers points.</figcaption></figure>

compared to low ones (PZT V150) since the seismic signal decays faster in time in Figure 11b compared to Figure 11a.

We show in Figure 12 the Hilbert amplitude envelope of the signal recorded on the lateral face with the 1 MHz PZT source (as in Figure 11b) for chosen central frequencies of 700 kHz ± 10 kHz (Figure 12a) and 1300 kHz ± 10 kHz (Figure 12b). We see that with the source of the central frequency 700 kHz (where the signal-to-noise ratio is relatively high), Figure 12a, we are able to recover the attenuation of the S-waves starting from the equipartition time $\tau_{PS}$ = 75 µs until the envelope reaches the noise level. On the contrary, for the central frequency of 1300 kHz (Figure 12b), the noise level is too high and prevents us from having a sufficient signal to determine the attenuation. In Figure 13, we study how the attenuation recovery may depend on the recording time of the signal. We see that we
<figure><img src="image_12.png" /><figcaption>**FIGURE 11** Experimental recorded seismogram and amplitude spectrum (normalized) on the lateral face of the Fontainebleau sandstone cylinder with (a) 250 kHz PZT source and (b) 1 MHz PZT source.</figcaption></figure>

need to record longer than the equipartition time $\tau_{PS}$ to obtain the appropriate shapes of the Hilbert envelopes and to obtain the decrease in amplitude needed to recover the appropriate attenuation. Indeed, the recorded time affects the sampling in frequency and consequently the Fourier transform from time to frequency domain, leading to a lack of information in frequency when we do not record for a sufficiently long time. This undersampling eventually affects the time-domain envelope when we perform the IFT after the Gaussian filtering used in our method (Figure 1).

<figure><img src="image_13.png" /><figcaption>**FIGURE 12** Comparison of Hilbert envelopes after Gaussian filtering on central frequency (a) 700 kHz ± 10 kHz and (b) 1300 kHz ± 10 kHz of the data recorded on the lateral face of Fontainebleau sandstone with the 1 MHz PZT source. [---] Equipartition time $\tau_{PS} = 75$µs. [.....] Estimated noise level.</figcaption></figure>

We show that the attenuation $\beta$ computed with our method increases with frequencies in Figure 14 (top), as well as the Q-factor in Figure 14 (bottom) with values from 200 to 1000. As both attenuation $\beta$ and Q-factor increase, the sandstone sample can be modeled as a Zener viscoelastic medium (Figure 4). We can then compute the
<figure><img src="image_14.png" /><figcaption>**FIGURE 13** Comparison of Hilbert envelopes after Gaussian filtering on central frequency 700 kHz ± 10 kHz depending on the recorded time $t_{max}$. [---] Equipartition time $\tau_{PS} = 75\mu s$.</figcaption></figure>

associated viscoelastic relaxation times (Carcione, 2007) for the S-waves and we obtain $\tau_{es} = 633$ ns and $\tau_{es} = 630$ ns using equations (6) and (11). The details of the computations for the Zener viscoelastic parameters are given in Appendix A.2. We further obtain a good approximation of the evolution of attenuation with ultrasonic frequencies for the dry Fontainebleau sandstone (Figure 14).

<figure><img src="image_15.png" /><figcaption>**FIGURE 14** Attenuation $\beta$ (top) and Q-factor (bottom) recovery on the Fontainebleau sandstone cylinder from measurements: [+] on the lateral face for the 250 kHz PZT source, [+] on the opposite face for the 250 kHz PZT source, [O] on the lateral face for the 1 MHz PZT source, [O] on the opposite face for the 1 MHz PZT source. [---] corresponds to the best fit of a Zener model corresponding to parameters $\tau_e = 0.633$ $\mu$s and $\tau_\sigma = 0.630$ $\mu$s.</figcaption></figure>
# 6 | DISCUSSION

We quantified the ultrasonic attenuation in a cube made of 5083C aluminum, recovering a Maxwell viscoelastic attenuation model to represent the attenuation between 60 kHz to 1.2 MHz, with quality factors from 300 to 7000. Therefore, the $Q_s$ factor is shown to increase within this frequency range. It is worth noting that, according to the authors, no other literature reports on frequency-dependent attenuation measurements at these frequencies. Nonetheless, McCann and Sothcott (1992) measured $Q_s > 700$ for 0.85 MHz frequency in an aluminum sample using the reflection method under a pressure of 60 MPa, and Hurley (1999) found the Rayleigh-waves attenuation at 10 MHz to be $1852 \text{ Np s}^{-1}$ in 6061-T651 an aluminum alloy sample; in conclusion, the attenuation values obtained in the current study are in the same range to those measured in both previous studies.

Concerning the dry Fontainebleau sandstone, we also demonstrated that $Q_s$ increases, from 200 to 1000 within the working frequencies, and its variation in frequency follows a Zener viscoelastic model. Bourbie and Zinszner (1985) measured in the same type of sandstone $Q = 100$ at a frequency of 9.5 kHz with a resonant bar method, while Pimienta et al. (2015a), with a forced oscillation method, measured a bulk Q-factor $Q > 50$; once more, the values of attenuation that we obtain for Fontainebleau sandstone are consistent with those obtained in other studies. Note that in the interpretation of our measurements developed in Section 5, we have used a single Zener model to fit an increasing attenuation $\beta_s$ and an increasing quality factor $Q_s$ throughout the range of measured frequencies. Although it would be possible to employ multiple viscoelastic systems in parallel, such as generalized attenuation models to cover wider frequency bandwidth or to model a constant quality factor (Carcione, 2007; Cao and Yin, 2014; Blanc et al., 2016), this was not within the scope of our study which focuses on the reconstruction of the ultrasonic frequency-dependent attenuation for low-attenuation media.

Finally, we have, both experimentally and numerically, observed the phenomenon of equipartition in media with attenuation. This phenomenon implies that regardless of the initial distribution of P- and S-wave sources within the system, the energy of the seismic signal is dominantly composed of S-waves after several seismic reflections/conversions at the boundaries. As a result, using our method which relies on analyzing the seismic envelope signal at a receiver location, we are only able to recover the attenuation factors associated with S-waves, due to the dominant nature of the S-waves in the seismic signal after equipartition occurs. The limitation is that we cannot access the P-waves attenuation values with our method. Nonetheless, we reconstruct the attenuation model, which characterizes both the P- and S-waves. Then, considering a small bandwidth in frequency, the attenuation of the P-waves could be obtained with other methods (such as spectral ratio methods), feeding in the attenuation model obtained with our approach.

Furthermore, in situations where the medium has high attenuation, the amplitude of the wavefield may diminish to the point of being indistinguishable from the noise level prior to achieving the equipartition state. Consequently, it is essential to initially calculate the time required for equilibrium state (Snieder, 2002) to guarantee the method's applicability to recorded seismograms. The accuracy of the reconstruction also depends on the choice of frequency windowing to obtain the appropriate frequency-dependent evolution of attenuation. The precision of the window certainly depends on the signal-to-noise ratio, sampling of the signal, and on the wavelength. This is part of ongoing investigations.

# 7 | CONCLUSIONS

We quantified the attenuation characteristics of low-attenuating samples of aluminum and dry Fontainebleau sandstone in the laboratory, using ultrasonic frequencies. Our approach uses multiple seismic wave reflections within the
samples. Prior to conducting laboratory measurements, we validated our method through 3D numerical simulations of wave propagation.

We have shown that with our method, we are able to retrieve the attenuation of the sample, for any geometry and regardless of or the receiver's location (i.e., within or at the boundaries of the sample). We have also emphasized, owing to the equipartition phenomenon, that the retrieved attenuation coefficient is associated with the S-waves. For both aluminum and Fontainebleau sandstone samples, our measurements provide new attenuation values at ultrasonic frequencies. We were able to derive the attenuation models from the recovery of $β_s$ or $Q_s$ as a function of frequency (between 60 kHz to 1.2 MHz) in aluminum and in the Fontainebleau sandstone.

More generally, our work leads the way to a better characterization of the attenuation properties in various media. In particular, the evolution of the representative attenuation as a function of frequencies in samples is a first step to the identification of structural heterogeneities, for instance, by using our findings as initial models for full waveform inversion.

## Acknowledgments

The authors thank the associate editor and reviewers for their helpful comments that have improved the manuscript. Experiments presented in this paper were carried out using open-source software *hawen* (https://ffaucher.gitlab.io/hawen-website/). We acknowledge the use of PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux and CNRS (and ANR in accordance to the programme d'investissements d'Avenir (see http://www.plafrim.fr/)). The authors wish to thank L. Pimienta for lending the Fontainebleau sandstone sample and for the scientific discussions.

## references

* Adam, L., Batzle, M., Lewallen, K. T. and van Wijk, K. (2009) Seismic wave attenuation in carbonates. *Journal of Geophysical Research: Solid Earth*, **114**.

* Agersborg, R., Johansen, T. A., Jakobsen, M., Sothcott, J. and Best, A. (2008) Effects of fluids and dual-pore systems on pressure-dependent velocities and attenuations in carbonates. *Geophysics*, **73**, N35–N47.

* Aki, K. and Richards, P. G. (2002) *Quantitative seismology*. University Science Books.

* Bardet, J. P. (1992) A Viscoelastic Model for the Dynamic Behavior of Saturated Poroelastic Soils. *Journal of Applied Mechanics*, **59**, 128–135.

* Best, A., Sothcott, J. and McCann, C. (2007) A laboratory study of seismic velocity and attenuation anisotropy in near-surface sedimentary rocks. *Geophysical Prospecting*, **55**, 609–625.

* Biot, M. A. (1954) Theory of stress-strain relations in anisotropic viscoelasticity and relaxation phenomena. *Journal of Applied Physics*, **25**, 1385–1391.

- (1956a) Theory of propagation of elastic waves in a fluid-saturated porous solid. i. low-frequency range. *The Journal of the Acoustical Society of America*, **28**, 168–178.

- (1956b) Theory of propagation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range. *The Journal of the Acoustical Society of America*, **28**, 179–191.
Blanc, E., Komatitsch, D., Chaljub, E., Lombard, B. and Xie, Z. (2016) Highly accurate stability-preserving optimization of the zener viscoelastic model, with application to wave propagation in the presence of strong attenuation. *Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society*, **205**, 427-439.

Bourbie, T. and Zinszner, B. (1985) Hydraulic and acoustic properties as a function of porosity in fontainebleau sandstone. *Journal of Geophysical Research: Solid Earth*, **90**, 11524-11532.

Bourbié, T. and Nur, A. (1984) Effects of attenuation on reflections: Experimental test. *Journal of Geophysical Research: Solid Earth*, **89**, 6197-6202.

Cao, D. and Yin, X. (2014) Equivalence relations of generalized rheological models for viscoelastic seismic-wave modeling. *Bulletin of the Seismological Society of America*, **104**, 260-268.

Carcione, J. M. (2007) *Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media*. Elsevier.

Dasgupta, R. and Clark, R. A. (1998) Estimation of q from surface seismic reflection data. *Geophysics*, **63**, 2120-2128.

Faucher, F. (2021) hawen: time-harmonic wave modeling and inversion using hybridizable discontinuous galerkin discretization. *Journal of Open Source Software*, **6**, 2699.

Faucher, F. and Scherzer, O. (2023) Quantitative inverse problem in visco-acoustic media under attenuation model uncertainty. *Journal of Computational Physics*, **472**, 111685.

Geertsma, J. and Smit, D. (1961) Some aspects of elastic wave propagation in fluid-saturated porous solids. *Geophysics*, **26**, 169-181.

Golden, J. M. and Graham, G. A. (1988) *Boundary value problems in linear viscoelasticity*. Springer Science & Business Media.

Gurevich, B. and Pevzner, R. (2015) How frequency dependency of q affects spectral ratio estimatespitfalls in q estimation. *Geophysics*, **80**, A39-A44.

Hennino, R., Trégourès, N., Shapiro, N. M., Margerin, L., Campillo, M., van Tiggelen, B. A. and Weaver, R. L. (2001) Observation of equipartition of seismic waves. *Phys. Rev. Lett.*, **86**, 3447-3450.

Hurley, D. C. (1999) Nonlinear propagation of narrow-band rayleigh waves excited by a comb transducer. *The Journal of the Acoustical Society of America*, **106**, 1782-1788.

Imperiale, A., Leymarie, N. and Demaldent, E. (2020) Numerical modeling of wave propagation in anisotropic viscoelastic laminated materials in transient regime: Application to modeling ultrasonic testing of composite structures. *International Journal for Numerical Methods in Engineering*, **121**, 3300-3338.

Jackson, I. (1993) Progress in the experimental study of seismic wave attenuation. *Annual Review of Earth and Planetary Sciences*, **21**, 375-406.

Johnston, D. H., Toksöz, M. N. and Timur, A. (1979) Attenuation of seismic waves in dry and saturated rocks: li. mechanisms. *Geophysics*, **44**, 691-711.

Knopoff, L. and MacDonald, G. J. (1958) Attenuation of small amplitude stress waves in solids. *Reviews of Modern Physics*, **30**, 1178.

Kosevich, Y. A. and Syrkin, E. (1985) Existence criterion and properties of deeply penetrating rayleigh waves in crystals. *Zhurnal Ehksperimental'noj i Teoreticheskoj Fiziki*, **89**, 2221-2229.

Lakes, R. S. (2009) *Viscoelastic materials*. Cambridge university press.
Leong, E.-C., Yeo, S.-H. and Rahardjo, H. (2004) Measurement of wave velocities and attenuation using an ultrasonic test system. *Canadian geotechnical journal*, **41**, 844-860.

Manthei, G., Eisenblatter, J. and Spies, T. (2006) Determination of wave attenuation in rock salt in the frequency range 1-100 khz using located acoustic emission events. *Progress in Acoustic Emission*, **13**, 129.

Margerin, L., Campillo, M. and van Tiggelen, B. (2000) Monte Carlo simulation of multiple scattering of elastic waves. *Journal of Geophysical Research: Solid Earth*, **105**, 7873-7892. URL: https://hal-insu.archives-ouvertes.fr/insu-03606316.

Margerin, L., Campillo, M., Van Tiggelen, B. and Hennino, R. (2009) Energy partition of seismic coda waves in layered media: theory and application to pinyon flats observatory. *Geophysical Journal International*, **177**, 571-585.

McCann, C. and Sothcott, J. (1992) Laboratory measurements of the seismic properties of sedimentary rocks. *Geological Society, London, Special Publications*, **65**, 285-297.

Molero, M., Segura, I., Aparicio, S., Hernández, M. G. and Izquierdo, M. (2010) On the measurement of frequency-dependent ultrasonic attenuation in strongly heterogeneous materials. *Ultrasonics*, **50**, 824-828.

Molyneux, J. B. and Schmitt, D. R. (2000) Compressional-wave velocities in attenuating media: A laboratory physical model study. *Geophysics*, **65**, 1162-1167.

Morochnik, V. and Bardet, J. (1996) Viscoelastic approximation of poroelastic media for wave scattering problems. *Soil Dynamics and Earthquake Engineering*, **15**, 337-346.

Müller, T. M., Gurevich, B. and Lebedev, M. (2010) Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks – a review. *Geophysics*, **75**, 75A147-75A164.

Norris, D. M. and Young, W.-C. (1970) Complex-modulus measurement by longitudinal vibration testing. *Experimental Mechanics*, **10**, 93-96.

Pandit, B. and King, M. (1979) The variation of elastic wave velocities and quality factor q of a sandstone with moisture content. *Canadian Journal of Earth Sciences*, **16**, 2187-2195.

Pang, M., Ba, J., Carcione, J. M., Picotti, S., Zhou, J. and Jiang, R. (2019) Estimation of porosity and fluid saturation in carbonates from rock-physics templates based on seismic q. *Geophysics*, **84**, M25-M36.

Pimienta, L., Borgomano, J. V. M., Fortin, J. and Guéguen, Y. (2017) Elastic dispersion and attenuation in fully saturated sandstones: Role of mineral content, porosity, and pressures. *Journal of Geophysical Research: Solid Earth*, **122**, 9950-9965.

Pimienta, L., Fortin, J. and Guéguen, Y. (2015a) Bulk modulus dispersion and attenuation in sandstones. *Geophysics*, **80**, D111-D127.

— (2015b) Experimental study of young's modulus dispersion and attenuation in fully saturated sandstones. *Geophysics*, **80**, L57-L72.

Pimienta, L., Schubnel, A., Violay, M., Fortin, J., Guéguen, Y. and Lyon-Caen, H. (2018) Anomalous vp/vs ratios at seismic frequencies might evidence highly damaged rocks in subduction zones. *Geophysical Research Letters*, **45**, 12,210–12,217. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018GL080132.

Quan, Y. and Harris, J. M. (1997) Seismic attenuation tomography using the frequency shift method. *Geophysics*, **62**, 895-905.

Rao, M., Sarma, L. and Lakshmi, K. (2002) Ultrasonic pulse-broadening and attenuation in volcanic rock-a case study. *IJPAP*.

Ryzhik, L., Papanicolaou, G. and Keller, J. B. (1996) Transport equations for elastic and other waves in random media. *Wave Motion*, **24**, 327-370.
Shapiro, N. M., Campillo, M., Margerin, L., Singh, S. K., Kostoglodov, V. and Pacheco, J. (2000) The Energy Partitioning and the Diffusive Character of the Seismic Coda. *Bulletin of the Seismological Society of America*, **90**, 655–665.

Shen, C., Brito, D., Diaz, J., Sanjuan, F., Bordes, C. and Garambois, S. (2022) Pulsed-laser source characterization in laboratory seismic experiments. *Geomechanics and Geophysics for Geo-Energy and Geo-Resources*, **8**, 1–26.

Snieder, R. (2002) Coda wave interferometry and the equilibration of energy in elastic media. *Phys. Rev. E*, **66**, 046615. URL: https://link.aps.org/doi/10.1103/PhysRevE.66.046615.

Soczkiewicz, E. (1997) The penetration depth of the rayleigh surface waves. *Nondestructive Testing and Evaluation*, **13**, 113–119.

Spencer Jr., J. W. (1981) Stress relaxations at low frequencies in fluid-saturated rocks: Attenuation and modulus dispersion. *Journal of Geophysical Research: Solid Earth*, **86**, 1803–1812.

Subramaniyan, S., Quintal, B., Tisato, N., Saenger, E. H. and Madonna, C. (2014) An overview of laboratory apparatuses to measure seismic attenuation in reservoir rocks. *Geophysical Prospecting*, **62**, 1211–1223.

Sánchez-Sesma, F. J., Pérez-Ruiz, J. A., Luzón, F., Campillo, M. and Rodríguez-Castellanos, A. (2008) Diffuse fields in dynamic elasticity. *Wave Motion*, **45**, 641–654. Special Issue on Selected Papers Presented at the International Symposium on Mechanical Waves in Solids - Volume 2.

Sánchez-Sesma, F. J., Victoria-Tobon, E., Carbajal-Romero, M., Rodríguez-Sánchez, J. E. and Rodríguez-Castellanos, A. (2018) Energy equipartition in theoretical and recovered seismograms. *Journal of Applied Geophysics*, **150**, 153–159.

Takeuchi, H. and Saito, M. (1972) Seismic surface waves. *Methods in computational physics*, **11**, 217–295.

Taner, M. T., Koehler, F. and Sheriff, R. (1979) Complex seismic trace analysis. *Geophysics*, **44**, 1041–1063.

Tisato, N. and Madonna, C. (2012) Attenuation at low seismic frequencies in partially saturated rocks: Measurements and description of a new apparatus. *Journal of Applied Geophysics*, **86**, 44–53. URL: https://www.sciencedirect.com/science/article/pii/S0926985112001322.

Tisato, N., Quintal, B., Chapman, S., Madonna, C., Subramaniyan, S., Frehner, M., Saenger, E. H. and Grasselli, G. (2014) Seismic attenuation in partially saturated rocks: Recent advances and future directions. *The Leading Edge*, **33**, 640–646.

Toksoz, M. N., Johnston, D. H. and Timur, A. (1979) Attenuation of seismic waves in dry and saturated rocks: I. laboratory measurements. *Geophysics*, **44**, 681–690.

Vinh, P. C. and Malischewsky, P. G. (2007) An approach for obtaining approximate formulas for the rayleigh wave velocity. *Wave Motion*, **44**, 549–562.

Winkler, K. W. and Nur, A. (1979) Pore fluids and seismic attenuation in rocks. *Geophysical Research Letters*, **6**, 1–4.

– (1982) Seismic attenuation: Effects of pore fluids and frictional-sliding. *Geophysics*, **47**, 1–15.

Winkler, K. W. and Plona, T. J. (1982) Technique for measuring ultrasonic velocity and attenuation spectra in rocks under pressure. *Journal of Geophysical Research: Solid Earth*, **87**, 10776–10780.
# A | APPENDIX: RECOVERY OF VISCOUS PROPERTIES

In this appendix we detail how to calculate the viscous parameters from the frequency-dependent Q-factor.

## A.1 | Maxwell model

Using the representation of the viscoelastic tensor in (7) and the constitutive law of the Maxwell model (4), we can express the complex-valued bulk modulus associated to the Maxwell attenuation model, $\mu_{Max}$, in terms of the elastic parameter $\mu_e$ and the viscous parameter $\eta_\mu$, such that,

$$ \mu_{Max} = \frac{\mu_e \eta_\mu^2 \omega^2 + i \omega \mu_e^2 \eta_\mu}{\eta_\mu^2 \omega^2 + \mu_e^2}. \quad (17) $$

The frequency-dependent Q-factor for the S-waves is given from equation (11) and writes as:

$$ Q_{sMax} = \frac{\mathrm{Re}(\mu_{Max})}{\mathrm{Im}(\mu_{Max})} = \omega \frac{\eta_{\mu}}{\mu_{e}}. \quad (18) $$

Therefore, the S-wave quality factor for the Maxwell model increases linearly with frequency. Since the elastic parameter $\mu_e$ can be recovered from the S-waves velocity (equation (8)), the remaining unknown in equation (18) is $\eta_\mu$, which is calculated by computing the linear function fitting the measured value of Q as a function of frequency (Figure 9).

## A.2 | Zener model

The complex parameter associated to the Zener model, $\mu_{Zener}$ is obtained from (6) such that,

$$ \mu_{Zener} = \mu_e \frac{1 + i\omega\tau_{es}}{1 + i\omega\tau_{es}}, \quad (19) $$

with $\mu_e = \frac{\mu_{e1}\mu_{e2}}{\mu_{e1} + \mu_{e2}}$, $\tau_{es} = \frac{\eta_\mu}{\mu_{e2}}$ and $\tau_{es} = \frac{\eta_\mu}{\mu_{e1} + \mu_{e2}}$.

We rewrite the Q-factor for the S-waves from equation (11) for the Zener viscoelastic model to obtain:

$$ Q_{sZener} = \frac{1 + \omega^2 \tau_{es} \tau_{ss}}{\omega(\tau_{es} - \tau_{ss})} = \frac{1}{\omega(\tau_{es} - \tau_{ss})} + \frac{\omega \tau_{es} \tau_{ss}}{\tau_{es} - \tau_{ss}} = \frac{a}{\omega} + bw, \quad (20) $$

with $a = \frac{1}{\tau_{es} - \tau_{ss}}$ and $b = \frac{\tau_{es} \tau_{ss}}{\tau_{es}}$.

We can compute the values a and b by fitting the function $Q_{sZener}(\omega) = \frac{a}{\omega} + bw$ to the measured Q values as a function of frequency (Figure 14). Once the values of a and b are retrieved, we can recover the parameters $\tau_{es}$ and
$\tau_{\sigma_s}$ by developing the system of equations:

$$
\left\{
\begin{array}{l}
a = \frac{1}{\tau_{e_s} - \tau_{\sigma_s}} \\
b = \frac{\tau_{e_s} \tau_{\sigma_s}}{\tau_{e_s} - \tau_{\sigma_s}}
\end{array}
\right.
\qquad
\Leftrightarrow
\qquad
\left\{
\begin{array}{l}
\tau_{e_s} = \frac{1}{a} + \tau_{\sigma_s} \\
b = \frac{\tau_{e_s} \tau_{\sigma_s}}{\tau_{e_s} - \tau_{\sigma_s}}
\end{array}
\right.
\qquad (21)
\\[1.5em]
\Leftrightarrow
\qquad
\left\{
\begin{array}{l}
\tau_{e_s} = \frac{1}{a} + \tau_{\sigma_s} \\
b = \tau_{\sigma_s} + a \tau_{\sigma_s}^2.
\end{array}
\right.
$$

We solve the second-order equation such that,

$$a\tau_{\sigma_s}^2 + \tau_{\sigma_s} - b = 0. \tag{22}$$

Using that $\tau_{e_s}$ and $\tau_{\sigma_s}$ are $\ge 0$ (Carcione, 2007) and following the same development as equation (21) for $\tau_{e_s}$, we have

$$
\begin{align}
\tau_{e_s} &= \frac{1 + \sqrt{1 + 4ab}}{2a}, \tag{23} \\
\tau_{\sigma_s} &= \frac{-1 + \sqrt{1 + 4ab}}{2a}. \notag
\end{align}
$$