# Inversion in viscoelastic media. Some effects of the memory kernel singularity on wave propagation and inversion in viscoelastic media, II: inversion

A. Ribodetti, A. Hanyga

► To cite this version:

A. Ribodetti, A. Hanyga. Inversion in viscoelastic media. Some effects of the memory kernel singularity on wave propagation and inversion in viscoelastic media, II: inversion. Geophysical Journal International, 2005, 158 (2), pp.426-442. 10.1111/j.1365-246X.2004.02337.x. hal-00407461

HAL Id: hal-00407461

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

Submitted on 19 Oct 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 CC BY 4.0 - Attribution - International License
# Some effects of the memory kernel singularity on wave propagation and inversion in viscoelastic media – II. Inversion

A. Ribodetti¹ and A. Hanyga²

¹IRD-GEOSCIENCES AZUR UMR6526, Observatoire Océanologique de Villefranche BP 48, Quai de la Darse 06235, Villefranche-sur-Mer Cedex France

²Department of Earth Sciences, University of Bergen, Allégaten 41, NO-5007 Bergen, Norway

Accepted 2004 April 4. Received 2004 March 24; in original form 2003 May 16

## SUMMARY

Viscoelastic inversion is developed for a realistic, simple, causal hereditary model with a power-law attenuation. In this class of models the energy travels with a delay behind the wave front because the effective speed of the seismic signal is lower than the infinite-frequency limit of the propagation speed which determines the wave front propagation. In inversion this discrepancy affects the correct positioning of scatterers and interfaces. Our inversion method is an extension of Ribodetti’s imaging method. Using a power-law model, propagation- and attenuation-related parameters are uncoupled for the acquisition geometries used in laboratory and in field seismic experiments. The method is applied to ultrasonic laboratory data where we have complete control of the acquisition parameters and the physical properties of the medium to be recovered are well known. An application to ultrasonic data for a wavefield scattered from a hollow PVC sample and a lava sample immersed in a water tank demonstrates that the proposed inversion method allows reliable parameter estimation in the power-law class of viscoacoustic models.

**Key words:** inversion, lava sample, seismic attenuation, ultrasonic laboratory data, viscoelasticity.

## NOTATION

$$ \hat{f}(\omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{i\omega t} f(t) dt $$

$$ \tilde{f}(s) = \int_{0}^{\infty} e^{-st} f(t) dt $$

## 1 INTRODUCTION

Realistic inversion of seismic wave data requires a model which is mathematically rigorous, physically acceptable and sufficiently flexible. Such a model has to include absorption and dispersion effects and must also satisfy the causality principle. Power laws of attenuation (Szabo 1994, 1995) are the simplest models of causal and attenuating media that include absorption and dispersion effects. Power-law models can be associated with a subset of the class of viscoelastic equations of motion with singular memory, studied in a related paper (Hanyga & Seredyńska 1999a). A more detailed mathematical study of this class of integro-differential equations can be found in (Lokshin & Rok 1978a; Hanyga & Rok 2000; Hanyga & Seredyńska 2002). Solutions of the associated initial-value problems are infinitely smooth at the wave fronts. Since the solution vanishes ahead of the wave front, an initial singular delta-spiked pulse sent by the source immediately spreads out and its maximum lags behind the wave front. The peaks of a smooth causal signal arrive with a fixed delay with respect to the wave front, determined by the properties of the signal only. In power-law models there is an additional delay depending on the attenuating properties of the medium and increasing with the propagation distance. Since the energy travels with an ever-increasing delay behind the wave front, the effective speed of the seismic signal is lower than the infinite-frequency limit of the propagation speed associated with the wave front. This has some importance for imaging scatterers such as interfaces. As shown by Hanyga & Seredyńska (1999a), for a given source, receiver and a singly scattered signal arrival time the location of the scatterer lies in a thin layer above the isochrone determined from wave front (or traveltime) considerations.

The principal objective of this paper is to present a method of viscoelastic inversion for power-law viscoacoustic models. Seismic inversion (or its first, imaging stage) is usually based on the Beylkin method of interface imaging (Beylkin 1985). It is, however, easy to see
that the original Beylkin method does not apply to attenuating media due to the failure of a symmetry of the phase function as a function of frequency. We shall therefore apply the method of viscoelastic inversion developed by Ribodetti et al. (2000).

In an attenuating medium the phase of a signal consists of an imaginary part, which is an odd function of frequency, and a negative real part, which is an even function of frequency. The Beylkin single-step inversion method (Beylkin 1985), commonly used in elastic imaging, implicitly assumes that the phase is an odd function of frequency. An extension of Beylkin's imaging method, adapted for dealing with viscoelastic media, was recently developed by Ribodetti et al. (2000) and applied to a non-causal constant-Q model. Ribodetti's method takes into account the lack of symmetry between positive and negative frequencies in attenuating media.

In this paper we apply the method of Ribodetti et al. (2000) to power-law attenuation. In causal hereditary models the appearance of an additional imaginary dispersive phase associated with the attenuation factor causes some additional complications in the transformation from the frequency variable to the composite wavenumber vector length, as required by Beylkin's method.

The inversion scheme is tested on data acquired from laboratory measurements of the acoustic wavefield scattered by a hollow PVC cylinder and by a lava cylinder immersed in a water tank. It is demonstrated that the inversion allows a complete parameter estimation in the power-law class of acoustic models.

The power-law model is presented in Section 2, while the inversion method is described in Sections 3 and 4. An application to laboratory data (a PVC and lava specimen immersed in a water tank) is discussed in the following sections.

## 2 MODEL

A scalar model of single-mode wave propagation with power-law attenuation can be represented by a frequency-domain expression:

$$ \hat{u}(\omega, \mathbf{x}) \propto f(\omega) \exp[i\omega T(\mathbf{x}) - (i\omega)^\alpha A(\mathbf{x})] / (4\pi|\mathbf{x}|) \quad (1) $$

with $A \ge 0$, $0 < \alpha < 1$, and a source spectrum $f(\omega)$. It is easy to see that the inverse Fourier transform $u(t, \mathbf{x})$ of $\hat{u}(\omega, \mathbf{x})$ vanishes for $t > T$ provided the function $f(\omega)$ is analytic in the upper complex half-plane (Hanyga & Seredyńska 2002). In fact, for $f = 1$ and $T = 0$ the right-hand side of (1) is a totally skewed $\alpha$-stable probability distribution, vanishing identically on the negative real axis (Uchaikin & Zolotarev 1999; Kreis & Pipkin 1986). Since we shall only need Green's functions for a homogeneous background medium, caustics are absent. Following the method of Appendix D in Hanyga & Seredyńska (2002) one can easily show that eq. (1) is the dominating term of the frequency domain solution

$$ \hat{u}(s, \mathbf{x}) = \frac{1 + \tilde{K}(s)}{4\pi r} \exp[-rs\phi(s)/c] \quad (2) $$

of the initial-value problem

$$ u_{tt} + K(t) * u_{tt} = c^2 \nabla^2 u + \delta(t) \delta(\mathbf{x}); \quad u(0, \mathbf{x}) = u_t(0, \mathbf{x}) = 0 \quad (3) $$

where $s = -i\omega$ is the Laplace-transform variable conjugate for $t$,

$$ \phi(s) = 1 + (\tau s)^{\alpha-1} $$

and

$$ 1 + \tilde{K}(s) = \phi(s)^2 \equiv 1 + 2(\tau s)^{\alpha-1} + (\tau s)^{2\alpha-2} \quad (4) $$

with a constant $\tau > 0$ of dimension $[T]$. The kernel $K(t)$, given by the formula

$$ \left[ \frac{2(t/\tau)^{-\alpha}}{\Gamma(1-\alpha)} + \frac{(t/\tau)^{1-2\alpha}}{\Gamma(2-2\alpha)} \right] / \tau \quad (5) $$

is singular at $t = 0$ and locally integrable. Eqs (3) with the kernel (5) were originally studied by Lokshin & Rok (1978a,b) (cf. also Lokshin & Suvorova 1982; Dautray & Lions 1992; Hanyga & Seredyńska 1999b, 2002; Hanyga & Rok 2000; Engler 1997). Since the kernel $\delta(t) + K(t)$ is invertible in the convolution algebra, it is possible to define an equivalent stress relaxation function at the expense of working with the Mittag–Leffler function. Consequently, the medium defined by eqs (3) and (5) is viscoacoustic.

Eq. (2) implies that complex-valued wave propagation velocity for the model is given by the formula

$$ v(\omega) = c/\phi(-i\omega) $$

while the quality factor

$$ Q^{-1} = 4\pi \cos(\pi\alpha/2)(\omega\tau)^{\alpha-1}. $$

Since the kernel $\delta(t) + K(t)$ is invertible in the convolution algebra, eq. (3) is equivalent to a viscoacoustic equation

$$ \rho u_{tt} = c^2 \nabla^2 u + c^2 f(t) * \nabla^2 u + [1 + f(t)] \delta(\mathbf{x}); \quad u(0, \mathbf{x}) = u_t(0, \mathbf{x}) = 0. \quad (6) $$

The stress response kernel $f(t)$ can be expressed in terms of the generalized Mittag–Leffler functions (Engler 1997). The creep function of the corresponding viscoelastic material is, however, quite simple (Engler 1997)

$$ J(t) = 1 + 2 \frac{(t/\tau)^{1-\alpha}}{\Gamma(2-\alpha)} + \frac{(t/\tau)^{2(1-\alpha)}}{\Gamma(3-2\alpha)}. $$

The causality and dissipative nature of power-law models for $0 < \alpha < 1$ is proved in Hanyga & Seredyńska (2003).
The inverse Laplace transform of $\tilde{u}(s, x)$ can be expressed in terms of Wright functions (Podlubny 1998). For $\alpha = 1/2, 1/3, 2/3$ the time-domain solution can be expressed in terms of easily computable functions (elementary functions, the complementary error function and the Airy functions, Hanyga & Seredyńska (2002)). For slightly more general kernels

$$K(t) = \left[ \frac{2a(\mathbf{x})(t/\tau)^{-\alpha}}{\Gamma(1-\alpha)} + \frac{b(\mathbf{x})(t/\tau)^{1-2\alpha}}{\Gamma(2-2\alpha)} \right] / \tau \quad (7)$$

the initial-value problems for eq. (3) can be approximately solved by the Born approximation (Hanyga & Seredyńska 1999a) or, for $\alpha = 1/2, 1/3, 2/3$, by ray-asymptotic methods (Hanyga & Seredyńska 2002). Efficient finite difference and pseudospectral methods can be based on an equivalent formulation in terms of fractional time derivatives (Podlubny 1998; Lu & Hanyga 2004). For applications in seismic inversion the Green's function for the background model should preferably be computable in the frequency domain. This leaves the following possibilities for the background model: it is either an arbitrary power-law model (3) with the kernel $K$ given by eq. (5) and constant parameters, or a possibly heterogeneous model for which a ray-asymptotic Green's function is available, defined by the time convolution kernels (7) with $\alpha = 1/2, 1/3, 2/3$.

### 3 BORN APPROXIMATION

Let $G_0(t, \mathbf{x})$ denote the Green's function for the background medium, i.e. the solution of eq. (3) with $c^2 = \kappa_0/\rho_0$, with positive constants $\kappa_0, \rho_0$,

$$K_0(t) = 2 \frac{(t/\tau)^{\alpha_0-1}}{\Gamma(\alpha_0)} + \frac{(t/\tau)^{2\alpha_0-1}}{\Gamma(2\alpha_0)} \quad (8)$$

a constant $0 < \alpha_0 < 1$ and a positive constant $\tau$.

The perturbed medium is characterized by the memory kernel $K(t, \mathbf{x}) = K_0(t) + \delta K(t, \mathbf{x})$, and the parameter $\kappa(\mathbf{x}) = \kappa_0 + \delta\kappa(\mathbf{x})$. The corresponding Green's function satisfies the equation

$$G_{tt} + K * G_{tt} = (1/\rho)\nabla \cdot [\kappa \nabla G] + \delta(t)\delta(\mathbf{x}); \quad u(0, \mathbf{x}) = u_t(0, \mathbf{x}) = 0 \quad (9)$$

or, equivalently

$$G_{tt} + K_0 * G_{tt} - (1/\rho_0)\nabla \cdot [\kappa \nabla G] = \delta(t)\delta(\mathbf{x}) \\ - \delta K * G_{tt} + \delta\rho^{-1}\nabla \cdot [\kappa \nabla G] + (1/\rho)\nabla \cdot [\delta\kappa \nabla G] \quad (10) \\ u(0, \mathbf{x}) = u_t(0, \mathbf{x}) = 0.$$

We shall consider the kernel

$$K(t, \mathbf{x}) = 2a(\mathbf{x}) \frac{(t/\tau)^{\alpha_0-1}}{\Gamma(\alpha_0)} + b(\mathbf{x}) \frac{(t/\tau)^{2\alpha_0-1}}{\Gamma(2\alpha_0)} \quad (11)$$

with $a(\mathbf{x})$, $b(\mathbf{x})$ close to unity and $\alpha \approx \alpha_0$.

The Green's function $G(t, \mathbf{x})$ of eq. (10) can be expressed in terms of $G_0(t, \mathbf{x})$:

$$G(\mathbf{y}, \mathbf{x}, t_1 - t_0) = G_0(\mathbf{y}, \mathbf{x}, t_1 - t_0) + \int G_0(\mathbf{y}, \mathbf{y}', t) \left[ -\delta K(t, \mathbf{y}') * G_{tt}(\mathbf{y}', \mathbf{x}, t) + \rho_0^{-1} \nabla \cdot [\delta \kappa(\mathbf{y}') \nabla G_0(\mathbf{y}', \mathbf{x}, t)] \right]_{t=t_1-t_0} d\mathbf{y}' \quad (12)$$

where the asterisk denotes the time convolution and

$$\delta K(t, \mathbf{x}) = 2[a(\mathbf{x}) - 1] \frac{(t/\tau)^{\alpha_0-1}}{\Gamma(\alpha_0)} + [b(\mathbf{x}) - 1] \frac{(t/\tau)^{2\alpha_0-1}}{\Gamma(2\alpha_0)} + 2 \left[ \ln(t/\tau) - (t/\tau)^{\alpha_0-1} \psi(\alpha_0) \right] (\alpha - \alpha_0) / \Gamma(\alpha_0)$$

where $\psi(z) = d \ln \Gamma(z)/dz$.

In the Born approximation

$$G(\mathbf{y}, \mathbf{x}, t_1 - t_0) \approx G_B(\mathbf{y}, \mathbf{x}, t_1 - t_0) := G_0(\mathbf{y}, \mathbf{x}, t_1 - t_0) \\ + \int G_0(\mathbf{y}, \mathbf{y}', t) \left[ -\delta K(t, \mathbf{y}') * G_{0,tt}(\mathbf{y}', \mathbf{x}, t) + \rho_0^{-1} \nabla \cdot [\delta \kappa(\mathbf{y}') \nabla G_0(\mathbf{y}', \mathbf{x}, t)] \right] d\mathbf{y}'. \quad (13)$$

In the frequency domain the Born approximation assumes the form

$$\widehat{G}_B(\omega, \mathbf{y}, \mathbf{x}) = \widehat{G}_0(\omega, \mathbf{y}, \mathbf{x}) + \int \widehat{G}_0(\omega, \mathbf{y}', \mathbf{x}) \{\omega^2 \widehat{G}_0(\omega, \mathbf{z}, \mathbf{y}') \delta \tilde{K}(\omega, \mathbf{y}') \\ + \delta \rho^{-1}(\mathbf{y}') \nabla \cdot [\kappa_0(\mathbf{y}') \nabla \widehat{G}_0(\mathbf{y}', \mathbf{x}, \omega)] + \rho_0(\mathbf{y}')^{-1} \nabla \cdot [\delta \kappa(\mathbf{y}') \nabla \widehat{G}_0(\omega, \mathbf{y}', \mathbf{x})]\} \Big|_{t=t_1-t_0} d\mathbf{y}'. \quad (14)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Experimental set-up: the water tank, the data acquisition system and the PVC cylinder are shown. The horizontal plane at a depth of z = 0.67 m beneath the water level containing source and receiver hydrophones and the horizontal section of the cylinder is plotted.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Acquisition geometry. (a) Schematic top view of the source-receiver system in a horizontal plane at z = 0.67 m. The target is the section of the cylinder. (b) Close-up of (a) centred on the target. The dimensions of the cylinder and the size of the numerical grid used for inversion are indicated.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Seventy-two seismic traces of the scattered wavefield from the PVC cylinder: red, recorded data; blue, reconstructed data after inversion for $\alpha = 0.5$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of recorded and synthetic seismograms for PVC for a guessed value of $\alpha = 0.4$.</figcaption></figure>
# 4 INVERSION

## 4.1 A generalized Beylkin-type inversion operator

We shall focus on the parameter $κ = λ + 2μ$ (which is bulk modulus for water; $λ$ and $μ$ denote the Lamé constants) and the attenuation parameter $a$ for water and for the scatterer. The parameter $b$, appearing in a smoother part of the perturbed memory kernel, will be neglected.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of recorded and synthetic seismograms for the PVC cylinder; guessed value $α = 0.6$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Recovered $δa$ and $δk$ sections for the PVC cylinder after one iteration for an assumed value of $α = 0.5$. The shape of the cylinder is clearly identifiable on both sections notwithstanding a contamination by the source signature and multiples.</figcaption></figure>
The parameter $\alpha$ is determined by the physics of the relaxation mechanism at the microlevel and therefore it is expected to be a material constant. Furthermore, it is unlikely to exhibit a significant continuous variation in a medium or part of it (a phase) consisting of a single material. Whenever the wave propagation occurs in one material it is convenient to hold $\alpha$ at a constant value while applying the inversion procedure; the value of $\alpha$ is determined a posteriori by optimizing the discrepancy between the recorded data and the synthetic data obtained from the parameters of the medium determined by inversion. In those cases where the scattered field is essentially reflected by an obstacle and there is no significant transmission through the scattering obstacle it is also possible to guess the appropriate value of $\alpha$ for the surrounding medium and for the scatterer before running inversion for the remaining parameters and check it a posteriori against alternative values.

In an experiment involving transmission through two media the exponent $\alpha$ is expected to assume two different values. In this case inversion can be applied to determine the value of $\alpha$ and the remaining parameters directly. We shall accordingly consider two- and three-parameter inversion, the latter involving the parameter $\alpha$ as an additional inversion target; it will be seen below that the local values of $\alpha$ can be used for imaging the target.

In view of eq. (1) the phase of the Green's function

$$ \Phi = iT(\mathbf{x}) + i\tau(-i\omega\tau)^{\alpha-1} A(\mathbf{x}) = \Phi_R + i\Phi_I \quad (15) $$

consists of an imaginary part

$$ \Phi_I = T(\mathbf{x}) + \operatorname{sgn}(\omega)\tau|\omega\tau|^{\alpha-1} \sin(\pi\alpha/2) A(\mathbf{x}) \quad (16) $$

and a real part

$$ \Phi_R = -\tau|\omega\tau|^{\alpha-1} \cos(\pi\alpha/2) A(\mathbf{x}). \quad (17) $$

$T$ denotes the traveltime, i.e. the arrival time of the wave front. While constructing the inversion operator the slowly varying factor $\exp(\omega\Phi_R)$ will be included in the amplitude, while the imaginary phase will be used in the construction of the inverse operator.

Let $\mathbf{x}(\xi)$, $\mathbf{z}(\xi)$ denote the set of sources and receivers parametrized by $\xi$ varying over a region $U \subset \mathbb{R}^2$ of the plane. The coordinates of the scatterers are denoted by $\mathbf{y}$. The total imaginary phase of a singly scattered wave in eq. (14) equals

$$ \Phi_{IT}(\mathbf{y}) = [T(\mathbf{x}, \mathbf{y}) + T(\mathbf{z}, \mathbf{y})] + \operatorname{sgn}(\omega)\tau(\omega\tau)^{\alpha-1} \sin(\pi\alpha/2) [A(\mathbf{x}, \mathbf{y}) + A(\mathbf{z}, \mathbf{y})]. \quad (18) $$

<figure><img src="image_8.png" /><figcaption>Figure 7. Rms for the PVC cylinder, $\alpha = 0.4, 0.5, 0.6$.</figcaption></figure>
The unperturbed Green's function $\hat{G}_0(\mathbf{z}, \mathbf{x}; \omega) = D(\mathbf{x}, \mathbf{y}) \exp(i\omega\Phi_1)$ where $D(\mathbf{x}, \mathbf{y}) = B(\mathbf{x}, \mathbf{y}) C(\mathbf{x}, \mathbf{y}); B(\mathbf{x}, \mathbf{y})$ denotes the geometric divergence factor (= $1/4 \pi |\mathbf{x} - \mathbf{y}|$ in a homogeneous background medium) and $C(\mathbf{x}, \mathbf{y}) = \exp[i\omega\Phi_R(\mathbf{x}, \mathbf{y})]$. After an asymptotic evaluation of the Laplacian the Born approximation can be expressed in a compact form

$$ \delta \hat{G}(\mathbf{z}, \mathbf{x}, \omega) := G_B(\mathbf{z}, \mathbf{x}, \omega) - G_0(\mathbf{z}, \mathbf{x}, \omega) = \int g(\mathbf{z}, \mathbf{x}, \omega; \mathbf{y}) \mathbf{M} \delta m(\mathbf{y}) \, d\mathbf{y} \quad (19) $$

with

$$ g(\mathbf{x}, \mathbf{z}, \omega; \mathbf{y}) = D(\mathbf{y}, \mathbf{x}) D(\mathbf{z}, \mathbf{y}) \exp(i\omega\Phi_{IT}(\mathbf{y})) \quad (20) $$

$$ \mathbf{M} = \left( -2(-i\omega)^{2-\alpha_0} \tau^{-\alpha_0}, \quad \omega^2/\rho_0(\mathbf{y}') [i\nabla\Phi_{IT}(\mathbf{y}) + \nabla\Phi_{RT}(\mathbf{y})]^2, \quad 2(-i\omega)^{2-\alpha_0} \tau^{-\alpha_0} [\ln(\omega\tau) - i\pi/2] \right) \quad (21) $$

and

$$ \delta m(\mathbf{y}) = \begin{pmatrix} \delta a(\mathbf{y}) \\ \delta \kappa(\mathbf{y}) \\ \delta \alpha(\mathbf{y}) \end{pmatrix}. \quad (22) $$

$\delta\alpha \equiv 0$ if the value of $\alpha$ is kept fixed; in this case the third components of the vector $\mathbf{M}$, $\delta m$, drop out. We denote by $\mathcal{G}$ the integral operator $(\mathcal{G}f)(\mathbf{x}, \mathbf{z}) = \int g(\mathbf{x}, \mathbf{z}; \mathbf{y}) f(\mathbf{y}) d\mathbf{y}$. Denoting by $\mathbf{k}$ the gradient of $\Phi_{IT}$ with respect to the scatterer coordinates we define the Jacobian $I = \partial(\omega, \xi)/\partial(k_1, k_2, k_3) = \partial(\omega, \xi)/\partial(k, \Omega)/(k^2 \sin\theta)$, where $(k, \Omega)$ are the spherical coordinates in the $\mathbf{k}$-space and $\theta$ is the complementary latitude in the $\mathbf{k}$-space.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the recorded trace 23 with the post-inversion synthetic trace 23 obtained by the ray-Born approximation for $\alpha = 0.1$ and for $\alpha = 0.9$, for PVC.</figcaption></figure>
If the acquisition geometry ensures that the mapping $\xi \to \Omega$ is surjective and the mapping $\omega \to k$ is monotone, then a change of the integration variables $(\omega, \xi) \to k$ yields a unit operator (Ribodetti et al. 2000). In addition to checking the invertibility of the mapping $\xi \to \Omega$ it is important to establish the monotonicity and range of the mapping $\omega \to k$.

In a homogeneous background medium the traveltime and the attenuation are proportional to the distance $d$ from the source, $T = d/c$, $A = d/(c\tau)$. Denoting by $D(y)$ the total distance $|x - y| + |z - y|$, we have

$$ \Phi_{IT}(y) = c^{-1}[1 + \sin(\pi\alpha/2)(\tau\omega)^{\alpha}]D(y). $$

Consequently

$$ k = c^{-1}\omega[\omega + \sin(\pi\alpha/2)(\tau\omega)^{\alpha-1}]\nabla D, $$

$\Omega$ does not depend on the frequency $\omega$, the Jacobian can be factored $1/I = (\partial\Omega/\partial\xi) dk/d\omega$, and $k$ is a monotonic function of $\omega$ ranging from 0 to $\infty$.

In an inhomogeneous background medium $k = \omega v + \omega^\alpha w$, and the two vectors $v = \nabla_y[T(x, y) + T(y, z)]$ and $w = \tau^\alpha \sin(\pi\alpha/2) \nabla_y[A(x, y) + A(y, z)]$ are in general not collinear. Suppose that $\gamma$ is the largest real number such that $v \cdot w \ge \gamma vw$ for all $\xi \in U$, where $-1 \le \gamma \le 0$ and $v, w$ denote the lengths of the vectors $v, w$. If $\gamma < -2\sqrt{\alpha}/(1+\alpha)$, then monotonicity breaks down for some low values of $\omega$. Indeed, the derivative on the left-hand side of the inequality

$$ (\omega/2)dk^2/d\omega \ge \omega^{2\alpha}[q^2 + \alpha w^2 + \gamma(1+\alpha)qw]; \quad q := \omega^{1-\alpha}v $$

is non-negative if $\gamma \ge -2\sqrt{\alpha}/(1+\alpha)$.

The operator $Y := G^\dagger W G$ with the weight $W = D[y, x(\xi)]^{-1} D[y', x(\xi)]^{-1} I(\omega, \xi)$ Z becomes

$$ 2\Re \int d\xi \int_0^\infty d\omega \bar{g}(\xi, \omega, y)g(\xi, \omega, y')(R+S) $$

where $M = R + S$, $R$ is the part of $M$ involving the highest powers of $\omega$, while $S$ contains lower powers of $\omega$. If $R$ is invertible and $Z = R^{-1}$, then

$$ Y = 2\Re \int d\xi \int_0^\infty d\omega I \exp[i\omega(\Phi_1(y) - \Phi_1(y'))][U + ZS] $$

where $U$ denotes the unit matrix. Following Beylkin (1985), we apply the approximation $\Phi_{TI}(y) - \Phi_{TI}(y') \approx k \cdot (y - y')$ in the oscillatory factor

<figure><img src="image_10.png" /><figcaption>Figure 9. Comparison of recorded and synthetic traces obtained in three-parameter inversion.</figcaption></figure>
$$Y \approx 2\Re \int d\mathbf{k} \exp[i\mathbf{k} \cdot (\mathbf{y} - \mathbf{y}')][\mathbf{U} + ZS] = 2(2\pi)^3 \mathbf{U}\delta(\mathbf{y} - \mathbf{y}') + X \quad (23)$$

where X is a pseudodifferential operator of a negative order. X maps elements of the model space to smoother objects while the first term of eq. (23) yields a sharp image of the model.

We shall choose R in such a way that it contains all the powers of $\omega$ which tend to $\omega^4$ when $\alpha \to 1$. Expanding $\mathbf{M}^\dagger \mathbf{M}$ we get the Hermitian matrix $\mathbf{R}$,

$$R_{11} = 4\omega^{4-2\alpha_0}\tau^{-2\alpha_0} \quad (24)$$

$$R_{12} = \overline{R_{21}} = -\rho_0^{-1} (\omega^2 \nabla T_T^2 - 2\omega^{1+\alpha_0} e^{i\pi(1+\alpha_0)/2} \nabla T_T \cdot \nabla A_T) \quad (25)$$

$$R_{22} = \rho^{-2} [\omega^4 (\nabla T_T)^4 + 4\omega^{2+2\alpha_0} \tau^{2\alpha_0} (\nabla T_T \cdot \nabla A_T)^2 + 2 \sin(\pi\alpha_0/2) \omega^{3+\alpha_0} \tau_0^\alpha (\nabla T_T)^2 \nabla T_T \cdot \nabla A_T] \quad (26)$$

$$R_{13} = \overline{R_{31}} = -4|\omega|^{4-2\alpha_0} \tau^{-2\alpha_0} [\ln(\omega\tau) - i\pi/2] \quad (27)$$

$$R_{23} = \overline{R_{32}} = 2(-i\omega)^{4-\alpha_0} \tau^{-\alpha_0} (\nabla T_T)^2 [\ln(\tau\omega) - i\pi/2] / \rho \quad (28)$$

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Inversion for the PVC cylinder. Results of the post-processing sequence applied to the tomographic model $\delta a$ to eliminate the source signature and to estimate the absolute values of the recovered parameters. The arrows show the influence of the inaccuracy of the experimental set-up on the tomographic model. (a) Original tomographic model. (b) Best-fitting convolved model. (c) Best-fitting impulsive model representing the geometry and dimensions of the hollow PVC cylinder. (d) Misfit for $\delta a$, the minimum represents the optimal absolute value.</figcaption></figure>
$$R_{33} = 4|\omega|^{4-2\alpha_0}[\ln^2(\tau\omega) + \pi^2/4]\tau^{-2\alpha_0} \quad (29)$$

with $T_T = T(x, y) + T(y, z)$, $A_T = A(x, y) + A(y, z)$, $\nabla = \nabla_y$. The inversion formula has the following form (Ribodetti et al. 2000):

$$\delta m = (G^\dagger W G)^{-1} G^\dagger W \delta G \quad (30)$$

where $\delta G(\xi) = G^{\text{obs}}(z(\xi), x(\xi)) - G_0(z(\xi), x(\xi))$. By construction it is an identity for $G^{\text{obs}}(z(\xi), x(\xi)) = G_B(z(\xi), x(\xi))$.

# 5 APPLICATION TO ULTRASONIC LABORATORY DATA: 1. A HOLLOW PVC CYLINDER

## 5.1 Two-parameter inversion

In two-parameter inversion the value of $\alpha$ is held fixed while the values of $a, \kappa$ are recovered. The value of $\alpha$ is then determined by optimizing the misfit between the recorded seismogram and the synthetic seismograms obtained for the recovered models.

The inversion method described in the previous section was tested on ultrasonic laboratory data collected in a water tank equipped with computer-based control and data acquisition systems (National Instruments system, Valero (1997)). The same experimental data were used in previous work based on the assumption of a constant quality factor $Q$ (Ribodetti et al. 2000). The inversion target was a PVC cylinder immersed in a vertical position in a water tank 2.0 m × 1.4 m × 1.5 m containing 5000 l water. The radius and the height of the PVC cylinder were 0.058 m and 1.30 m, respectively. The axis of the cylinder was positioned at $(X, Y) = (0.0925 \text{ m}; -0.1925 \text{ m})$. The density of PVC is reported as 1.39 g cm⁻³, while the P-wave velocity in PVC is 2380 m s⁻¹ (e.g. in the tables of www.ondacorp.com/tables). The roles of the source and the receiver were played by two hydrophones with the frequency bands 1–110 kHz. The source and the receiver

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Inversion for PVC. Results of the post-processing sequence applied to the tomographic model $\delta\alpha$ in (a). (b) Best-fitting convolved model. (c) Best-fitting impulsive model representing the geometry and dimensions of the hollow PVC cylinder. (d) Misfit for $\delta\alpha$.</figcaption></figure>
were located on a horizontal plane at a depth of 0.67 m beneath the water level and rotated simultaneously around a fixed vertical axis $X = 0, Y = 0$ (Fig. 1). The azimuthal offset between source and receiver was kept constant at 15° during the acquisition. The step $d\theta$ between two consecutive source azimuths was 5°, which results in 72 seismograms per common-offset gather. The dominant frequency of the source signal was 100 kHz, corresponding to a wavelength of 0.015 m in water. Waveforms were digitized with sampling interval of $2 \times 10^{-6}$ s.

In inversion the background medium was assumed to be pure water. The velocity of sound in water was assumed to be $1520 \text{ m s}^{-1}$ and the $Q$ factor was assumed as 210 000, corresponding to the value reported for the frequency of 100 kHz (Fujii & Masui 1993). The target and the acquisition are shown in Fig. 1, the geometry is explained in Fig. 2. Following Ribodetti et al. (2000) we focus on the diffracted wavefield (Fig. 3). The main arrival ($(4-5) \times 10^4 \text{ s}$) is the reflection from the external edge of the cylinder. The second arrival is the reflection from the second cylinder surface. The arrival around $9 \times 10^4 \text{ s}$ is a multiple reflection from inside the cylinder. The flat arrival is a reflection from the bottom of the water tank plus a reflection from the water level. The inversion was performed for several hypothetical values of the exponent $\alpha = 0.1, 0.5, 0.4, 0.6, 0.9$. Comparisons of recorded and synthetics for $\alpha = 0.5, 0.1, 0.4$ are presented in Figs 3, 4 and 5.

The images of the perturbations $\delta a$ and $\delta k$ of $a$ and $k$ recovered after one iteration are shown in Fig. 6. The gross contours of the cylinder both on the attenuation $\delta a$ and the compressibility $\delta k$ are clearly identifiable on both images in Fig. 6. The centre of the cylinder image lies at $(X, Y) = (0.09 \text{ m}; -0.19 \text{ m})$, which is in good agreement with the actual centre of the cylinder. The nearly vertical trace in the $\delta a$ and $\delta k$ sections is due to a missing source-receiver pair during the acquisition of the data. The rms for $\alpha = 0.5, 0.4, 0.6$ are shown in Fig. 7. The discrepancy is markedly lower for $\alpha = 0.5$. Because of the good quality of the fit between observed and synthetics the inversion procedure was limited to a single iteration.

For $\alpha = 0.1$ and $\alpha = 0.9$ the fit is very bad (Fig. 8). For $\alpha = 0.9$ the model is close to the constant-$Q$ model used in the work of (Ribodetti et al. 2000). For these extreme cases we note a deterioration of the fit, particularly noticeable in the phase shift.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Inversion for the PVC cylinder. Results of the post-processing sequence applied to the tomographic model $\delta k$ in (a). (b) Best-fitting convolved model. (c) Best-fitting impulsive model representing the geometry and dimensions of the hollow PVC cylinder. (d) Misfit for $\delta k$.</figcaption></figure>
## 5.2 Three-parameter inversion

In the three-parameter inversion $κ, a$ and $α$ are recovered simultaneously. For a background value $α = 0.5$ the recovered perturbation $δα$ is very small. In Fig. 9 the synthetic seismic traces obtained by the two- and three-parameter inversion are compared. There is hardly a noticeable difference because of the very small variation of the parameter $δα$ for the PVC material, as we will see below.

In order to obtain parameter estimates for the target the deconvolution method described in Ribodetti et al. (2000) was applied. The method is based on a matching of a boxcar family (model space) with a fixed amplitude and width and with a variable position. These boxcars mimic the radius of a cylinder. The amplitude of the boxcar represents the amplitude of the perturbations ($δa, δα, δk$). The width of the boxcar represents the thickness of the cylinder section. To build the predicted tomographic model, the boxcars are converted from space to time using the velocity of the background medium and are convolved with the source wavelet. To estimate the source wavelet, several traces centred on the direct arrival are stacked. In the next step the convolved boxcar is converted back from time to space for comparison with the trace extracted from the tomographic image. The procedure is performed by a systematic exploration of the model space for each azimuth choosing the best-fitting model in a least-squares sense. Once all the azimuths are explored, the best-fitting is convolved and impulsive boxcars are put together in the (X, Y) plane to build 2-D synthetic images of the cylinder section. The results of this post-processing are shown in Figs 10, 11 and 12 for $δa, δα, δk$ perturbations.

The shape of the recovered cylinder section on the $δa, δα, δk$ synthetic images (Figs 10c, 11c and 12c) is acceptable. It is more difficult to recognize the shape from the $δa, δα$ images related to the viscoelastic attenuation. The recovered shape of the cylinder in all the images is slightly oval. This may result from a combination of inaccuracies in the experimental set-up (denoted by arrows on the models) rather than from the processing itself. The post-processing failed locally to fit the exact shape due to the sensitivity of the parameters to the amplitude. To estimate the absolute value of the perturbations in the cylinder, we applied the post-processing for several values of $δa, δα, δk$ perturbations (namely, for several amplitudes of the boxcars). We computed the $l_2$ misfit for each tested amplitude. The misfit as a function of the perturbation is displayed in Figs 10(d), 11(d) and 12(d). A minimum is well identified for all the parameters: for $δa$ the minimum lies around $0.8 \times 10^{-4}$ corresponding to a perturbation $δQ \approx -50000$ of the $Q$ factor with respect to the reference $Q$ value, equal to the experimental value at the central frequency of 100 kHz. This result is in agreement with previous results obtained in Ribodetti et al. (2000), since the PVC presents essentially an elastic response with a very high $Q$. The perturbation $δα$ has a minimum around 0.05, which indicates a small variation of the water around the reference value $α = 0.5$. An exact value $α = 0.5$ is possibly attributable to purely thermal dissipation of the water (Kelbert & Chaban 1986). $δk$ has a minimum of around $5.2 \times 10^{-6} \text{ N m}^{-2}$ corresponding to a velocity reported at www.ondacorp.com/tables.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison of recorded data (in red) and synthetic seismograms (in blue) for the lava specimen.</figcaption></figure>
## 6 APPLICATION TO ULTRASONIC LABORATORY DATA: 2. A CYLINDRICAL LAVA SPECIMEN

A cylindrical Neapolitan Yellow Tuff sample (pyroclastic rock) from a drilling experiment performed near the Chiaiano (Campi Flegrei, Italy) of height 0.5 m and diameter 0.06 m was immersed in the centre of the water tank described in the previous section. The source and receiver were identical were the same as in the previous section. They were placed on a horizontal plane orthogonal to the sample at a depth of 0.67 m below the water level and rotated synchronously around the axis of the cylinder at a distance of 0.469 m. The azimuthal offset between the source and the receiver was again 15°, the initial azimuth of the source was 7.55° and the azimuth step was 55°. The central frequency of the source is around 100 kHz. The diffracted wavefield is presented in Fig. 13. The main arrivals between 0.5 and 0.6 ms are respectively the reflections from the external and internal edge of the rock sample. The ray-Born synthetics obtained at the first iteration using the perturbation in Figs 14, 15 and 16, are superimposed with the observed traces. A good fit is obtained. The tomographic models for $\delta k$, $\delta a$ and $\delta \alpha$ are shown in Figs 14(a), 15(a) and 16(a). The best-fitting convolved model shown in parts (b) is obtained by convolution with the source wavelet shown in

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Inversion for the lava sample. (a) Tomographic model for $\delta k$. (b) The best-fitting convolved model results of the post-processing sequence applied to (a). The source wavelet is plotted in (c). (d) The best-fitting impulsive model showing the geometry and location of contour of the section of the rock sample. (e) Misfit function representing the perturbations of $\delta k$ and the minimum.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Inversion for the lava sample. (a) Tomographic model for $\delta a$. (b) The best-fitting convolved model from the post-processing sequence applied to (a). (c) The best-fitting impulsive model. (d) The misfit function representing the perturbations of $\delta a$ and the minimum.</figcaption></figure>

Fig. 14(c), while the model shown in parts (d) is obtained by convolution with the best-fitting impulsive model. The influence of an imperfect acquisition geometry is indicated by arrows. Because of the ray-Born approximation and the constant background model assumed during inversion, only the physical properties along the contour of the rock sample can be estimated. Concerning the geometry, a good location of the sample and its dimensions are obtained successfully (see Figs 14c, 15c and 16c). The shape of the recovered section is acceptable on all three images. Nevertheless, we note that the recovered shape of the cylinder is slightly oval. This may result from a combination of inaccuracies in the experimental set-up.

On the $\delta a$ and $\delta\alpha$ synthetic images the post-processing failed locally to fit the exact shape due to the sensitivity of the parameters to the amplitude variations related to some local high heterogeneity of the sample.

The absolute quantitative estimation of the propagation parameter $\delta k$ and attenuation parameters $\delta a$ and $\delta\alpha$ are obtained during the post-processing of the tomographic models. The misfit in the $l_2$ sense between the tomographic model and the best-fitting convolved model is also determined. The misfit for $\delta k$ presented in Fig. 14(e) shows a minimum at $3.908 \times 10^6$ N m$^{-2}$, corresponding to a velocity perturbation $\delta v = 1318$ m s$^{-1}$, then the velocity of the tuff sample is $v_{tuff} = v_0 + \delta v \approx 3800$ m s$^{-1}$; for $\delta a$ the misfit presents a minimum at $1.689 \times 10^{-4}$, corresponding to a $\delta Q$ around $-209880$ estimated at the central frequency of 100 kHz. The attenuation $Q$ of the sample is $Q_{tuff} = Q_0 + \delta Q \approx 120$. The estimates of the velocity and the $Q$ factor are in agreement with results obtained by Vanorio & Nur (2000). For the parameter $\delta\alpha$ (Fig. 16d) the minimum of the misfit function is obtained at $\delta\alpha = 0.15$ corresponding to $\alpha_{tuff} = \alpha_0 + \delta\alpha = 0.5 + 0.15 = 0.65$. The estimated value of $\alpha$ is attributable to wave attenuation in a porous water-saturated sample. A similar value was obtained by Hanyga & Rok (2000) from the Gurevich-Lopatnikov model of fast wave attenuation in a thin-layered porous medium, but there is not enough evidence for validating this model for the lava specimen in our experiment.

## 7 CONCLUSIONS

We have demonstrated that the viscoelastic inversion developed by Ribodetti et al. (2000) yields fairly sharp and reliable parameter estimates in the power-law model class. The power model class is more realistic than the constant-$Q$ model class because dispersion and absorption
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Inversion for the lava sample. (a) Tomographic model for $\delta\alpha$. (b) The best-fitting convolved model results of the post-processing sequence applied to (a). (c) The best-fitting impulsive model. (d) The misfit function representing the perturbations of $\delta\alpha$ and the minimum.</figcaption></figure>

effects are simultaneously taken into account and the causality principle is respected. For ultrasonic laboratory data the optimal value of the exponent of the power law for wave propagation in water is almost exactly 0.5, which can be attributed to thermal rather than viscous dissipation.

Our results indicate that the model parametrization is physically meaningful, in contrast to the commonly used approximation of viscous memory by a superposition of exponentials (Day & Minster 1984; Emmerich & Korn 1987). The latter approach requires an arbitrary choice of relaxation times and therefore cannot yield a unique result in seismic inversion. Furthermore, the value of the power-law exponent can be linked to a specific physical mechanism, such as viscous or thermal boundary layers (Kelbert & Chaban 1986; Kelbert & Sazonov 1996).

Low-frequency behaviour of materials is incorrectly represented by the power-law model, since according to the power-law model the wave speed vanishes at zero frequency and stress relaxes to zero after an infinite time. This might not be relevant for wave propagation, but the Cole–Cole model (Hanyga 2003) and its generalizations (Rossikhin & Shitikova 2001) can be used to correct this shortcoming. The Cole–Cole model allows for a different power law in the low- and high-frequency regions. More complicated models can accommodate several transitions in the frequency region.

The main implication of the theoretical analysis of the inversion algorithm is that when the power-law model is used in the linear inversion schema, the parameter related to the attenuation factor and the parameter related to the propagation are decoupled for the acquisition geometries used in laboratory and in field seismic experiments (the Hessian matrix is not singular).

The essential practical result is that when the power of the model $\alpha$ is 0.5 (in the case of water as a reference medium) the viscoelastic inversion runs successfully on ultrasonic laboratory data. In many respects the results of inversion in power-law media and in the constant-$Q$ class are in agreement. The choice of a power law-medium in the inversion procedure, however, yields a good fit between observed and synthetics in only one iteration. This fact is important because of its implication for the reduction of the computation time.

Our results of inversion of ultrasonic laboratory data for a hollow PVC cylinder and for a lava sample show that the adopted method can be useful for estimating rock properties in laboratory experiments; in particular the parameter $\alpha$ governs the frequency dependence of the relaxation functions and can be an interesting indicator of the degree of porosity of the material.

The CPU time on one processor of a SUN E450 computer required for one iteration was 45 min for the three-parameter inversion.
ACKNOWLEDGMENTS

The research was made possible by a grant of PAI and the Norwegian Royal Society in the framework of the programme Aurore of the Franco-Norwegian scientific exchanges 2002–2003. The authors thank H. P. Valero, G. Saracco, F. Conil and D. Gibert for the ultrasonic data collected in the water tank at the University of Rennes, France. The authors also acknowledge stimulating comments of one reviewer on the mathematical challenges of a rigorous definition of the operators applied in the inversion. Publication number 648 of Géosciences-Azur.

REFERENCES

Beylkin, G., 1985. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform, *J. Math. Phys.*, **26**, 99–108.

Dautray, R. & Lions, J.-L., 1992. *Mathematical Analysis and Numerical Methods for Science and Technology*, Vol. 5, Springer, Berlin.

Day, S.M. & Minster, J.B., 1984. Numerical simulation of wavefields using a Padé approximant method, *Geophys. J. R. astr. Soc.*, **78**, 105–118.

Emmerich, M. & Korn, M., 1987. Incorporation of attenuation into time-domain computation of seismic wavefields, *Geophysics*, **52**, 1252–1264.

Engler, H., 1997. Similarity solutions for a class of hyperbolic integro-differential equations, *Differential Integral Eqns*, **10**, 815–845.

Fujii, K. & Masui, R., 1993. Accurate measurements of the sound velocity in pure water by combining a coherent phase-detection technique and a variable path-length interferometer, *J. Acoust. Soc. Am.*, **93**, 273–282.

Hanyga, A., 2003. An anisotropic Cole–Cole viscoelastic model of seismic attenuation: well-posedness and numerical methods, *J. Comput. Acoust.*, **11**, 75–90.

Hanyga, A. & Rok, V.E., 2000. Wave propagation in micro-heterogeneous porous media: a model based on an integro-differential equation, *J. acoust. Soc. Am.*, **107**, 2965–2972.

Hanyga, A. & Seredyńska, M., 1999a. Some effects of the memory kernel singularity on wave propagation and inversion in poroelastic media, I: forward modeling, *Geophys. J. Int.*, **137**, 319–335.

Hanyga, A. & Seredyńska, M., 1999b. Asymptotic ray theory in poro- and viscoelastic media, *Wave Motion*, **30**, 175–195.

Hanyga, A. & Seredyńska, M., 2002. Asymptotic wave front expansions in hereditary media with singular memory kernels, *Q. Appl. Math.*, **LX**, 213–244.

Hanyga, A. & Seredyńska, M., 2003. Power law attenuation in acoustic and isotropic anelastic media, *Geophys. J. Int.*, **155**, 830–838.

Kelbert, M.Y. & Chaban, I.Y., 1986. Relaxation and propagation of pulses in fluids, *Izv. Ak. Nauk, Ser. Mech. Fluids Gases*, **5**, 153–160.

Kelbert, M. & Sazonov, I., 1996. *Pulses and Other Wave Processes in Fluids: an Asymptotical Approach to Initial Problems*, Kluwer Academic, Dordrecht.

Kreis, A. & Pipkin, A.C., 1986. Viscoelastic pulse propagation and stable probability distributions, *Q. Appl. Math*, **44**, 353–360.

Lokshin, A.A. & Rok, V.E., 1978a. Fundamental solutions of the wave equation with delayed time, *Doklady Akad. Nauk SSSR*, **239**, 1305–1308.

Lokshin, A.A. & Rok, V.E., 1978b. Automodel solutions of wave equations with time lag, *Russ. Math. Surveys*, **33**, 243–244.

Lokshin, A.A. & Suvorova, Y.V., 1982. *Mathematical Theory of Wave Propagation in Media with Memory*, Moscow University, Moscow (in Russian).

Lu, J. & Hanyga, A., 2004. Numerical modeling of wave propagation in a linear viscoelastic medium with singular memory, Submitted to *Geophys. J. Int.*

Podlubny, I., 1998. *Fractional Differential Equations*, Academic Press, San Diego.

Ribodetti, A., Operto, S., Virieux, J., Lambaré, G., Valero, H.P. & Gibert, D., 2000. Asymptotic viscoacoustic diffraction tomography of ultrasonic laboratory data: a tool for rock properties analysis, *Geophys. J. Int.*, **140**, 324–340.

Rossikhin, Y.A. & Shitikova, M.V., 2001. Analysis of dynamic behaviour of viscoelastic rods whose rheological models contain fractional derivatives of two different orders, *Z Angew. Math. Mech.*, **6**, 363–378.

Szabo, T.L., 1994. Time domain wave equations for lossy media obeying a frequency power law, *J. acoust. Soc. Am.*, **96**, 491–500.

Szabo, T.L., 1995. Causal theories and data for acoustic attenuation obeying a frequency power law, *J. acoust. Soc. Am.*, **97**, 14–24.

Uchaikin, V.U. & Zolotarev, V.M., 1999. *Chance and Stability: Stable Distributions and their Applications*, VSP, Utrecht.

Valero, H.P., 1997. Endoscopie sismique, *PhD thesis*, Institut de Physique du Globe de Paris, Paris.

Vanorio, T., Prasad, M., Patella, D. & Nur, A., 2000. Asymptotic viscoacoustic diffraction tomography of ultrasonic laboratory data: a tool for rock properties analysis, *Geophys. J. Int.*, **140**, 324–340.