## Abstract

Efforts at modelling the propagation of seismic waves in half-spaces with continuously varying properties have mostly been focused on shear-horizontal waves. Here, a sagittally polarized (Rayleigh type) wave travels along a symmetry axis (and is attenuated along another) of an orthotropic material with stiffnesses and mass density varying in the same exponential manner with depth. In contrast to what could be expected at first sight, the analysis is very similar to that of the homogeneous half-space, with the main and capital difference that the Rayleigh wave is now dispersive. The results are illustrated numerically for (i) an orthotropic half-space typical of horizontally layered and vertically fractured shales and (ii) for an isotropic half-space made of silica. In both the examples, the wave travels at a slower speed and penetrates deeper than in the homogeneous case. In the second example, the inhomogeneity can force the wave amplitude to oscillate as well as decay with depth, in marked contrast with the homogeneous isotropic general case.

## 1. Introduction

Love (1911) showed that an *inhomogeneous* half-space, consisting of an elastic layer covering a semi-infinite body made of a different elastic material, can sustain the propagation of a linearly polarized (shear horizontal) surface wave. The Love wave is faster than the elliptically polarized (vertical) Rayleigh (1885) wave and it has been recorded countless times during earthquakes or underground explosions. Another recorded phenomenon is that Rayleigh waves are *dispersive*, a characteristic which is incompatible with the context of a homogeneous half-space given by Rayleigh (1885). Love showed that his layer/substrate configuration could also support a two-partial, vertically polarized, surface wave. Since this configuration introduces a new characteristic length, the layer thickness *h* (say), a dispersion parameter is now *kh* where *k* is the wavenumber, and that surface wave is dispersive.

Subsequent analyses introduced more and more layers to refine the model, until it was considered practical to view the inhomogeneity of the half-space as a *continuous* variation of the material properties (Ewing *et al*. 1957). The chief among these continuous variations is the one for which the elastic stiffnesses and the mass density vary exponentially with depth, all in the same manner, proportional to a common factor exp(−2*αx*_{2}) say, where *α* is the inverse of a inhomogeneity characteristic length, and *x*_{2} is the coordinate along the normal to the free surface, so that here a dispersion parameter is now *α*/*k* for instance. Hence, Wilson (1942), Deresiewicz (1962), Dutta (1963) and Bhattacharya (1970) and many others studied the propagation of surface waves in such inhomogeneous media; however, they were interested in shear-horizontal waves (Love-type). The literature on Rayleigh-type surface waves in that type of media is quite scarce, probably because the difficulty exposed below is encountered quite early in the analysis.

In an anisotropic elastic body with continuously variable properties, the general equations of motion read(1.1)where * u* is the mechanical displacement, and

*C*

_{ijkl}and

*ρ*are the elastic stiffnesses and the mass density, respectively. Now consider the propagation of an inhomogeneous plane wave with speed

*v*and wavenumber

*k*in the

*x*

_{1}-direction and with attenuation in the

*x*

_{2}-direction,(1.2)in a half-space

*x*

_{2}≥0 made of an orthotropic1 material with an exponential depth profile,(1.3)Here, the

*x*

_{1},

*x*

_{2},

*x*

_{3}directions are aligned with the axes of symmetry,

*α*is a real number and the and

*ρ*° are constants. In addition,

*° is a constant vector and*

**U***q*a complex number so that the attenuation factor is

*k*Im(

*q*). Then the equations of motion (1.1) yield(1.4)At

*α*=0, the material is

*homogeneous*, and the associated determinantal equation—the propagation condition—is a real

*quadratic*in

*q*

^{2}which can be solved exactly (Sveklo 1948).

At *α*≠0, the propagation condition is seemingly a *quartic* in *q* with complex coefficients, whose analytical resolution might appear to be a daunting task and to preclude further progress towards the completion of a boundary value problem. (Note that it remains a quartic even when the material is isotropic.) Hence, Das *et al*. (1992) and Pal & Acharya (1998) stopped their analytical study of the problem at that very point. In fact, the transformation of the quartic to its canonical form reveals that it is a *quadratic* in *q*+i(*α*/*k*), with *real* coefficients. Biot (1965) seems to be the only one who has recognized this simplification (in the context of incremental static deformations). The present paper shows that the Stroh (1962) formulation of this problem, combined with a change of unknown functions, leads naturally to the biquadratic in question. Then the propagation condition can be solved exactly, and the general solution of form (1.2) to the equations of motion follows. In particular, the resolution of the dispersive Rayleigh wave boundary value problem poses no particular difficulty after all. Section 2 exposes this analysis and §3 applies it to two types of exponentially graded half-spaces: one which would be made of orthotropic shales if *α*→0 and another which would be made of silica (isotropic). There, it is seen for both the examples that the influence of the inhomogeneity is more marked upon the wave speed (rapidly decreasing with *α*/*k*) than upon the attenuation factors (slowing increasing with *α*/*k*). It is also found that the attenuation factors for the displacement amplitudes are *distinct* from those for the traction amplitudes, and that the amplitudes can decay in an *oscillating* manner for the isotropic silica. These two features are unusual and are clearly due to the inhomogeneity.

The overall aim of the paper is to show that simple, analytical and exact results can be obtained for seismic Rayleigh wave propagation in an anisotropic, inhomogeneous Earth. Of course, it is unlikely any ‘real’ inhomogeneity can be such that the stiffnesses and the mass density, all vary in the same manner as in equation (1.3), because it then leads to bulk wave speeds (proportional to the square root of stiffnesses divided by the density) which are constant with depth. The analysis of more realistic models must turn to numerical simulations such as those based on the finite difference technique or on the pseudospectral technique or on techniques with Fourier or other function expansions (e.g. Tessmer 1995). However, these methods encounter difficulties for the implementation of accurate boundary conditions and of strong heterogeneity. The spectral element method seems to alleviate those difficulties but, as stressed by Komatitsch *et al*. (2000), it must be validated against analytical solutions. Such a solution validation procedure is indeed a crucial necessity of numerical simulations in geophysics, where different software packages can give widely different predictions (Hatton 1997).

## 2. The dispersion equation

Consider the propagation of a Rayleigh wave, travelling with speed *v* and wavenumber *k* in the *x*_{1}-direction, in an inhomogeneous half-space *x*_{2}≥0 made of the orthotropic material presented in §1. The associated mechanical quantities are the displacement components *u*_{j} and the traction components *σ*_{j2} (*j*=1,2). They are now taken in the form(2.1)where the *U*_{j}, *t*_{j2} (*j*=1,2) are yet unknown functions of *x*_{2} alone, to be determined from the equations of motion and from the boundary conditions.

The equations of motion: *σ*_{ij,j}=*ρu*_{i,tt}, can be written as the second-order differential system (1.1), or as the following first-order differential system,(2.2)Here *N*_{1}, *N*_{2} and *K* are the usual constant matrices of Stroh (1962), given by(2.3)where and . With the new vector function * ξ*, defined as(2.4)the system (2.2) becomes(2.5)Hence, the apparently anodyne change of unknown functions (2.4) transforms the differential system with variable coefficients (2.2) into one with constant coefficients.

Now solve the differential system (2.5) with a solution in exponential evanescent form,(2.6)where * ζ* is a constant vector,

*p*is a scalar, and the inequality ensures that(2.7)because by (2.4) and (2.6)

_{1},

*(*

**u***x*

_{2}) behaves as exp

*k*(i

*p*+

*α*/

*k*)

*x*

_{2}and

*(*

**t***x*

_{2}) behaves as exp

*k*(i

*p*−

*α*/

*k*)

*x*

_{2}. Note in passing that, in sharp contrast to the homogeneous case, the displacement field and the traction field have

*different*attenuation factors: for

*it is*

**u***k*[Im(

*p*)−

*α*/

*k*]; for

*it is*

**t***k*[Im(

*p*)+

*α*/

*k*].

Then * ζ* and

*p*are solutions to the eigenvalue problem,

*N*=

**ζ***p*. The associated determinantal equation is the

**ζ***propagation condition*, here a

*biquadratic*(and not a quartic as equation (1.4) suggested),(2.8)where(2.9)Let

*p*

_{1}and

*p*

_{2}be the two roots of equation (2.8) satisfying inequality (2.6)

_{2}. This pair may be in one of the two forms:

*p*

_{1}=i

*b*

_{1},

*p*

_{2}=i

*b*

_{2}or

*p*

_{1}=−

*a*+i

*b*,

*p*

_{2}=−

*a*+i

*b*, where

*b*,

*b*

_{1}and

*b*

_{2}are positive. In both the cases,

*p*

_{1}

*p*

_{2}is a real negative number and

*p*

_{1}+

*p*

_{2}is a purely imaginary number with positive imaginary part. It follows in turn that(2.10)The associated eigenvectors

**ζ**^{1},

**ζ**^{2}are determined from

*N*

**ζ**^{j}=

*p*

_{j}

**ζ**^{j}, as(2.11)where the non-dimensional quantities

*e*

_{0},

*f*

_{1},

*f*

_{0}appearing in the displacement components are given by(2.12)and the quantities

*g*

_{1},

*g*

_{0},

*h*

_{0}(dimensions of a stiffness) appearing in the traction components are given by(2.13)Now construct the general solution to the equations of motion (2.5) as(2.14)where the constants

*γ*

_{1},

*γ*

_{2}are such that the surface

*x*

_{2}=0 is free of tractions,

*(0)=*

**t****0**or equivalently,

*(0)=[*

**ξ***(0),*

**U****0**]

^{t}. This condition leads to a homogeneous linear system of two equations for the two constants whose determinant must be zero. After factorization and use of equation (2.10), the dispersion equation follows as(2.15)This equation is fully explicit (

*X*is the sole unknown) because

*P*and

*S*are given in equation (2.9) and

*g*

_{1},

*g*

_{0},

*h*

_{0}are given in equation (2.13), and it is clearly

*dispersive*due to the multiple appearance of the dispersion parameter

*α*/

*k*. At

*α*=0 (homogeneous substrate), it simplifies to(2.16)the classic (non-dispersive) secular equation for Rayleigh waves in orthotropic solids.

## 3. Examples: exponentially graded shales and silica

As two examples of application, consider in turn that the half-space is made of a material with exponentially variable properties which is (i) with orthotropic symmetry, and (ii) isotropic.

In example (i), the starting point is a model proposed by Shoenberg & Helbig (1997), accounting for the vertical fine stratification and the vertical fractures found in many shales. In their numerical simulations, they used the following orthotropic elastic stiffness matrix,(3.1)Note that here the matrix is density-normalized so that its components have the dimensions of squared speeds, expressed in (km s^{−1})^{2} (Shoenberg 1994). Schoenberg and Helbig remark that ‘the rock mass behaves as if it contains systems of parallel fractures increasing the compliance in some directions’; integrating this information, *α* is assumed positive here. In addition, equation (3.1) is assumed to be the elastic stiffness matrix on the free surface *x*_{2}=0.

In example (ii), the half-space is assumed to be made of an exponentially graded material such that at the boundary, , (10^{10} N m^{−2}) and *ρ*°=2203 kg m^{−3} as in silica (Royer & Dieulesaint 2000). Here too, *α* is taken positive.

If the half-spaces were homogeneous, then the Rayleigh wave would travel with speed where *X* is given by equation (2.16), that is *v*°=1.412 km s^{−1} for shales and *v*°=3409 m s^{−1} for silica. For any given dispersion parameter *α*/*k*, the dispersion equation (2.15) in the inhomogeneous half-spaces gives a unique root *X*. In both the examples, it has then been checked that for *X*, the propagation condition (2.8) gives two roots such that the inequality (2.6)_{2} is always satisfied. Thus, the surface wave exists for arbitrary value of *α*/*k*, and it travels with speed . Although this state of affair is acceptable mathematically, it seems reasonable to limit the range of *α*/*k* to values where the wave amplitude decays faster than the inhomogeneity. Since the amplitudes of the tractions * t* decay as exp−

*k*[Im(

*p*)+

*α*/

*k*], they always decrease faster than exp−2

*αx*

_{2}by inequality (2.6)

_{2}. On the other hand, the amplitudes of the displacements

*decay as exp−*

**u***k*[Im(

*p*)−

*α*/

*k*]; thus, they decrease faster than the inhomogeneity as long as Im(

*p*)>3

*α*/

*k*. In example (i), it turns out that this latter inequality is verified for

*α*/

*k*<0.107, and in example (ii), for

*α*/

*k*<0.274.

Figure 1 shows the variation of the wave speed (decreasing) and of Im(*p*_{1}), Im(*p*_{2}) (increasing) in example (i) over the range 0≤*α*/*k*≤0.1. It has also been checked there that the attenuation factors for both the amplitudes of displacements (*k*[Im(*p*)−*α*/*k*]) and tractions (*k*[Im(*p*)−*α*/*k*]) also increase. In conclusion, the surface wave travels at a slower speed in the inhomogeneous shales than in the homogeneous shales, and it is less localized.

Figure 2 shows the variation of the wave speed (decreasing) and of Im(*p*_{1}), Im(*p*_{2}) (increasing) in example (ii) over the range 0≤*α*/*k*≤0.25. It has been checked again that the attenuation factors for both the amplitudes of displacements (*k*[Im(*p*)−*α*/*k*]) and the tractions (*k*[Im(*p*)−*α*/*k*]) also increase. Here again, the surface wave travels at a noticeably slower speed in the inhomogeneous case than in the homogeneous case, and it is slightly less localized. A most interesting phenomenon occurs at *α*/*k*≈0.211 where the nature of the roots changes from the form: *p*_{1}=i*b*_{1}, *p*_{2}=i*b*_{2}, to the form: *p*_{1}=−*a*+i*b*, *p*_{2}=*a*+i*b*, so that the amplitudes switch from decaying in a real exponential manner to decaying in an exponential oscillating manner. This latter situation *never* arises in a *homogeneous* isotropic half-space.

## Acknowledgments

I thank F. Watchorn for stimulating discussions.

## Footnotes

↵An anisotropic material belongs to the orthotropic symmetry class when it possesses three mutually orthogonal planes of mirror symmetry.

- Received May 21, 2006.
- Accepted September 5, 2006.

- © 2006 The Royal Society