Scattering of spherical elastic waves from a small-volume spherical inclusion

Summary

A theory is presented for the scattering of spherical elastic waves by a small, spherical inclusion. The domain of validity of the treatment is ka′« 1, where k is the wavenumber and a′ is the radius of the spherical inclusion. The treatment is an extension of that of Gritto, Korneev & Johnson (1995), who derived exact expressions for the scattered wavefield from a plane P wave incident on a spherical inclusion in the low-frequency limit. Since neither treatment assumes small perturbations, they are also valid for an arbitrary material property contrast between the spherical inclusion and the background medium. For spherical elastic wavefields propagating through an infinite, homogeneous medium and interacting with a single small, spherical inclusion, the scattered wavefield is obtained directly as a sum of spheroidal motions of degrees 0, 1 and 2 in scatterer-centred coordinates, with nine independent amplitude coefficients which depend on the contrast in material properties and are proportional to the volume of the inclusion. If the inclusion is embedded in a stratified and bounded medium, then these scattering coefficients provide a valid local description of the scattered wavefield in the vicinity of the inclusion. Exploiting this fact, source equivalents in the frequency domain (six moment tensor elements and three single forces) are derived from the scattering coefficients. To first order in material property contrasts, equivalence is established between these source equivalents and the coupled mode theory developed for surface waves and long-period body waves using perturbation theory.

1 Introduction

There are three main approaches to handling seismic wave propagation through a heterogeneous elastic medium. The first addresses the problem of elastic wave scattering by a spherical inclusion embedded in an infinite, homogeneous medium and involves the solution of a boundary value problem on the boundary of the spherical inclusion. The second approach is based on the expansion of the elastic wavefield into a sum of normal modes, and mode-coupling theory describes the interaction of the wavefields with 3-D, distributed heterogeneities throughout the Earth. The third approach uses the approximation of ray theory, in which the scale of heterogeneity is assumed to be much broader than the wavelength of the propagating seismic wave which samples it. Other approaches are possible (for example Mie scattering, WKBJ approximation or boundary integral methods), but these three adequately represent the main mathematical tools used in current applications. The first two approaches can be implemented to obtain exact estimates of the scattered wavefield generated by an incident wavefield on a suitable reference earth model.

Several studies have been devoted to comparing seismic wavefields calculated using the second and third approaches, for both surface waves and body waves, and the second and third approaches have been widely used in observational studies because they can naturally handle the effects of 1-D stratification in the elastic properties of the reference medium and the free surface of the Earth. Of the first two approaches, the first is more expedient because it does not require a decomposition of the scattered wavefield into normal-mode components, but additional work is required in order to handle the effects of stratification and boundedness of the medium. This problem can be solved by deriving source equivalents for the scattered wavefield, and the calculation of the scattered wavefield is then primarily a matter of choosing a convenient form for the incident wavefields and choosing an appropriate reference earth model.

Using the first approach, source equivalents for the scattering of incident plane P and S waves by a small, spherical inclusion embedded in an infinite, isotropic, homogeneous background medium were presented by Wu & Aki (1985) using perturbation theory. Their results require that |δμ/μ| « 1, |δκ/κ| « 1 and |δρ/ρ| « 1, where μ, κ and ρ are the shear modulus, bulk modulus and density, respectively, and δμ, δκ and δρ denote the contrast between the inclusion and the background medium. Korneev & Johnson (1993a,b) gave exact expressions for the scattered wavefields generated by plane P waves incident on a finite sphere embedded in an infinite, homogeneous medium and compared the finite-sphere results with equivalent results for the point-scattering treatment of Wu & Aki (1985) as well as with Mie scattering. They showed that a point-scattering description of the effects of interaction with a finite sphere is strictly valid only when the quantity ka′ (wavenumber times radius of spherical inclusion) is less than a critical number, generally referred to as the Rayleigh limit. Korneev & Johnson (1996) present a detailed review of their earlier work and include results for the scattered wavefields resulting from an incident S wave.

Gritto, Korneev & Johnson (1995) extended the treatment of Wu & Aki (1985) to allow for low-frequency scattering from a spherical inclusion with arbitrary contrast in material properties between the inclusion and the background medium. Their treatment took the lowest-order terms in the analytic solution presented by Korneev & Johnson (1993a), and it is valid provided that ka′« 1. Their results are for incident P waves but can be adapted to the case of incident S waves. Derivation of source equivalents for this case could then be done straightforwardly and would yield a description of the wavefield scattered from incident plane waves in terms of equivalent sources (a moment tensor and single forces); this would extend the source equivalents presented by Wu & Aki (1985) to first order in material property contrasts.

The applications discussed above have almost exclusively been devoted to the case of incident plane P or S waves. It is worthwhile, however, to address the problem of the scattering of incident spherical waves by a spherical inclusion in the domain of validity of Rayleigh scattering, ka′« 1, the same domain as investigated by Gritto . (1995). Taking an r-θ-φ spherical coordinate system, incident spherical waves have the form

where u_o(r, ω) is the displacement spectrum at frequency ω observed at r = (r, θ, φ), and are potential functions defined on the unit sphere satisfying for a given wavenumber k. The radial functions and represent P-SV- (or Rayleigh-) wave propagation, and the function represents SH- (or Love-) wave propagation. The strongest advantage of working with spherical waves is that all seismic wavefields of interest on both a global and a regional scale can be expressed as a superposition of waves of the form of eq. (1). The radial functions , and are almost always chosen to satisfy stress-free boundary conditions at the Earth's surface, so that waves of this form are directly comparable with observations made at the Earth's surface. We may point out that the scattering theory necessary to handle spherical waves of various characters has been previously worked out. This has been done in the context of free oscillations (Woodhouse 1980; Mochizuki 1986; Park 1987; Romanowicz 1987; Dahlen 1987; Lognonne & Romanowicz 1990), surface waves (Snieder & Romanowicz 1988) and body waves (Li & Tanimoto 1993; Li & Romanowicz 1995; Marquering, Snieder & Nolet 1996; Zhao & Dahlen 1996). Mode-mode coupling and perturbation theory are built into all of these treatments; the treatment of Zhao & Dahlen (1996) goes further, in asymptotically evaluating the summations over mode branches and spherical harmonic degrees for rays characterized by a particular ray parameter. These methods provide a powerful framework for seismic wavefield computation when the Earth can be adequately represented in terms of small perturbations from a spherically symmetric reference model, particularly when signals of very long duration and relatively low frequency are required.

The purpose of this paper is to derive source equivalents for the scattering of an incident spherical wave by a spherical inclusion in the domain of Rayleigh scattering. Since wave-number k= frequency/velocity, the condition ka′« 1 is satisfied in either the low-frequency limit or the small-inclusion-radius (or volume) limit We shall focus particularly on the small-inclusion-volume limit, which not only satisfies the condition of Rayleigh scattering but also justifies the description of the scattered wavefield in terms of equivalent sources, even when the inclusion is embedded in a stratified and bounded background medium. The first step of this new approach consists of formulating and solving a boundary value problem on the surface of a spherical inclusion subjected to an incident spherical wavefield. Source equivalents are then derived by identifying the outgoing scattered wavefields generated at the boundary of the spherical inclusion with the well-known wavefields generated by particular seismic sources. As tests of the derived expressions for the source equivalents, I compare the theoretical expressions for the scattered wavefields obtained by other methods with the corresponding wavefields obtained from the source equivalents. These tests are, further, intended to illuminate the physical connection between the present approach and surface-wave-coupling theory.

The philosophy adopted here is similar to the T-matrix method of Bostock (1991) and Bostock & Kennett (1992), who investigated surface-wave scattering under similar conditions and whose results also do not depend on the limits of perturbation theory. The present paper considers the problem of scattering from a small embedded sphere in the presence of large material property contrasts, with the aim of providing a framework for more sophisticated treatments designed to handle arbitrarily large, strongly heterogeneous volumes. It allows an exact description of the scattered wavefield from a single small, spherical inclusion embedded in a stratified background medium in terms of equivalent sources. Even in the presence of small structural perturbations, such a description is of high general utility because it allows the calculation of the scattered wavefield using the same tools that have been developed to calculate seismic wavefields in layered earth models from ordinary seismic sources. Such a treatment would facilitate waveform modelling of long-period body waves, where it may be necessary to incorporate both amplitude and phase information and where distortions of the wavefield may be produced by wave propagation through heterogeneous material in the crust and upper mantle. Direct observational evidence for scattered long-period body waves in regional networks (Neele & Snieder 1991) and local networks (e.g. Frankel & Vidale 1992) strongly encourages such an approach. The approach in this paper may contribute to analytic methods of calculating seismic wavefields in the presence of strong structural perturbations, because perturbation theory is ill-suited to deal with such phenomena as fault-zone guided waves (Li & Vidale 1996) or amplification and resonance within sedimentary basins. In the latter case, a fair measure of success has been obtained by analytical techniques which employ exact solutions for wave propagation within both a sedimentary basin and the surrounding crustal material and which match solutions at a semi-spherical boundary (Trifunac 1971) or an irregular boundary (Kawase 1988). These approaches are obviously designed for highly specialized geometries, but they are good examples of what types of problems can be attacked by semi-analytical methods without using perturbation theory. In an investigation also related to the present approach, Cormier (1995) studied D″ structure by representing it in terms of numerous finite spheres with large contrast to the background medium and subjected to incident wavefields calculated using dynamic ray tracing. The D″ region may contain packets of heterogeneous material with scale lengths of 100 m and very large contrasts (Knittle & Jeanloz 1991), although longer-wavelength averages yield contrasts of only 0.5–2 per cent (Bataille, Wu & Flatte 1990; Cormier 1998). The possibility that the D″ layer may be composed of an enormous number of strong, small-scale heterogeneities indicates that calculation of the complex wavefield effects in that region may be better understood by applying a treatment based on scattering by small spheres with large contrasts.

2 Preliminaries

2.1 Geometry and formulation of the problem

Referring to Fig. 1, we assume a spherical earth geometry (r, θ, φ) centred on O, and denote the laterally homogeneous earth model by A. A homogeneous sphere embedded in A and with generally different material properties from A is denoted by B. It is taken to have a centre O′, located at r₀ with respect to O. The radii of the earth and of the spherical inclusion are a and a′, respectively. Spherical shells at radii and are defined as the top and bottom radii, respectively, of region B with respect to O. It is assumed that the material outside B is homogeneous between the spherical shells of radii r₀ and r₁.

We first consider the problem of determining the displacement field excited by a seismic source in the laterally homogeneous model A. We take a seismic source at r_s, with coordinates (r=r_sθ= 0). It is assumed that the bulk modulus κ(r) shear modulus μ(r) and density ρ(r) within A are laterally homogeneous outside of B, depending only on radius. For simplicity, we assume that the corresponding bulk and shear moduli and density within B are constant, and denote them by κ′, μ′ and ρ′.

The displacement field for a point dislocation is obtained as a summation of spherical-harmonic wavefields. Let represent the fully normalized spherical harmonic of total degree l and azimuthal order number m (Edmonds 1960). Explicitly, for positive m we define

where is the associated Legendre polynomial and where · denotes complex conjugation. We define the position vector r = (r, θ, φ). In the frequency domain the displacement field may be represented in terms of a basis set of spheroidal and toroidal wavefields, with

for spheroidal motion and

for toroidal motion. In eqs (3) and (4), ▿₁ is the surface gradient operator

The corresponding normal and shear tractions on the spherical shell of radius r may be expressed in the form

for spheroidal motion and

for toroidal motion, where T denotes the stress tensor in A, which is assumed to be an isotropic medium. The total displacement and traction fields are given by

In the time domain, the equation of motion in the coordinate system centred on O has the form

where M is a moment tensor and F₀ a single force vector acting at r = r_s. Defining the Fourier transform of a function f(t) as

the transform of eq. (9) is

This is subject to the boundary conditions (where a is the earth's radius) and the condition that the solution must be regular at the origin r= 0. Anticipating subsequent manipulations with l and m, we let k denote the pair (l, m). We define the displacement-stress vector via

for spheroidal motion and

for toroidal motion. Eq. (10) possesses unique solutions of the form (8a,b), composed of a sum of spheroidal- and toroidal-motion wavefields. With the definitions (11) and (12) it is convenient to rewrite the total displacement and traction fields (8a,b) in the more compact form

where the various differential operators appearing above are row vectors of the form

The application of operators and to a displacement-stress vector yields the corresponding displacement vectors; the application of operators and yields the corresponding tractions.

Henceforth we shall assume a fixed ω and an implicit dependence of all displacement-stress vectors upon this frequency. For the laterally homogeneous earth model A, eq. (10) may be rewritten as a system of six first-order ordinary differential equations, with a source term which is discontinuous at r=r_s:

where A is a 4 × 4 matrix for spheroidal motion and a 2 × 2 matrix for toroidal motion. These equations have been given by Takeuchi & Saito (1972). The form of the source term f and a stable method of numerically integrating eqs (16) have been given by Friederich & Dalkolmo (1995). If eq. (16) is solved in terms of a stack of layers, then within each layer or may generally be written as the sum of four or two independent solutions, respectively. In particular, in the region r₀ < r < r₁ we may write

where the solutions and correspond to basis solutions which are regular at the origin r= 0, and the solutions and correspond to basis solutions which are singular at the origin r= 0 and represent outgoing waves. Explicit forms for these basis solutions are required for future manipulations. They are given in Appendix A and were derived from eqs (6.2)–(6.4) of Ben-Menahem & Singh (1981) by noting the following correspondence:

where

and the vector spherical harmonics N^±, L^± and M^± are defined by eq. (2.85) of Ben-Menahem & Singh (1981).

Within the spherical inclusion (medium B), let u′ and T′ represent displacement and stress. We also denote the corresponding displacement-stress vector by y′ and use representations of displacement, shear and normal stresses in the forms of eqs (3) and (4) or (6) and (7), using O′-centred coordinates r′= (r′, θ′, φ′) instead of O-centred coordinates r = (r, θ, φ). The equation of motion in the spherical inclusion may then be written in the form

In terms of the displacement-stress vectors, eq. (20) may be written as the system of differential equations

In eq. (21) it is understood that A′ is A defined with κ′, μ′ and ρ′ replacing κ(r), μ(r) and ρ(r). Because homogeneous elastic parameters are assumed in the spherical inclusion, the solutions in B must be of the form

With the above remarks and definitions, the strategy for solution of eq. (10) in the presence of the spherical inclusion is as follows.

• (1)

  Start with an appropriate solution of eq. (10) in the laterally homogeneous earth model A centred on O (in most applications this solution will be a wavefield excited by a seismic source and will be referred to as the incident wavefield), and represent this in terms of a superposition of solutions of eq. (10) centred on O′.

• 2

  Match this solution with the solution of eq. (20), valid in B, on the boundary r′=a′ by introducing an outgoing scattered wave into A from the boundary of B.

• 3

  In the small-volume limit, only a finite number of outgoing waves are generated and this solution may be described exactly in terms of equivalent sources.

In Section 2.2 we describe steps (1) and (2) in greater detail, and in Sections 3.1–3.3 we discuss step (3) in terms of the structure of the scattered wavefields and the equivalent sources.

2 Wavefield matching at boundary of spherical inclusion

To transform the incident wavefield from the O-centred into the O-centred spherical coordinate system requires manipulation of solutions of eq. (10) under translation and rotation of the origin of the spherical coordinate system. This transformation is derived in Appendix B. As given in eqs (B6) and (B7), we may transform the incident wavefield, which is a solution of (10) (i.e. eq. 17), from O-centred to O′-centred coordinates by means of a superposition of the form

where χ denotes the pair (ν, μ), and the differential operators and have the same structure as and in eq. (15) but are defined in the (r′, θ′, φ′) coordinate system:

where

Exterior to the inclusion, the displacement-stress vectors which appear on the right-hand side of eq. (23), and , must have the form

Exact expressions for the coefficients could be obtained using the transformations contained in Appendix J3 of Ben-Menahem & Singh (1981) and would involve a countably infinite number of coefficients. In the following section I will show how only a finite subset of these coefficients is required in the small-volume limit.

The parametrization (26), applied to each displacement-stress vector component which appears on the right-hand side of eq. (23), allows us to express the solution of eq. (10) with respect to O in terms of solutions of eq. (10) with respect to O′. In general it is necessary to match the solutions of eq. (20) with solutions of eq. (10) on the spherical boundary r′=a′. This requires us to introduce an outgoing scattered wave from the spherical inclusion. As explained below, we may assume that the outgoing scattered wave propagates through a homogeneous medium outside B, with material properties equal to those of medium A in the volume r₀ < r < r₁. In this boundary value problem we require that the sum of the incident wavefield, currently specified by the coefficients and in eq. (26), and an outgoing scattered wavefield equals a solution of eq. (20) on r′=a′. We arbitrarily replace the dummy variable χ on the right-hand sides of eqs (23) and (26) by κ. Since the outgoing waves outside B must be represented by the functions or , we have

Because all of the displacement-stress vectors in eq. (27) are defined with respect to the scatterer centre O′, it decouples into independent 4 × 4 or 2 × 2 linear systems, one for each component k.

The terms proportional to and on the left-hand side of eq. (27) represent the outgoing scattered spheroidal or toroidal wavefield evaluated at r′=a′ generated by interaction of the incident wavefield component with the inclusion; the terms proportional to and represent the spheroidal or toroidal wavefield within the inclusion evaluated at r′=a′. If the inclusion were embedded in an infinite, homogeneous medium, then the left-hand side would represent exactly the total wavefield in the background medium (i.e. for r′ > a′). Since medium A will generally have stratification and possess internal boundaries and a free surface, however, the scattered wavefield will interact with these boundaries and produce a new wavefield which should, for completeness, be treated as an additional incident wavefield. In the small-volume limit, however, such internal boundaries (or significant variations due to smoothly varying stratification) will be several inclusion radii away from the boundary of the inclusion. Since we are in the Rayleigh scattering regime, where the wavelength is much greater than the scatterer size, it suffices to consider the near-field dependence of the scattered wavefield on the distance r′ from the centre of the scatterer. It is shown in Section 3.1 that the scattered wavefield decays rapidly away from the inclusion boundary, with an r′⁻¹ or r′⁻² dependence [this result also follows immediately from eqs (9) and (15) of Gritto . 1995]. Thus the wavefield away from the inclusion varies as (r′/a′)⁻¹ or (r′/a′)⁻² relative to the scattered wavefield at the inclusion boundary. Provided that the distance between the scatterer and an internal boundary or the free surface is sufficiently large, (r′/a′)⁻¹« 1. and the additional wavefields generated by reflection from internal boundaries or the free surface will have negligible interaction with the inclusion boundary; in most applications this condition should be easily met. It is instructive to note, however, that this condition would not generally be satisfied if several closely spaced spherical inclusions were employed to represent a large heterogeneous volume. In that case, the distance r′ between a spherical inclusion and its closest neighbours would be comparable with its radius a′, and to obtain accurate results it would be necessary to consider the multiple interactions of these spheres with one another, in the sense investigated by Bostock & Kennett (1992) for interactions among closely spaced cylinders.

3 Scattering From A Small-Volume Inclusion

We assume that an incident spherical wavefield prescribed by and in eq. (1), in spherical coordinates (r, θ, φ), interacts with a scatterer centred on r₀. We further assume that the expansion coefficients or of this incident wavefield in scatterer-centred coordinates (i.e. eqs 23 and 26) are available. It then remains to determine the scattered wavefield by solving the boundary value problem on the surface of the spherical inclusion (eq. 27) for each component of this expansion. Corresponding results for more complicated incident wave-fields (i.e. those of the form of eq. 8) could be obtained by linear superposition.

For a small inclusion it is appropriate to consider the outgoing waves generated at the boundary r′=a′ in the limit a′→ 0, and we require solutions of eq. (27) in this limit. In the small-a′ case we take advantage of the fact that the displacement–stress vectors y in O′-centred coordinates can be written in terms of asymptotic expressions for the various Bessel functions. We first seek to identify which wavefields in scatterer-centred coordinates are excited by interaction of the incident wavefield with a point scatterer. We then obtain exact expressions for the required scattering coefficients and derive source equivalents for the scattering interaction in terms of the incident spherical wavefield.

3.1 Structure of wavefields generated by scattering from a small inclusion

For spheroidal motion, in order to retain linearly independent solutions of the equations of motion for small r′ we must take linear combinations of the displacement-stress vectors within both the inclusion and the background medium. We rewrite eq. (27) in the form

where the vectors are defined by eq. (A4). The incident wavefield is parametrized by the coefficients and via eq. (26). As r′ approaches zero, we may take advantage of asymptotic expressions for and in powers of r′ via eqs (A6). These expressions were derived under the assumptions

where k^α=ω/α and and are the wavenumbers defined in the background medium, and and are the wavenumbers defined within the inclusion. These conditions constitute Rayleigh scattering. They are less restrictive than the Born approximation, in which solutions for the scattered wavefield are linearized with respect to the contrasts in material properties δλ=λ′–λ, δμ=μ′–μ and δρ=ρ′–ρ between the inclusion and the local background medium. This linearization requires that

A necessary, but not sufficient, condition for the validity of the Born approximation is that the scattered wavefield evaluated just exterior to the inclusion is small compared with the incident wavefield. As pointed out by Gritto . (1995), it is a valid approximation when the conditions for both Rayleigh scattering (eq. 29a) and small contrasts in material properties (eq. 29b) are satisfied. The present treatment, like that of Gritto . (1995), requires only the Rayleigh approximation (eq. 29a), and not the linearization requirement (eq. 29b). Under these more general conditions the scattered wavefield just outside the inclusion boundary may be of the same order of magnitude as the incident wavefield.

We may now state the simplification resulting from adopting the condition of Rayleigh scattering. If the outgoing wavefield is expanded in a Taylor series in powers of a′, then the lowest-order term is ∼a′³, only spheroidal-motion wavefields of degrees 0, 1 and 2 in inclusion-centred coordinates contribute to these lowest-order terms, and the higher-order terms are negligible. In order to demonstrate the structure of the scattered wavefields, we rewrite eq. (28) in the form

Eq. (30a) is directly applicable for l≥2 but must be further manipulated for l= 0 and l= 1.

In the case k= (0,0) (l= 0, m= 0), the horizontal displacement and shear stress vanish, the vertical displacement and normal stress associated with are trivial, and the vertical-displacement and normal-stress components of and are linearly dependent. Arbitrarily eliminating , (30a) thus reduces to and

The subscripts 1 and 2 here refer to the elements occupying the first or second row respectively of the given column vector.

In the case k= (1,0) or k= (1, ± 1) (l= 1), the leading terms in the series expansion of the left-hand and right-hand sides of eq. (30a) correspond to a rigid translation of the spherical inclusion and the background medium, as may be seen from the asymptotic forms of the functions and in eq. (A6). This requires that where is of order a′² relative to . Removal of this translation component from eq. (30a) yields, to next leading order in a′,

where and are given by eqs (A7); b and b′ on the right-hand side of eq. (30c) are evaluated with the material properties of the background medium and spherical inclusion, respectively.

The matrix elements required on the left-hand sides of eqs (30a–c), consisting of elements of the column vectors (background medium) or (inclusion) evaluated at the radius r′=a′, are readily supplied from eq. (A6). The matrix elements required on the right-hand sides of eqs (30a–c), consisting of elements of the column vectors and (background medium) or (inclusion) evaluated at r′=a,′ are specified by eqs (26), (A6) and (A7).

It can be shown that the determinants of the matrices on the left-hand sides of eqs (30a–c) are always greater than zero (up to multiplication by a power of a′). Using the asymptotic expressions (A6) for the column vectors which appear on the left-hand side of eq. (30a), we find that for l≥ 2 there are solutions and with proportionalities (relative to and and . Letting V denote the volume of the spherical inclusion (V= (4π 3)a′³), and thus and , it is immediately clear that spheroidal wavefields of order 2, but not of higher order, may be generated by interaction with a small inclusion, and that the corresponding scattered wavefield is proportional to . In the case l= 0 it follows from eqs (30b) and (A6) that , so that scattered waves proportional to are generated. Similarly, in the case l= 1, eqs (30c) and (A6) show that and , so that scattered waves proportional to are generated. We may add that the structure of the scattered wavefields described above immediately reveals its near-field behaviour. From the asymptotic properties of or in eq. (A6), it follows that the near-field displacement field associated with the scattered wavefield varies with distance r′ from the scatterer as r′⁻¹ for the l= 1 component and as r′⁻² for the l= 0 and l= 2 components.

For toroidal motion, the second of eqs (27) may be written as

The required matrix elements in this case, (background medium) or (inclusion), evaluated at r′=a′, are provided by eq. (A8).

From the asymptotic forms of in eq. (A8), it follows that as a′→ 0, . A non-trivial scattered wave-field is then only possible for l= 1 (since toroidal wavefields of degree 0 do not exist) and would correspond to rigid rotation of the spherical inclusion. In the case l= 1, however, substitution of the required matrix elements from eq. (A8) into eq. (31) yields the result . Thus there is no scattered toroidal-motion wavefield generated by interaction of an elastic wavefield with a small-volume inclusion in scatterer-centred coordinates.

3.2 Explicit solutions for the scattered wavefield

The solutions of eq. (30) for the leading scattering coefficients are obtained by analytic inversion of the 4 × 4 linear systems (30a) and (30c) or the 2 × 2 linear system (30b). The results are as follows:

l= 0:

l= 1:

l= 2:

Here δλ=λ′–λ, δμ=μ′–μ and δρ=ρ′–ρ are the contrasts in material properties between the inclusion and the local background medium. Eqs. (32) allow us to characterize the scattered wavefield in terms of nine scattering coefficients which depend on the incident wavefield and , their spatial derivatives evaluated at the centre of the inclusion, the volume of the inclusion and the contrast in material properties. The dependence on the incident wavefield is implicitly contained in the expansion coefficients and of the incident wavefield in scatterer-centred coordinates, as shown in Appendix B. With these scattering coefficients we may express the scattered displacement or traction vector field in the vicinity of the inclusion as the finite sum

Since all of the matrix elements on the left-hand and right-hand sides of eq. (30) are subject to formal errors of order (ka′)², where k is the appropriate wavenumber for either the background medium or the inclusion, the scattering coefficients listed in eq. (32) and all of the quantities derived from them in the remainder of this paper implicitly contain formal errors of relative order . That this quantity should be close to unity is a restatement of the condition (29b) for the validity of Rayleigh scattering. This condition can always be satisfied for sufficiently small a′, provided that all of the elastic moduli are finite. One class of problems which fails to satisfy this condition is that of scattering from a fluid-filled inclusion. In that case, μ′= 0, and the product k′_βa′ is infinite for all a′.

Since the inclusion is not embedded in an infinite, homogeneous medium, the spatial domain over which eq. (33) is valid depends on several complicated factors such as the gradients in elastic stratification near the inclusion, the locations of discontinuities in elastic stratification within the earth and the presence of the free surface. Provided that the inclusion is not located directly on a discontinuity, however, the domain of validity will always satisfy r′/a′» 1 for sufficiently small a′. That is, eq. (33) is a valid description of the scattered wavefield in the vicinity of the inclusion. A similar argument is made in Section 4.5.4 of Ben-Menahem & Singh (1981) to justify the derivation of source jumps, and here it justifies comparisons of wavefields of the form (33) with the wavefields resulting from particular seismic sources; we shall now derive source equivalents using such comparisons.

3 Source equivalents

We next obtain a description of the outgoing scattered waves in the region r′ > a′ in terms of equivalent sources. The orientation of the spherical coordinate system (r′, θ′, φ′) centred on r₀ is defined by three Cartesian axes :

and these axes have been defined such that and that points along the minor arc connecting and in the direction away from (Fig. 2). From the structure of the nine different scattered wavefields determined in Section 3.2, combined with Table 4.4 of Ben-Menahem & Singh (1981), I find that the equivalent sources are moment tensor elements or single forces of the following four classes:

• 1

  l= 0 centre of compression (M₃₃=M₁₁=M₂₂=M);

• 2

  l= 1: concentrated force (F₀ in , or direction);

• 3

  l= 2, m= 0: compensated linear vector dipole ;

• 4

  l= 2, m=±1: double couple (M₁₃, M₂₃ non-zero);

  l= 2, m=±2: double couple (M₁₂ or M₁₁−M₂₂ non-zero).

I have derived the non-zero concentrated forces and moment tensor elements listed above from eqs (A2), (A3) and (A4), and from Table 4.4 of Ben-Menahem & Singh (1981), using the correspondences in eq. (18) to rewrite numerous expressions involving vector spherical harmonics contained in Ben-Menahem & Singh (1981) as equivalent expressions using the displacement-stress vectors defined in this paper. The displacement or traction vector corresponding to the k= (l, m) component of the scattered wavefield is given by a single term of eq. (33). For each k, this may be equated with a corresponding expression for the displacement or traction vector generated by a particular source component extracted from Table 4.4 of Ben-Menahem & Singh (1981). For example, in class (1) above, for k= (l= 0, m= 0), we find

leading to the expression for M given in eq. (36). In class (2), for k= (l= 1, m=±1), defining , we have

where the Legendre function is defined by . Equating real and imaginary parts in the above expression and noting that k_α/k_β=β/α yields eq. (38) below. The other source terms may be derived through similar manipulations, and I summarize the results in terms of the scattering coefficients or :

Eqs (36)–(41) are the main result of this paper and allow us to characterize the scattered wavefield in terms of equivalent sources M(ω) and F₀(ω) in the frequency domain. Since the scattering coefficients which appear on the right-hand sides of these equations depend on the incident wavefields their spatial derivatives evaluated at the centre of the inclusion, the volume of the inclusion and the contrast in material properties, the equivalent sources also depend on these factors. In typical applications these source equivalents could be used to calculate the scattered wavefield within the entire earth by simply solving eq. (10) with them, using well-established methods (e.g. Friederich & Dalkolmo 1995 and references therein).

In many applications it may be useful to have concrete expressions for source equivalents in the frequency domain under the assumption of small material property contrasts. These source equivalents may be derived by substituting the expressions (B16) and (B17) for the nine incident-wavefield expansion coefficients into eqs (32), linearizing with respect to δ_κ, δ_μ and δ_ρ, and then substituting the resulting scattered-wavefield coefficients into eqs (36)–41) to obtain the nine equivalent source terms. For general reference, we summarize in Table 1 the six moment tensor and three force elements for each incident wavefield type in the scatterer-centred coordinates depicted in Fig. 2.

4 Comparison With Previous Results

In this section we compare our expressions for the scattered wavefields generated by particular incident wavefields with those determined by other methods. These comparisons serve the main purpose of validating the scattering coefficients (eq. 32) (from which the source equivalents in general form, eqs 36–41. were derived), as well as the source equivalents in linearized form (Table 1). They also show that a broad class of problems is amenable to analysis under the fairly general conditions of Rayleigh scattering.

4.1 Plane P-wave incidence

By applying the operator through eq. (33) the scattered-wavefield displacement vector can be written in the form

The subscripts 1 and 3 here refer to the elements occupying the first or third row respectively of the given column vector. From eqs (A2)–(A4) it is seen that vectors represent waves travelling with the P-wave velocity, and the vectors are combinations of waves travelling with S-wave and P-wave velocities. The normalized coefficients and depend only on the material properties of the background medium and the inclusion and not on the radius of the inclusion.

With eq. (42), we can compare our results with those of Gritto . (1995) by constructing an incident P wave and directly comparing the scattering coefficients obtained in that paper with those obtained here. We construct an incident plane P wave propagating in the direction (through the pole 0′= 0 of the spherical coordinate system) through an infinite, homogeneous background medium. From eq. (2.141) of Ben-Menahem & Singh (1981), the potential exp(–ik·r′) with can be expanded as (assuming the spherical-harmonic conventions of this paper)

The incident plane P-wave displacement is

From eqs (44a) and (A3) it follows that the incident wavefield is a superposition of -type displacement–stress vectors. We prescribe a spherical inclusion at the origin r= 0 with volume and material property contrasts with the background medium such that the conditions for Rayleigh scattering are satisfied (i.e. eq. 29a). For point scattering we need only consider the following non-zero incident-wavefield expansion coefficients, which can be derived by rewriting eq. (44a) in the form

and retaining only the l≤ 2 terms:

Substitution of eqs (45) into eqs (32) yields the required scattered-wavefield coefficients. It is found that these coefficients are related to the coefficients a₁, for the P-wave component of the scattered wavefield derived by Gritto . (1995) by

The P-wave components of the spatial factors and which multiply these scattering coefficients are, from eq. (A4), proportional to the displacement components of the displacement-stress vectors and , respectively. Comparison of eq. (3) of Gritto . (1995) with eq. (42) of the present paper shows that the corresponding spatial factors involving the a₁, are in the ratios given by eq. (46), which proves that the scattering coefficients given in eq. (32) agree exactly with those derived by Gritto . (1995) for the scattered P wave. It also follows from the definition of the z vectors in eq. (A4) that the coefficients b₁, of the S-wave components of the scattered wavefield in eq. (3) of Gritto . (1995) should obey the relation

This relation is consistent with the a₁, and b₁, values given in eq. (4) of Gritto . (1995), showing that the S-wave components of the scattered wavefield in the two treatments are also in agreement.

4.2 Surface-wave incidence

The steps of the comparison in this section are, first, to write a suitable general expression for the displacement spectrum of a particular surface-wave mode branch resulting from excitation by a point source with a moment tensor spectrum M(ω) and single-force spectrum F₀(ω); second, to define a spheroidal- or toroidal-mode wavefield Ψ₀(r, θ, θ) incident upon a scatterer located at r₀; third, to make use of the source equivalents listed in Table 1 to characterize the scattered wavefield; fourth, to utilize our general expression for the displacement spectrum to obtain concrete expressions for the scattered wavefield in terms of Ψ₀; and, fifth, to compare these expressions with corresponding expressions derived by other methods.

Let us define an r-δ-θ′ spherical coordinate system to have δ= 0 passing through the scatterer epicentre, r defined as the distance from the centre of the Earth and θ′ defined as the azimuth (Fig. 2). Here θ′ plays the role of a scattering angle measured counter-clockwise from the forward-scattering direction. In this geographic coordinate system, the scattered wavefield will have both spheroidal- and toroidal-mode components. (This should not be confused with the result that no scattered toroidal-motion wavefields are generated in scatterer-centred coordinates in the vicinity of the spherical inclusion, as discussed in Section 3.1.) A complete proof of the equivalence between the derived source equivalents and normal-mode coupling theory would require the consideration of both incident spheroidal- and toroidal-mode wavefields, as well as the spheroidal- and toroidal-mode scattered wavefields which they generate. Since the cases of spheroidal → spheroidal and toroidal → spheroidal cross-branch coupling suffice for a test of the validity of the source equivalents, I will restrict attention to the scattered spheroidal-mode component generated by interaction of an incident spheroidal- or toroidal-mode wave-field with a spherical inclusion. Anticipating application to the scattered surface wavefield, the response of a laterally homogeneous earth to a point source with moment tensor spectrum M and single-force spectrum F₀ for the first-arriving surface waves on mode branch S′ will be denoted by the potential δΨ_s′n′, where n denotes an incident spheroidal (S) or toroidal (T) mode branch. At fixed ω, the corresponding first-order scattered spheroidal-mode wavefield can be written in the form

In eq. (48), S′ denotes the mode branch number, and and are evaluated for that mode l_s′ on the given branch corresponding to frequency ω. For spheroidal-mode wavefields (defined by and , with l_S corresponding to the mode on a particular spheroidal mode branch at frequency ω) and for toroidal-mode wavefields (defined by , with l_T corresponding to the mode on a particular toroidal mode branch at frequency ω), we shall assume the normalization of the radial eigenfunctions adopted by Dahlen (1980):

We define ν=l_S′+ 1/2. An expression for δΨ_S′,n′ is readily available from eqs (19) and (35) of Snieder & Nolet (1987), which give the displacement field in the form of eq. (48). The results of Snieder & Nolet (1987) depend on the far-field approximation to the Legendre functions , but they can be generalized to use a uniformly valid asymptotic approximation. In Appendix C I derive a uniformly valid asymptotic expression which, further, assumes the e^iωt sign convention of the Fourier transform adopted in this paper. From eqs (C10) and (C11) the S′th spheroidal wavefield potential is evaluated as

where C_S′=ω v is the phase velocity, U_S′ is the group velocity, the Σ functions are given by eq. (C4) and depend implicitly on the incident-mode branch index n, and G_m is defined in eq. (C12).

Comparison with surface-wave mode-coupling theory can be made by substituting the source equivalents listed in Table 1 into eq. (C4) to determine the Σ functions and then evaluating eq. (50). The resulting theoretical expressions for δΨ_S′n can then be compared with corresponding expressions obtained in other studies.

For spheroidal → spheroidal mode scattering (n=S), we first consider the isotropic (azimuthally symmetric) scattered wave-field generated by an incident spheroidal-mode wavefield, with k=l_S+ 1 2 corresponding to the wavenumber at frequency ω on the incident spheroidal-mode branch. Let and be the corresponding radial displacement functions and the potential which describes this branch (e.g. eq. 1), and let the centre of the scatterer be located at r₀ with material property contrasts δ_μ, δ_κ and δ_ρ with the local background medium. In the case of isotropic scattering we require the contributions from the three equivalent moment tensor components M₁₁, M₂₂ and M₃₃ and the force component Substitution of these components from Table 1 into the first of eqs (C4) yields

where

If the volume V is written as V= (r²dr)(dA). where dA denotes an area element on the unit sphere, then substitution of eq. (51) into eq. (50) yields

where

Eq. (52) is identical to eqs (43), (82) and (87) of Pollitz (1994) for isotropic surface-wave scattering, using the kernel values for given in eqs (A1)(A3) of that paper [terms and in place of δκFF′ should be included in the right-hand side of eq. (A1) in that paper].

Similar agreement may be demonstrated for those scattering terms which depend on m= 1 (sin φ′, cos φ′) or m= 2 (sin 2φ′, cos 2φ′). For the m= 1 scattered wavefield we require the contributions from the two equivalent moment tensor components M₁₃ and M₂₃ and force component . Substitution of these components from Table 1 into the second of eqs (C4) yields

Substitution of eq. (54) into eq. (50) yields

where

For the m= 2 scattered wavefield we require the contributions from the three equivalent moment tensor components M₁₁, M₂₂ and M₁₂. Substitution of these components from Table 1 into the third of eqs (C4) yields

The same forms for δΨ_S′S in eqs (55) and (58) are obtained from eqs (43), (82), (88) and(89) of Pollitz (1994).

For toroidal → spheroidal mode scattering (n=T), we may consider the scattered spheroidal-mode wavefield generated by an incident toroidal-mode wavefield, with k=l_T+ 1/2 corresponding to the wavenumber at frequency ω on the incident toroidal mode branch. Let be the corresponding radial displacement function and the potential which describes this branch. In this case substitution of the source equivalents from Table 1 into eq. (C4) shows that Σ₁= 0; that is, the isotropic component of the scattered wavefield is identically zero. The m= 1 (sin φ′, cos φ′) components of the scattered wavefield may be derived by evaluating the azimuthal factor Σ₂(φ′) with M₁₃, M₂₃ and F₀ taken from Table 1, with the result

The m= 2 (sin 2φ′, cos 2φ′) components of the scattered wavefield may be derived by evaluating the azimuthal factor Σ₃(φ′) with M₁₁, M₂₂ and M₁₂ taken from Table 1, with the result

As before, eqs (61) and (64) are found to agree with the corresponding results in Pollitz (1994) (eqs 43 and 90–94 of that paper).

It may also be shown that eqs (52), (55), (58), (61) and (64) are in exact agreement with an independent analysis by Friederich (1998), who derived expressions for the first-arriving scattered surface wavefield potentials by simplifying the integral expressions of Woodhouse (1980) for the perturbed potential from specified mode branch–branch coupling by means of generalized spherical harmonics (Phinney & Burridge 1973). After converting Friederich's normalization convention for the mode eigenfunctions to match that of eq. (49), his expressions take the form

where is the second-order Hankel function of degree m, and

Equivalence between eqs (52), (55), (58), (61) and (64) and eq. (66) is established by noting

and the direct relation (eq. C12) between G_m and the Hankel functions. Eqs (68) appear in equivalent form in Appendix D of Friederich (1998).

The preceding analysis demonstrates a subtle connection between the theory of scattering from a small spherical inclusion developed in Section 3 and the theory of coupled normal modes. The scattered wavefield generated by interaction with the spherical inclusion arises from the contrast in material properties between the inclusion and the background medium. To first order, this contrast represents, more precisely, a perturbation with respect to a 1-D spherically stratified reference model. For a sufficiently small single inclusion, the scattered wavefield may be represented by equivalent seismic sources M(ω) and F₀(ω). Solution of eq. (10) in the frequency domain with these equivalent sources yields the scattered displacement field u(r, ω) throughout the earth. The asymptotic forms for the displacement spectrum (eqs C10 and C17) employ a sum of normal modes to solve eq. (10) in a spherically symmetric reference earth model using M(ω) and F₀(ω) (only eq. C10 was actually needed in this section, because we examined only scattered spheroidal modes). When this normal-mode sum is used to represent u(r, ω), the result is equivalent to normal-mode dispersion branch coupling theory in a laterally heterogeneous earth. This connection, demonstrated for small perturbations, further suggests a natural extension to strong contrasts for which normal-mode coupling theory based on small perturbations is inadequate. The entire development undertaken in this section could be repeated by simply substituting the wavefield expansion coefficients (eqs B16 and B17) into eq. (32), evaluating the general form of the source equivalents (eqs 36–41) for strong contrasts with the scattering coefficients of eq. (32), using them to derive the equivalent source excitation functions Σ_j defined in eqs (C4) and (C16), and then using these to represent the scattered wavefield in terms of the potentials given by eqs (C11) and (C18) on spheroidal and toroidal mode dispersion branches, respectively. This could form the basis for efficient multiple-scattering calculations to handle distributed strong lateral heterogeneities in the crust, extending similar methods designed for global scattering calculations in the presence of small perturbations (e.g. Friederich 1998).

The preceding analysis validates all of the source equivalents listed in Table 1. It is also possible to verify the linearized source equivalents by substituting a spherical wavefield of the form of eq. (1) into eqs (4) and (5) of Wu & Aki (1985). Their eq. (4) is identical to the force elements listed in Table 1.

Their eq. (15) reads

where the e_tk denote the elements of the strain tensor. For an incident spheroidal-motion wavefield we have

and for an incident toroidal-motion wavefield we have

As an example, the moment tensor component M₁₂ in the geometry of Fig. 2 is evaluated as

Substituting into eq. (72) the strain component from eq. (8.36) of Aki & Richards (1980),

yields the M₁₂ components in Table 1.

5 Conclusions

The problem of elastic wave scattering from a spherical inclusion has been reconsidered in order to give a complete mathematical description of spherical-elastic-wave interaction with an embedded spherical inclusion in the small-volume limit. This limit satisfies the conditions for Rayleigh scattering, in which the quantity ka′ is much smaller than unity, where k is the wavenumber and a′ is the radius of the spherical inclusion. The treatment yields expressions for the scattered wavefield in terms of (1) vector spherical harmonics centred on the inclusion, which provide a locally valid description even when the inclusion is embedded in a stratified and bounded medium. and (2) equivalent sources. which provide a globally valid description of the scattered wavefield. The treatment does not require any assumption about the contrast in material properties between the inclusion and the background medium, the only exception being the case of interaction with a fluid-filled inclusion. The scattered wavefield in the vicinity of the inclusion involves only spheroidal motion of degrees 0, 1 and 2 m inclusion-centred spherical coordinates. From the nine scattering coefficients describing these vector spherical harmonics, source equivalents (six moment tensor elements and three single forces) are derived. The treatment presented here agrees exactly with that given by Gritto . (1995) for the case of plane P waves incident on a spherical inclusion embedded in an infinite homogeneous medium, and it is further adaptable to applications involving incident S waves or surface waves. To first order in material property contrasts, the treatment as applied to seismic surface waves incident on a small inclusion embedded in a spherically stratified background medium is equivalent to surface-wave scattering theory based on coupling of normal-mode dispersion branches. In addition to yielding another perspective on the physics of surface-wave scattering, it has a high potential for facilitating scattering calculations of body waves for applications in global and regional tomography. As a first step, the first-order source equivalents summarized in Table 1 could form the basis for calculating global scattered wavefields from an embedded heterogeneity using the Born approximation, complementing related efforts (Dalkolmo & Friederich 1996). Further applications could take advantage of the fact that the scattering coefficients are valid for strong material property contrasts. The source equivalents in the general form (eqs 36–41) could be combined with the normal-mode methods of Section 4.2 to yield a ‘kind of’ normal-mode approach in which the assumption of small perturbations has been dropped.

Acknowledgments

This work has been stimulated by discussions with Wolfgang Friederich and Roland Gritto. Comments from them and John Hudson are gratefully acknowledged. This paper benefited from constructive criticisms by Gerhard Müller and an anonymous reviewer.

References

APPENDIX A: DISPLACEMENT-STRESS VECTORS AND THEIR NEAR-FIELD ASYMPTOTIC PROPERTIES

Here we give expressions for the displacement–stress vectors and , defined in eqs (11) and (12), which represent solutions of the equation of motion (10) in a layer where the material properties are taken as constant. Thus we assume that the shear velocity , compressional velocity and density ρ are constant in a given layer. Since these solutions depend explicitly only on l we replace the subscript k with l below. Eqs (6.2)–(6.4) of Ben-Menahem & Singh (1981) yield two independent solutions for toroidal motion,

and four independent solutions for spheroidal motion,

In these equations, is a spherical Bessel function of the first kind and is a spherical Hankel function of the second kind. A prime denotes differentiation with respect to the argument, k_α=ω/α and k_β=ω/β. The ‘+’ solution is regular at the origin r= 0, and the ‘−’ solution represents outgoing waves.

As r (or ω) tends towards zero, then k_βr« 1 and k_αr« 1. We may utilize asymptotic expressions for the Bessel functions in order to obtain small-r expansions of the displacement–stress vectors given in eqs (A1)–(A3). In the case of spheroidal motion, however, in order to retain linearly independent solutions for small r we must take an appropriate linear combination of these vectors:

The asymptotic expansions of the spherical Bessel functions of the first and second kinds, j_l(x) and n_l(x) respectively, for |x| « 1, are

where (2l+ 1)!!= (2l+ 1) × (2l− 1) × (2l− 3) ×… 3 × 1. With and , we obtain the following asymptotic expansions of the spheroidal-motion displacement–stress vectors and , in the limits k_βr« 1 and k_αr« 1:

For the case l= 1, we require additionally the functions resulting when the rigid-translation component is removed from or :

Note that if k denotes the pair (l, m) then may be replaced with defined above, since the do not depend on m.

Asymptotic expressions for the toroidal-motion displacement–stress vectors may be derived similarly from eqs (A1) and (A5) under the condition k_βr« 1:

APPENDIX B: EXPANSION COEFFICIENTS OF INCIDENT SPHERICAL WAVEFIELD IN SCATTERER-CENTRED COORDINATES

For a given incident spherical wavefield defined in the geographic (r, θ, φ) coordinate system, we derive here explicit expressions for the spherical-harmonic components of this wavefield in the scatterer-centred (r′θ′, φ′) coordinate system.

The description of the incident wavefield in scatterer-centred coordinates is necessary to formulate the boundary value problem for the scattered wavefield beginning with eqs (26) and (27).

At fixed frequency ω the displacement field associated with an incident wavefield in a r–θ–φ spherical coordinate system can generally be written in the form prescribed by eq. (8):

with a corresponding traction

For a seismic source located at the pole of the spherical coordinate system (θ= 0) in a laterally homogeneous earth, the m summation goes from −2 to 2. The superscript S or T above denotes either the spheroidal or the toroidal wavefield. In the following we wish to devote attention to a single component of this summation and will work with a wavefield of the form

for spheroidal motion, or

for toroidal motion, such that

If the wavenumber k is a half-integer, then eq. (B3) or (B4) represents simply one component of the summation (B1). Appendix C shows that surface waves can always be written as a superposition of wavefields of the form (B3) and (B4), and incident surface wavefields of more complicated character can, for a single dispersion branch, always be written in those simple forms (Friederich, Wielandt & Stange 1994).

Let a small heterogeneous region be located at r₀= (r₀, θ₀, φ₀) in a coordinate system in which corresponds to the pole θ= 0. We wish to obtain the equivalent displacement-stress vectors in r′= (r′, θ′, φ′) coordinates centred on r₀. We may write

In eqs (B6) and (B7), the subscripts 1, 2, 3 and 4 denote evaluation of the corresponding element of the given column matrix, is given by eq. (25) and denotes the pair (ν, μ). For a scattering analysis in the small-volume limit, we require the expansion coefficients with subscripts (0,0), (1,0), (1, ± 1), (2,0), (2, ±1) and (2, ± 2). In the following we note that the arc φ′= 0 is defined to be directed along the minor arc connecting and in the direction away from (Fig. 2). Define the surface integral

which is to be evaluated on the sphere of radius a′ centred on r₀. Only the coefficient combinations on the right-hand side of eqs (32) are required for further analysis. From the asymptotic expressions for the various y and z vectors which appear in Appendix A, combined with the orthogonality relations for the spherical harmonics, it can be shown that the required coefficient combinations are given by

The integrals which appear in the numerators of eqs (B9)–(B11) may be evaluated by taking the Taylor expansion

For an incident spheroidal-mode wavefield of the form (B3), the derivatives which appear in eq. (B12) are

The orientation of the small sphere centred on r_n has been chosen such that

Using the values for the spherical harmonics

and substitution of eqs (B12)(B15) into eqs (B9)(B11). we obtain, after some algebra, the following coefficient combinations for spheroidal-mode incident wavefields:

A similar analysis for an incident toroidal-mode wavefield of the form (B4) yields the following coefficient combinations:

where

APPENDIX C: UNIFORMLY VALID ASYMPTOTIC EVALUATION OF A SUM OF NORMAL MODES

Here we seek a uniformly valid expression for the seismic displacement spectrum in terms of a sum over normal-mode dispersion branches by decomposing the wavefield into a sum of travelling waves and then analytically evaluating the contribution of the first-arriving surface waves. We shall provide a complete derivation for the case of spheroidal modes and then quote the result for toroidal modes.

Referring to the geometry of Fig. 2, we shall assume a seismic source with moment tensor spectrum M(ω) and force vector spectrum F₀(ω) located at r₀ with spherical coordinates (r=r₀, Δ= 0). The displacement spectrum for spheroidal modes in a laterally homogeneous, non-rotating earth model can generally be written in the form

where (Dahlan 1980)

In eq. (C3)ω_sl and α_sl are the frequency and attenuation parameters, respectively, of the spheroidal mode of degree l on the Sth dispersion branch. The Σ functions in eq. (C2) are given by eq. (13) of Dahlen (1980) for excitation by a moment tensor; when modified to include excitation by a single force they read

where and the functions and are taken to depend implicitly on and ω. In employing these formulae it should be borne in mind that they depend in the present application on the source spectra M(ω) and F₀(ω) [for impulsive sources of the form M(t) =MH(t) or F₀(t) =F₀H(t) the two forms are related by M(ω) =M/(iω) and F₀(ω) = F₀/(iω)]. Application of the Poisson sum formula (Aki & Richards 1980) to (C1) yields

For the first-arriving surface waves we retain only the term s= 0 in eq. (C5):

In order to obtain a uniformly valid asymptotic expression which is valid in both the near field and far field it is expedient to employ an expression for which is uniformly valid in the interval 0 < δ > π–ε, where vε» 1. This is provided by Hub's formula (eq. 60 of Dahlen 1960):

where J_m is the Bessel function of order m. We now decompose J_m into Hankel functions via and further note the relation Combining these properties with the fact that and are even functions of v and the product is an odd (m= 0,2) or even (m= 1) function of v allows us to rewrite eq. (C6) as

Owing to the form of the resonance function , the right-hand side of eq. (C8) has poles in the second and fourth quadrants of the complex plane at

where l^S > 0 denotes the degree of that mode on the Sth dispersion branch on the equivalent non-dissipative reference earth model such that . To prove this we note that on the dissipative earth at fixed l_S the poles of u(r, ω) move from to . We are situated on the dispersion branch at (ω, v_S) rather than , and we therefore have , where U_S is the group velocity at frequency ω on mode branch S. Thus . Since we evaluate eq. (C8) by completing a contour in the lower half of the complex plane, only poles at contribute residues. The final result is

where C_S=ω/v_S is the horizontal phase velocity and

To obtain the corresponding results for toroidal modes we may write the displacement spectrum for toroidal modes in a laterally homogeneous non-rotating earth model in the form

where (Dahlen 1980)

In eq. (C15), ω_Tl and α_Tl, are the frequency and attenuation parameters, respectively, of the toroidal mode of degree l on the Tth dispersion branch. The Σ functions in eq. (C14) are given by eq. (13) of Dahlen (1980) for excitation by a moment tensor: when modified to include excitation by a single force they read

where The uniformly valid asymptotic displacement spectrum is

where l_t is the degree of the mode on the Tth toroidal-mode branch at frequency ω, and v_T is defined by , where is the group velocity at frequency ω on mode branch T).

It is useful to note that the function G_m in eq. (C12) can also be written in the form

where G(v, Δ) is the Green's function for the scattering of first-arriving propagating surface waves on a spherical membrane derived in Appendix B of Pollitz (1994). If the far-field expression for G(v, Δ) given by eq. (B10) of Pollitz (1994) is substituted into eq. (C19) and then into eqs (C10) and (C11) or (C17) and (C18), then the result is identical to the complex conjugate of the far-field result given in eq. (19) of Snieder & Nolet (1987) for excitation by a single-force spectrum (the complex conjugate arises from the difference in Fourier transform convention between that paper and the present paper).