## Abstract

An asymptotic Green's function in homogeneous anisotropic viscoelastic media is derived. The Green's function in viscoelastic media is formally similar to that in elastic media, but its computation is more involved. The stationary slowness vector is, in general, complex valued and inhomogeneous. Its computation involves finding two independent real-valued unit vectors which specify the directions of its real and imaginary parts and can be done either by iterations or by solving a system of coupled polynomial equations. When the stationary slowness direction is found, all quantities standing in the Green's function such as the slowness vector, polarization vector, phase and energy velocities and principal curvatures of the slowness surface can readily be calculated.

The formulae for the exact and asymptotic Green's functions are numerically checked against closed-form solutions for isotropic and simple anisotropic, elastic and viscoelastic models. The calculations confirm that the formulae and developed numerical codes are correct. The computation of the *P*-wave Green's function in two realistic materials with a rather strong anisotropy and absorption indicates that the asymptotic Green's function is accurate at distances greater than several wavelengths from the source. The error in the modulus reaches at most 4% at distances greater than 15 wavelengths from the source.

## 1. Introduction

The theory of waves propagating in viscoelastic media has become increasingly important and intensively developed in recent years (Auld 1973; Caviglia & Morro 1992; Carcione 2001). This concerns studying the basic attributes of homogeneous and inhomogeneous plane waves such as slowness and attenuation vectors (Červený & Pšenčík 2005), deriving formulae for linearized attenuation coefficients (Zhu & Tsvankin 2006, 2007), studying acoustic axes (Shuvalov & Scott 1999; Shuvalov 2001), energy flux (Carcione & Cavallini 1993; Deschamps *et al.* 1997; Boulanger 1998; Červený & Pšenčík 2006) and reflection/transmission coefficients at a plane interface between two homogeneous viscoelastic halfspaces (Borcherdt 1977, 1982; Krebes 1983; Richards 1984; Wennerberg 1985; Nechtschein & Hron 1996; Brokešová 2001; Daley & Krebes 2004). In contrast to plane-wave propagation, considerably less attention has been paid to the theory of waves excited by point or finite sources in isotropic and anisotropic viscoelastic media (Carcione 1994, 2001). So far, such wave fields are usually computed using numerical methods that directly solve the wave equation (Carcione 1990; Saenger & Bohlen 2004; Moczo *et al.* 2007).

In this paper, I will study the waves generated by point sources. I generalize the formula for the exact Green's function in homogeneous anisotropic elastic media (Buchwald 1959; Burridge 1967; Willis 1980; Norris 1994; Wang & Achenbach 1994) to viscoelastic media and present closed-form solutions of the exact viscoelastic Green's function for isotropy and three simple types of anisotropy. Further, I derive the asymptotic Green's function valid for high-frequency signals at distances far from the source and define wave fronts and other asymptotic quantities such as ray velocity and ray attenuation. In numerical modelling, I verify the correctness of the analytical, exact numerical and asymptotic Green's functions, and demonstrate their accuracy and efficiency.

In the following, real and imaginary parts of complex-valued quantities are denoted by superscripts R and I, respectively. A complex-conjugate quantity is denoted by an asterisk. The vector direction of a complex-valued vector ** v** is calculated as (the normalization condition is not used because it complicates generalizing some of the real-valued equations to complex-valued ones). If any complex-valued vector is defined by a real-valued direction, it is called homogeneous, and if defined by a complex-valued direction, it is called inhomogeneous. In formulae, the Einstein summation convention is used for repeated subscripts.

## 2. Anisotropic viscoelastic medium

Let us assume a time-harmonic plane wave described by the displacement vector(2.1)where ** x** is the position vector;

*U*is the amplitude;

**is the unit polarization vector;**

*g**ω*is the circular frequency; and

*t*is the time. Vector

**is the slowness vector defined as**

*p***=**

*p***/**

*n**c*, where

**is the slowness direction and**

*n**c*is the phase velocity. The plane wave propagates in the viscoelastic medium defined by the following stress–strain relation:(2.2)where

*σ*

_{ij}are the components of the stress tensor;

*e*

_{kl}are the components of the strain tensor; is the time derivative of

*e*

_{kl}; and and

*η*

_{ijkl}are the components of the elasticity and viscosity tensors, respectively. These tensors must be positive definite to ensure that strain energy

*W*and dissipation energy rate

*R*are positive,(2.3)Taking into account that is expressed for the time-harmonic wave (2.1) as(2.4)the stress–strain relation (2.2) can be formally written in the frequency domain in a form identical with that for elastic media, where the real-valued elasticity tensor is substituted by the complex-valued frequency-dependent viscoelasticity tensor ,(2.5)This tensor is used to construct the complex-valued density-normalized viscoelasticity tensor , and subsequently the complex-valued Christoffel tensor

*Γ*

_{jk}, defined alternatively in terms of slowness direction

**,(2.6)or slowness vector**

*n***,(2.7)Slowness direction**

*p***is real valued for homogeneous waves, but complex valued for inhomogeneous waves.**

*n*The Christoffel tensor *Γ*_{jk} has three eigenvalues and three eigenvectors, which define the properties of three plane waves propagating in anisotropic media. Eigenvalues *G*(** n**) and

*G*(

**) read(2.8)(2.9)and define phase velocity**

*p**c*and slowness vector

**as a function of slowness direction**

*p***. The eigenvectors define polarization vectors**

*n***. The polarization vectors are normalized so that**

*g***.**

*g***=1. Note that when applying condition**

*g***.**

*g*

*g*^{*}=1, equations (2.8) and (2.9) would have a different form also containing complex-conjugate quantities.

Using the eigenvalue *G*(** p**), we define the energy velocity vector as(2.10)which is called the group velocity vector in elastic media. Vectors

**and**

*v***are related by the equation**

*p***.**

*v***=1. The slowness vector, phase velocity, energy velocity and the polarization vectors are in general complex valued.**

*p*## 3. Exact Green's function

### (a) Elastic medium

The exact elastodynamic Green's function in unbounded, homogeneous, anisotropic, perfectly elastic media can be expressed in the frequency domain as the sum of regular and singular terms, and , as follows (Norris 1994, eqn. (3.22); Wang & Achenbach 1994, eqns (17)–(21)):(3.1)(3.2)(3.3)where the superscript *M*=1, 2 and 3 denotes the type of wave (*P*, *S*1 and *S*2); ** g**=

**(**

*g***) is the unit polarization vector;**

*n**c*=

*c*(

**) is the phase velocity;**

*n**ρ*is the density of the medium;

*ω*is the circular frequency;

**=**

*x*

*N**r*is the position vector;

*r*is the distance of the observation point from the source;

**is the ray vector;**

*N***is the slowness direction; and**

*n**S*(

**) is the unit sphere. The regular term is integrated over a hemisphere, defined by**

*n***.**

*n***>0. The singular term is integrated over a unit circle defined by**

*x***.**

*n***=0. The reduction of the surface integral to a line integral in the singular term is due to the Dirac delta function**

*x**δ*(

**.**

*n***) in the integrand. Since the singular term does not depend on frequency**

*x**ω*, it physically corresponds to the static Green's function.

All quantities in (3.1)–(3.3) are real valued, including slowness direction ** n**. Hence, the regular part of the exact Green's function is calculated as a superposition of homogeneous plane waves. This is the essential advantage of formulae (3.1)–(3.3) compared with other formulae for the Green's function which also involve inhomogeneous waves specified by a complex-valued slowness direction

**(e.g. the Green's function in isotropic media calculated using the Weyl integral; see Aki & Richards 2002, §6.1).**

*n*Formulae (3.1)–(3.3) can be derived in several alternative ways. Probably, the most straightforward way is to solve the elastodynamic equation in the slowness-frequency domain to obtain . Subsequently, is calculated by applying the inverse three-dimensional Fourier transform (see Mura 1987; Norris 1994; Červený 2001), which physically corresponds to a superposition of homogeneous plane waves propagating in an arbitrary direction with an arbitrary phase velocity. Using the residuum theorem, the volume integral is further reduced to a surface integral representing a superposition of homogeneous plane waves propagating in an arbitrary direction, but with the phase velocity satisfying equation (2.8).

The integration in (3.2) and (3.3) can also be performed over slowness surface *S*(** p**). Taking into account the relation between surface elements d

*S*(

**) and d**

*n**S*(

**) (see Burridge 1967),(3.4)and the relation between the phase and energy velocities**

*p**c*and

*v*,(3.5)we can put for a particular wave (

*P*,

*S*1 or

*S*2)(3.6)(3.7)The complete Green's function is then a sum of contributions of all three waves.

### (b) Viscoelastic medium

The Green's function in viscoelastic media is derived in a quite analogous way as in elastic media. The Green's function is calculated by applying the inverse three-dimensional Fourier transform, which physically corresponds to a superposition of homogeneous plane waves propagating in an arbitrary direction with an arbitrary (real-valued) phase velocity. Using the residuum theorem, the volume integral is reduced to a surface integral representing a superposition of homogeneous plane waves propagating in an arbitrary direction, but with complex-valued velocities satisfying equation (2.8) and with polarization vectors satisfying the Christoffel equation. Hence, the Green's function is expressed by exactly the same formulae as formulae (3.1)–(3.3), (3.6) and (3.7) derived for elastic media. The only difference is that some quantities in the integrands become complex valued and frequency dependent. This concerns phase velocity *c*, slowness vector ** p**, energy velocity

*v*and polarization vector

**. Slowness direction**

*g***, ray direction**

*n***, position vector**

*N***and density remain real valued and do not depend on real-valued frequency**

*x**ω*. Since some quantities in (3.3) and (3.7) are frequency dependent, the singular term of the Green's function becomes also frequency dependent

Emphasize that ** p** is a homogeneous complex-valued vector; hence its direction

**in (3.1)–(3.3) remains real valued. This is in agreement with the above statement that the exact Green's function in viscoelastic media is calculated as a superposition of homogeneous plane waves, similarly as in elastic media.**

*n*## 4. Asymptotic Green's function

Since the exact Green's function is expressed in terms of complicated integrals, it is advantageous to evaluate the Green's function asymptotically. Physically, it corresponds to studying the high-frequency Green's function in the far field, i.e. at distances far from the source with respect to the predominant wavelength of the wave field. The asymptotic Green's function significantly simplifies, because the singular term (3.7) is negligible for high frequencies, and the regular term (3.6) can be calculated using the stationary phase method for the elastic medium or the steepest descent method for the viscoelastic medium.

### (a) Elastic medium

Formula (3.6) can be expressed as(4.1)where(4.2)Functions and *θ*(** p**) are assumed to be smooth, and the parameter

*λ*=

*ωr*is large and positive. Phase

*θ*(

**) is positive inside the integration area and zero at its boundary. Expanding phase function**

*p**θ*(

**) near stationary point**

*p*

*p*_{0}, defined by(4.3)yields(4.4)where

*K*

_{1}and

*K*

_{2}are the principal curvatures of the slowness surface at

*p*_{0}, and

*s*

_{1}and

*s*

_{2}are the local coordinates defined on the slowness surface. The coordinate axes have their origin at

*p*_{0}and point along the principal curvatures of the slowness surface at

*p*_{0}. The third coordinate axis

*s*

_{3}is normal to the slowness surface at

*p*_{0}and parallel to

**.**

*N*Applying the stationary phase method to integral (4.1), we obtain the asymptotic Green's function in the following form (Buchwald 1959; Burridge 1967; Every & Kim 1994; Červený 2001, eqn. (2.5.75)):(4.5)where is the Gaussian curvature of the slowness surface, and *σ*_{0} is defined as(4.6)Both the principal curvatures *K*_{1} and *K*_{2} are positive for a convex surface, negative for a concave surface, and of opposite signs for a saddle-shaped surface. All quantities dependent on ** p** in (4.5) are taken at stationary point

*p*_{0}. Formula (4.5) works for all directions except for singularities on the slowness surface and cusp edges on the wave front, where a more involved approach is required (Vavryčuk 1999, 2002; Gridin 2000).

### (b) Viscoelastic medium

Since phase velocity *c* is complex valued in viscoelastic media, the slowness vector ** p** and functions and

*θ*(

**) in (4.1) are also complex valued. Parameter**

*p**λ*is real, large and positive. The integration in (4.1) is over complex-valued slowness surface

*S*(

**). The stationary point**

*p*

*p*_{0}on

*S*(

**) is defined by(4.7)where**

*p**s*

_{1}and

*s*

_{2}are the complex-valued local coordinates defined on

*S*(

**). The third complex-valued coordinate axis**

*p**s*

_{3}is normal to the slowness surface at

*p*_{0}and parallel to

**. Since**

*N***is a real-valued vector, the**

*N**s*

_{3}-axis is homogeneous. Consequently, the energy velocity vector,(4.8)must also be homogeneous. Expanding phase function

*θ*(

**) near stationary point**

*p*

*p*_{0}yields(4.9)where

*K*

_{1}and

*K*

_{2}are the complex-valued principal curvatures of the slowness surface at

*p*_{0}. Applying the steepest descent method (Ben-Menahem & Singh 1981, appendix E) to integral (4.1), we obtain(4.10)and finally(4.11)whereAngles and define the phases of the principal curvatures

*K*

_{1}and

*K*

_{2}, and is the modulus of the Gaussian curvature . All quantities dependent on

**in (4.11) are taken at stationary point**

*p*

*p*_{0}. Position vector

**, distance**

*x**r*, frequency

*ω*, phase angles and and density

*ρ*are real valued, but polarization vector

**, Gaussian curvature**

*g**K*, principal curvatures

*K*

_{1}and

*K*

_{2}, energy velocity

*v*and slowness vector

*p*_{0}are complex valued. Similarly as for elastic media, formula (4.11) fails in singularities.

Decomposing the energy velocity *v* into real and imaginary parts, *v*^{R} and *v*^{I}, the exponential term in (4.11) becomes(4.12)where(4.13)(4.14)The asterisk in (4.13) and (4.14) indicates a complex-conjugate quantity. Consequently, the Green's function reads(4.15)Since and control the propagation velocity and attenuation along a ray, I will refer to them as the ray velocity and ray attenuation, respectively.

## 5. Numerical procedure

### (a) Stationary point

Determining stationary point *p*_{0} is the crucial and the most complicated step in calculating the asymptotic quantities in anisotropic viscoelastic media. Although the exact Green's function (3.1)–(3.3) is calculated as a superposition of homogeneous plane waves, the asymptotic Green's function is calculated at a stationary point which corresponds to an inhomogeneous plane wave. Hence, the slowness direction at the stationary point is no longer real valued, but it is generally complex valued.

The complex-valued stationary point *p*_{0} can be determined either by iterations or by solving a system of polynomial equations. The iterative procedure is fast and works well provided the wave front is free of triplications. The procedure is based on seeking a complex-valued slowness direction *n*_{0}, for which energy velocity ** v** is homogeneous and points along a specified ray vector

**. This represents an inversion for four unknown angles: two angles define the real part and two others define the imaginary part of**

*N*

*n*_{0}. We can consider the modulus of the complex-valued deviation between the given and predicted ray vectors as the misfit function. The inversion is nonlinear and can be performed using standard methods (Press

*et al.*2002). If the wave front displays triplications, determination of stationary slowness

*p*_{0}is more involved. When using iterations, we have to be careful to find all slowness vectors corresponding to a given ray, which may sometimes be tricky. The other possibility is to follow Vavryčuk (2006

*a*) and solve a system of polynomial equations for the unknown components of

*p*_{0}.

### (b) Principal curvatures of the slowness surface

The principal curvatures of the slowness surface *K*_{1} and *K*_{2} can conveniently be calculated using the wave metric tensor (Vavryčuk 2003, eqn. (9)),(5.1)which is expressed for the *P* wave as follows (Vavryčuk 2003, eqn. (17)):(5.2)(5.3)(5.4)where superscripts *M* and *N* denote the type of wave (*P*, *S*1 or *S*2). The wave metric tensors of the *S*1 or *S*2 waves are analogous (see Vavryčuk 2003, eqns (18) and (19)). The diagonal form of reads (Klimeš 2002)(5.5)hence principal curvatures *K*_{1} and *K*_{2} can directly be calculated from the eigenvalues of . Subsequently, Gaussian curvature *K* is calculated as (Klimeš 2002, eqn. (45))(5.6)The formulae for the principal and Gaussian curvatures work safely in all directions except for singularities, in which or . In such directions and in their close vicinity, equation (5.2) fails.

## 6. Examples

In this section, the derived formulae are validated on numerical examples. First, the formulae are checked against simple closed-form solutions. Second, the formulae are applied to exemplify the properties of the Green's function in realistic anisotropic viscoelastic materials.

### (a) Test models

For isotropic and some simple anisotropic media, the Green's function can be expressed in closed form. The existing closed-form solutions have been derived for elastic media, but they can simply be generalized to be applicable also to viscoelastic media (see appendices A and B). Subsequently, the viscoelastic closed-form Green's functions can be used to verify the derived formulae and to test numerical codes.

I adopted two viscoelastic and two elastic test models. The first viscoelastic model is isotropic with parameters in the Voigt notation: , , and with density ; and the second viscoelastic model is transversely isotropic with parameters: , , , , , and with density . The elastic models are described by real parts of the above viscoelastic parameters. The circular frequency *ω* of the wave field equals 1.

For the four models, I calculated the exact Green's function by direct numerical calculation of integrals (3.2) and (3.3), an asymptotic Green's function (4.11) and the exact and far-field analytical Green's functions using (A 1) and (B 3). Figures 1 and 2 show the Green's functions in both the elastic and viscoelastic models. The exact solution obtained using the numerical integration coincides with the analytical solution within the width of the line. Similarly, the asymptotic solution perfectly coincides with the analytic far-field approximation. This validates the formulae derived as well as the codes developed.

### (b) Real models

The Green's function is calculated in two anisotropic viscoelastic materials: clay shale and carbon-epoxy composite (see table 1). The materials display anisotropy in velocities as well as in attenuation. The values for the clay shale were taken from Carcione & Cavallini (1995) and those for the carbon-epoxy composite from Hosten *et al.* (1987). The values were slightly modified to make the models transversely isotropic (see table 2).

Figure 3 shows the basic properties of the *P* waves propagating in the materials studied: the left-hand plots show the directional dependence of ray velocity and real-valued phase velocity , and the right-hand plots show the directional dependence of ray attenuation and of phase attenuation . The ray velocity and ray attenuation characterize the propagation velocity and attenuation along a ray, while the phase velocity and phase attenuation characterize the propagation velocity and attenuation along the normal to the plane wave front. The magnitudes of the phase quantities are calculated from the complex-valued slowness vector ** p** as(6.1)and their direction is along the normal to the plane wave front,(6.2)The magnitudes of the ray quantities are defined by (4.13) and (4.14) and they point along ray vector

**. Emphasize that slowness vector**

*N***is stationary for both the ray and phase quantities; hence it is, in general, an inhomogeneous complex-valued vector. The corresponding energy velocity vector**

*p***must be homogeneous, and ray vector is real valued.**

*v*Figure 4 shows a comparison of the exact and asymptotic *P*-wave Green's functions for the two materials studied. The figure indicates that the asymptotic Green's function approximates the exact Green's function quite efficiently. As expected, the fit improves with increasing distance from the source. At distance of 15 wavelengths from the source, the error in the modulus of the asymptotic Green's function is approximately 4% or less for all four examples in figure 4.

## 7. Conclusion

The exact Green's function in homogeneous anisotropic viscoelastic media is calculated as in elastic media by numerically integrating surface and line integrals. The surface integral defines the regular part of the Green's function and physically corresponds to a superposition of homogeneous plane waves. The line integral defines the singular part of the Green's function. While all quantities in the integrals are real valued in elastic media, some of them become complex valued and frequency dependent in viscoelastic media. However, this difference is minor and does not pose complications to computing the Green's function. The complex-valued integrals can be calculated either directly using numerical integration or asymptotically using the steepest descent method. In some very simple cases, the integrals can also be calculated analytically yielding closed-form solutions.

The calculation of the asymptotic Green's function is more involved in viscoelastic media than in elastic media. Although the exact Green's function in viscoelastic media is calculated as a superposition of homogeneous plane waves, the asymptotic Green's function is calculated at a stationary point which corresponds to an inhomogeneous plane wave. In elastic media, we usually start with setting a slowness vector and then calculate the corresponding ray direction. The slowness vector is always homogeneous. In viscoelastic media, we first define the ray direction and seek the slowness vector that predicts the vector of the energy velocity parallel to the ray. Since the ray direction is real, the energy velocity vector is homogeneous. Consequently, the stationary slowness vector is, in general, inhomogeneous and its computation involves finding two independent unit vectors, which specify the directions of its real and imaginary parts. This can be done iteratively in a simple way, provided the wave surface is free of triplications. If not, we have to take care of getting all slowness vectors corresponding to a given ray, or we have to solve a system of polynomial equations.

## Acknowledgments

I thank Colin Thomson for discussions on the subject, and Klaus Helbig and two anonymous referees for their comments. This work was supported by the grant Agency of the Czech Republic, grant no. 205/05/2182 and by the Consortium Project SW3D ‘Seismic Waves in Complex 3-D Structures’.

## Footnotes

- Received March 7, 2007.
- Accepted May 4, 2007.

- © 2007 The Royal Society