6(r3

TR diss

2034

DIRECTIONAL DECOMPOSITION OF TRANSIENT ACOUSTIC WAVE FIELDS

DIRECTIONAL DECOMPOSITION OF TRANSIENT ACOUSTIC WAVE FIELDS

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, prof. drs. P.A. Schenck, in het openbaar te verdedigen ten overstaan van een commissie aangewezen door het College van Dekanen op donderdag 27 februari 1992 te 14.00 uur

door

Maarten Valentijn de Hoop

geboren te Delft
...

natuurkundig doctorandus

Delft University Press

Dit proefschrift is goedgekeurd door de promotoren prof. dr. ir. A.T. de Hoop en prof. dr. ir. P.M. van den Berg

toegevoegd promotor: dr. ir. J.T. Fokkema

Published by Delft University Press Stevinweg 1 2628CN Delft the Netherlands

ISBN 90-6275-747-2/CIP

Copyright @ 1992 by M.V. de Hoop

All rights reserved. No parts of this publication may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the permission of the Koninklijke/Shell Exploratie en Produktie Laboratorium.

In memory of my mother

The author is employed by the Koninklijke/Shell Exploratie en Produktie Laboratorium, Rijswijk, the Netherlands, and gratefully acknowledges the support for the research reported here and the permission to publish this monograph.

Contents

1 Introduction 2 Scalar space-time waves in their spectral-domain second-order Thiele approximations 2.1. Introduction . first- and

19

20
23 27 37

2.2. The exact wave equation and its solution in the space-time and spectral domains . 2.3. The first-order Thiele approximation 2.4. The second-order Thiele approximation 2.5. Arrival times as functions of horizontal offset . 2.6. The local behavior of the waves near the wave fronts in the vertical direction .. 2.7. Concluding remarks.

51
56 58

3 Propagation of acoustic waves in fluid media in the first- and second-order Thiele approximations 61 3.1. Introduction. ........ 62 63

3.2. The wave-matrix formalism

3.3. The transform-domain solution for an unbounded homogeneous fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

3.4. The exact solutions in the space-time domain 3.5. The first-order Thiele-approximated solutions in the space-time domain

· 76 76

3.6. The second-order Thiele-approximated solutions in the spacetime domain . . . . . . . . . . . . . . . . . . . . · 83 3.7. Exact, and first- and second-order Thiele-approximated reflection and transmission coefficients ... . . . . . . . . . . . .. 3.8. Exact, and first- and second-order Thiele-approximated reflected waves in a two-medium configuration 3.9. Discussion of the results . 91 97 .119

Propagation of acoustic waves in an elastic solid in the parabolic approximation 121 4.1. Introduction 4.2. The elastodynamic wave-matrix formalism 4.3. The parabolic approximation. ... · 122 · 124 · 138 · 141

4.4. Jump conditions at the source level

4.5. The transform-domain particle velocity in an unbounded homogeneous solid due to a point force , . 143 4.6. The elastodynamic Green's tensor in the space-time domain 4.7. The tensorial "one-way" wave operator 4.8. Discussion of the results Directional ing media ........ in smoothly . 145 . 151 . 157 vary159 .160 .162

acoustic wave field decomposition

5.1. Introduction 5.2. Decomposition of the acoustic scattering process

5.3. The dimensionally reduced scattering problem

.... . . . . . · 169

5.4. The system of one-way wave equations in the time-Laplace. . transform domain . .... .177

..

5.5. The generalized vertical slowness 5.6. The Green's functions of the one-way wave operators 5.7. The Bremmer coupling series 5.8. Evaluation of the symbols 5.9. Discussion of the results 6 Iterative acoustic inverse scattering theorem of time-correlation type 6.1. Introduction. .......... based on the reciprocity

.179 .183 · 189 · 193 · 198

205 .205

6.2. Operator formalism of the basic equations of the acoustic wave problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 6.3. Contrast-source formulation of the inverse scattering problem .. 210 6.4. Reciprocity theorems and the weak formulation of the inverse problem . . . . . . . . . . . . . . . . . . 211 6.5. Green's functions and their properties. 6.6. Inverse scattering .. . 219 . 221 . 227 . 229 approximations 231

6.7. The iterative scheme 6.8. Discussion . . . . . . A Thiele's continued-fraction

B Closed-form expressions for the second-order modified Cagniard contour 243 C Interrelation between the wave amplitudes at an interface with

surface sources D The first-order Thiele-approximated relation to the retarded potential E Evaluation of the first-order namic Green's tensor F Calculus of pseudo-differential Green's function and its

251

255 elastody259

Thiele-approximated

operators

263 269 273 287

G An inversion scheme with variable source distributions Bibliography Samenvatting

Chapter 1 Introduction
Directional decomposition of acoustic wave fields is a mathematical tool to analyze seismic field data in exploration geophysics [1,2,3]. In such a decomposition, a direction of preference in which the subsurface properties in the geological framework are assumed to vary more rapidly than in the family of planes perpendicular to it, plays a crucial role. The latter property reflects itself in a certain manner in the structure of the acoustic wave fields that carry the information about the distribution of acoustic properties of the subsurface. The extraction of the latter distribution from the seismic field data is in itself an inverse scattering problem. The purpose of the present thesis is to analyze the interrelation that exists between different degrees of approximation carried out in the directional decomposition procedures. Each of the approximations leads to an accompanying reduction in the number of operations to be carried out in the pertaining numerical algorithm. To place the different aspects in the whole of seismic processing, the directional decomposition is linked to a novel view on seismic inverse scattering theory and depth migration. Each aspect is covered in a separate chapter. The chapters are written such that they can be read independently; the interrelation between the different chapters is indicated below.

INTRODUCTION

Inverse scattering and depth migration

The ultimate purpose is to incorporate the directional wave field decomposition into an (approximate) seismic inverse scattering procedure (depth migration). Conceptually, depth migration algorithms in seismic prospecting, and reconstruction algorithms in acoustic imaging in general, can be envisaged as a "remote sensing problem" where, through the intermediate of judiciously chosen computational methods, the physical properties of a remote part of the configuration in which the probing (in our case acoustic) waves propagate are to be reconstructed from the measured data. Mathematically, such a type of interaction is described by a reciprocity theorem [4, 5, 6, 7, 8] through which different acoustic states are interrelated. In acoustic wave motion there are two reciprocity theorems: one of the time-convolution type and one of the time-correlation type. The reciprocity theorem of the time-convolution type leads, via the introduction of point-source solutions (Green's functions) to the wave problem, to integral representations for the quantities that characterize the wave motion. These representations are the ones that Kirchhoff introduced in the theory of light waves and form the basis for the "Kirchhoff" depth migration algorithms [9, 10, 11, 12]. In Chapter 6, the properties of the reciprocity relation of the time-correlation type are investigated with a view of their giving rise to a different type of reconstruction algorithm. In its linearized version one can recognize from it the algorithms that are known under the names "wave-equation" or "finitedifference" [1, 2, 3, 13] and "phase-shift-plus interpolation" [14, 15J pre-stack depth migration, while the observations needed on the boundary enclosing the unknown objects are taken into account and replaced by computational fields. An important aspect of taking the reciprocity relation of the time-correlation type as the point of departure is that it introduces a proper norm in the function space of the pertaining wave fields and thus leads, via the iterative minimization according to the corresponding error criterion, to a "minimumnorm" gradient technique with a decreasing error in the data fit [16], which technique (locally) accounts for the non-linearity of the inverse problem in the medium parameters. The general aspects of the latter method are discussed,

to which end an appropriate operator formalism is introduced. Through the parametrization of the medium, a well-defined hierarchy of its complexity is incorporated in the successive steps of the scheme. In the process, the wave fields in the relevant states are directionally decomposed. Employing this decomposition, the computational effort in the numerical calculations is at the same time reduced.

Wave field decomposition

In Chapter 5 an operator formalism is developed to expand the transient acoustic wave field in a smoothly varying medium, generated by a source localized in space and time, into a superposition of constituents each of which can be interpreted as having traveled up and down with respect to a direction of preference a definite number of times. This is a generalization of the Bremmer coupling series [17J. Both the existence and the convergence of the relevant series are discussed. The pseudo-differential operator calculus that is the appropriate tool for this leads to a natural generalization of the concept of slowness surface to multi-dimensionally smoothly varying media. Let the "vertical" direction be the preferred one, then the operator associated with the corresponding generalized vertical slowness induces the full "one-way" wave operator. The latter is loosely denoted as the "parabolic" wave operator [1, 18J. In its exact form the decomposition reduces the scattering problem in n dimensions to a family of scattering problems in n - 1-dimensional hyperplanes, after which the corresponding results are used in a Neumann series to solve the remaining scattering phenomenon in the direction of preference. The decomposition introduced applies to media with smoothly varying properties. It is not known whether the operators involved can be defined in media with discontinuous properties. The theory leads to integro-differential equations that generate the different constituents (terms in the series). The kernel associated with the relevant integral operator can be approximated in such a way that the matrix representation in the corresponding numerical scheme becomes sparse.

----------

..

--

INTRODUCTION

Compression of the relevant operators in a fluid medium

The relevant operators in the directional decomposition procedure can, in principle, be compressed with the aid of a multi-scale transformation [19]. In Chapters 2-4 a two-step approach is followed, viz., a rational expansion of the

_, horizontal offset

Figure 1.1: A two-dimensional salt model with compressional wave speed as a function of position. relevant operators in terms of powers of the horizontal-space Laplacian is considered; the sparse matrix representations of standard elliptic operators like the

SSP
->

horizontal offset

sC" ~
0

.." :e
1

..
.".

0,500

1,000

, .SOQ

Figure 1.2: Snapshot at t = O.750sof the full acoustic pressure due to a vertical point force at {O,O}in the model of Figure 1.1.

Laplacian are well known (e.g., employing finite-difference techniques [20, 21], see Figures 1.1-1.3, or using orthonormal compactly supported wavelet bases

[19]).
Thus, for a simple and sparse matrix representation of the desired so-called "one-way" wave operator we introduce an intermediate partial differential operator. The Thiele approximation, when carried out on the left symbols of the associated pseudo-differential operators about a preferred direction of propagation, yields such a tool. The resulting partial differential equations,

10

INTRODUCTION

SSP
-->

horizontal offset

<J.500

S
~
Ol

~
0

1l
1.000

.~ 1:: :"
\ .500

Figure 1.3: Snapshot at t = 0.750s of the second-order Thiele-approximated acoustic pressure due to a vertical point force at {O,O}in the model of Figure 1.1 (leading term of the Bremmer coupling series).

11

which are denoted as one-way wave equations, are usually solved with the aid of finite-difference techniques carried out in the time Fourier-transform domain. In this respect, it is noted that the stable two-level implicit finite-difference integration of a one-way wave equation along the direction of preference leads to the problem of solving a matrix equation. To fully exploit the freedom in the structure of this matrix equation a rational approximation of the vertical slowness as a function of the horizontal slowness is desired, which is provided by Thiele's formula. Since the directional decomposition does not separate waves that travel from right to left and vice versa in the horizontal plane, the theory is also able to handle guided waves along the direction of preference. Figures 1.4-1.8 illustrate this. The Thiele expansion of the left symbol is a way to approximate the operator for interactions in the family of horizontal planes inside a cone centered on the axis of preference, and as such differs from the standard asymptotic expansion of symbols [22]. In Chapter 2, for the simple case of the scalar wave motion generated by a point source in an unbounded homogeneous medium, it is investigated what the consequences of the first- and second-order Thiele approximations of the vertical slowness in the spectral domain, i.e., after having transformed the actual space-time wave motion to the time Laplace transform domain and the spatial Fourier transform domain as far as the horizontal coordinates are concerned, are in the space-time domain. To this end, the corresponding spectral-domain Green's function is transformed back to the space-time domain with the aid of the modified Cagniard method [23]. The exact solution to the problem is a spherical wave with the same wave shape as the source signature and a single wave front. The first-order Thiele, or parabolic, approximation has a wave front in the shape of a double oblate spheroid, combined with a headwave-like precursor having a cylindrical wave front. The second-order Thiele approximation contains two wave fronts: one associated with a fast body wave and the other with a slow body wave, in combination with a head-wave-like precursor, the latter again having a cylindrical wave front. The first- and second-order Thiele approximants are basic III the sense that any odd- or even-order approximation can be written as a summation of

12

INTRODUCTION

Figure 1.4: A three-dimensional salt model (range in the horizontal directions is 2.048kmj range in the vertical direction is 1.6km).

13

SSP

0.000

s :':'.,
0

~ ~
'"il .~
t:

(J.500

~
1
1.000

1.500

Figure 1.5: A vertical slice of the model of Figure 1.4.

14

INTRODUCTION

s c

.. ~
~
0

0,500

~
1
1.000

Figure 1.6: The same vertical slice as in Figure 1.5 of the snapshot at t = Is of the second-order Thiele-approximated acoustic pressure due to a vertical point force at {O,0, O}in the model of Figure 1.4 (leading term of the Bremmer coupling series).

15

Figure 1.7: The same vertical slice as in Figure 1.6 of the wave front at t evaluated through ray tracing in case the salt dome is non-reflecting.

= Is

16

INTRODUCTION

Figure 1.8: The same vertical slice as in Figure 1.6 of the wave front at t = Is evaluated through ray tracing in case the salt dome is reflecting.

17
terms that resemble the low-order approximants. This implies that the waves associated with the higher-order approximations applied in a homogeneous medium are equivalent to the waves in a quasi-multilayer configuration in which the vertical slowness in each individual layer is given by a first- or second-order-like approximant. In view of the applications, special attention is paid to a comparison with asymptotic ray theory along the vertical direction. The arrival-time curves and the snapshots show that with the Thiele expansion the accuracy of asymptotic ray theory is achieved along the vertical ray only; in oblique directions of propagation there is a discrepancy. In oblique directions a better approximation is found if the vertical slowness is expanded about the actual direction of the ray trajectory connecting the source and receiver. Thus, Thiele approximations reproduce amplitudes along the "vertical ray" at least as accurately as the asymptotic ray theory does. However, in contrast to asymptotic ray theory, the Thiele-approximated wave theory has the advantage that it automatically initializes the ray amplitudes in the immediate neighborhood of the source. In Chapter 3 the acoustic reflection problem is analyzed with regard to the approximations of the reflection and transmission coefficients and the physical head wave. The analysis reveals the general consequences that spectral-domain Thiele approximations for the vertical wave slowness have for the corresponding space-time domain wave motion. The most striking feature is that Thiele approximations lead to a number of non-physical artefacts that do not occur in the exact wave motion. First of all, for even-order approximants, a singularity in the wave motion is observed right above and below the source, which is absent in both the exact wave motion and in odd-order approximants. Further, for horizontal offsets that are large compared to the vertical offset, non-physical head-wave-like precursors occur. In addition, any Thiele approximant of order higher than the first introduces "slow" waves, which again are absent in the exact wave motion. In the odd-order approximations, however, they never occur along a direction of propagation which lies interior to a certain cone centered on the vertical axis, which is a desired feature. At this stage, the third-order Thiele approximation has been analyzed as well [24] and indeed, inside a cone centered on the direction of preference no artificial arrivals do

18

INTRODUCTION

occur. On the other hand, the physical head wave cannot be approximated with the aid of the Thiele expansion. Through an adapted form of Thiele's formula the theory can be extended to anisotropic fluids. This is discussed in Appendix A.

Compression of the relevant operators in an isotropic solid

In Chapter 4, the directional decomposition and the associated compressions are extended to perfectly elastic isotropic solids. A fundamental feature of the latter is that the interaction of waves that travel up and down with respect to the direction of preference must be separated from the interaction between polarization states. The extension is of particular importance in the case of a medium with localized high contrasts in the wave speeds, where the conversions between compressional and shear waves are significant.

Concluding remark
The main purpose of the thesis is to arrange in coherence and investigate in depth the underlying mathematical structure (Chapters 5 and 6) of seismic inversion and to show in relatively simple cases (Chapters 2, 3 and 4) what consequences particular approximate techniques have for the interpretation of seismic data.

19

Chapter 2 Scalar space-time waves in their spectral-domain first- and second-order Thiele approximations
For the simple case of the scalar wave motion generated by a point source in an unbounded homogeneous medium, it is investigated what the consequences of the first- and second-order Thiele approximations in the spectral domain are in the space-time domain. To this end, the corresponding spectral domain Green's function is transformed back to the space-time domain with the aid of the modified Cagniard method. The exact solution to the problem is a spherical wave with the same wave shape as the source signature and a single wave front. The first-order Thiele, or parabolic, approximation has a wave front in the shape of a double oblate spheroid, combined with a head-wave-like precursor having a cylindrical wave front. The second-order Thiele approximation contains two wave fronts: one associated with a fast body wave and the other with a slow body wave, in combination with a head-wave-like precursor, the latter again having a cylindrical wave front. From the results, it can be concluded in which regions below or above the source the approximations have sufficient accuracy for their application in inversion and "true amplitude"

20

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

depth migration procedures in geophysical prospecting.

2.1.

Introduction

The first-order Thiele continued-fraction (or parabolic) approximation was introduced in the analysis of wave phenomena by Leontovich and Fock [25];their field 'of application was the propagation of electromagnetic waves close to the surface of a convex electrically conducting body [26]. Higher-order approximations were considered by Bremmer [27]in the asymptotic evaluation of the integrals that occur in the theory of diffraction of waves by an aperture in a screen. Since then, the parabolic approximation has been applied to various wave problems, amongst which are the propagation of waves in random media [28], the downward continuation of acoustic waves in seismic prospecting [1], and underwater acoustics [18, 29]. In recent years there is an increasing interest in the use of higher-order Thiele-type continued-fraction approximations (which include the parabolic and Pade approximations) in the application of seismic modeling, migration and inversion techniques. We mention the papers by Ristow [30], Ma [31], Lee and Suh [32],Kelamis and Kjartansson [33],and Graves and Clayton [34]. The motivation of the Thiele approximation is the following. For one-dimensional wave propagation it is known that the wave operator can be factorized into two partial differential operators that can be identified as to apply to waves propagating in opposite directions. For wave propagation in more dimensions such a factorization can again be carried out, be it that the factors contain pseudo-differential operators. The latter have no simple discrete counterparts that are sparse in their matrix representation. For a simple and sparse matrix representation of the desired so-called "one-way" wave operator we need an intermediate partial differential operator. The Thiele approximation, when carried out on the slowness in the spectral domain, i.e., after having transformed the actual space-time wave motion to the time Fourier transform (or time Laplace transform) domain and the spatial Fourier transform domain as

INTRODUCTION

21

far as the horizontal coordinates are concerned, about a preferred direction of propagation, yields such a tool. The preferred direction will be identified as the vertical direction. (A Taylor-type approximation leading to such a factorization has been discussed by Corones and Krueger [35].) The resulting partial differential equations, which are denoted as one-way wave equations, are usually solved with the aid of finite-difference techniques carried out in the time Fourier-transform domain [34]. In this respect, it is noted that the stable two-level implicit finite-difference integration of a one-way wave equation along the direction of preference leads to the problem of solving a matrix equation. To fully exploit the freedom in the structure of this matrix equation, a rational approximation of the vertical slowness as a function of the horizontal slowness is needed, which is provided by Thiele's formula. Thiele approximations have been either arrived at with the aid of a perturbation theory applied to the wave equation in a comoving frame of reference [33] or carried out on the wave function in the spectral domain [36, 37]. The properties of the waves in the first-order approximation have, to a certain extent, been studied by Joly [38] and Bamberger et al. [39]; higher-order approximations have been considered by Bamberger et al. [40]. Deviations, due to the approximation, in the phase and group velocities of normal modes in underwater acoustic ducts have been analyzed by McDaniel [29]. Approximations of the vertical slowness were further employed to obtain asymptotic expressions in space-time for the wave constituents associated with the conical points in an anisotropic elastic medium having cubic symmetry [41]. Although the one-way wave equations are strictly valid only in a horizontally shift-invariant medium, they have nonetheless been applied in arbitrarily inhomogeneous media. Viewed from this perspective, a Thiele approximation is recognized as a high-frequency approximation that does not suffer from the singularities at which asymptotic ray theory fails (caustics and foci). Proper extensions of the one-way wave equations to fully inhomogeneous media have been discussed by Palmer [42],who also includes a source term in the analysis, and by Fishman and McCoy [43], who employ pseudo-differential operators to account for the mutual coupling between the different spatial Fourier components. A related approach has also been developed by Gazdag [14]and Gazdag

22

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

and Sguazzero [15]. The Thiele approximants have, in the spectral domain, been compared with approximants of the vertical slowness obtained through least-squares fitting by Lee and Suh [32],through interpolation at a number of values of the horizontal slowness by Halpern and Trefethen [44], and through the replacement of it by an integral representation by Zhang Guan-quan et al. [45]. The relation between the solutions of the exact (Helmholtz) and the firstorder approximated (Schrodinger) equations for propagation of sound in an acoustic wave guide has been analyzed by Polyanskii [46] and De Santo [47]to investigate what consequences the approximation of the equation has for the solutions. The purpose of the present chapter is to investigate, for the simple case of the scalar wave motion generated by a point source in an unbounded homogeneous medium, what the consequences of the approximations are in the space-time domain. In this case all calculations can be carried out exactly in closed form. We determine the first- and second-order Thiele approximations of the wave function in the spectral domain and transform the results back to the space-time domain with the aid of the modified Cagniard method as it has been developed by A.T. de Hoop (cf. [23, 48] and [49, pages 298301], [50, pages 302-314], [51, pages 243-253] and [52]). (An alternative but equivalent method has been given by Petrowsky [53].) The exact solution to the problem is a spherical wave having the same shape as the source signature and a single wave front. The first-order Thiele approximation turns out to have a composite wave front consisting of an oblate spheroid associated with the body wave, in combination with a head-wave-like precursor having a cylindrical wave front. The second-order Thiele approximation turns out to have two non-spherical wave fronts associated with two body waves containing cusps, again accompanied by a head-wave-like precursor having a cylindrical wave front. The results enable one to judge in which region of space-time the relevant approximations are accurate enough for use in seismic inversion or migration procedures. In view of these applications, special attention is paid to a comparison with asymptotic ray theory along the vertical direction.

EXACT

SOLUTION

23

The first- and second-order Thiele approximants are basic in the sense that any odd- or even-order approximation can be written as a summation of terms that resemble the low-order approximants. This implies that the waves associated with the higher-order approximations applied in a homogeneous medium are equivalent to the waves in a quasi-multilayer configuration in which the vertical slowness in each individual layer is given by a first- or second-order-like approximant. With the aid of the general form of Thiele's formula the theory can be extended to anisotropic fluids (Appendix A).

2.2.

The exact wave equation and its solution in the space-time and spectral domains

To locate a point in space we use the coordinates {Zl' Z2, Z3} with respect to the orthogonal Cartesian reference frame with origin eJ and the three mutually perpendicular base vectors {i1' i2, i3} of unit length each. In the indicated order, the base vectors form a right-handed system. The subscript notation for Cartesian vectors and tensors applies and for repeated subscripts the summation convention holds. Lowercase Latin subscripts are used for this purpose; they are to be assigned the values {I, 2, 3}. The position vector corresponding to {Xl, X2, X3} is therefore denoted by :Il = Zm im• The time coordinate is denoted by t. The symbol 8m indicates differentiation with respect to Xm; 8t is a reserved symbol for differentiation with respect to time. The scalar wave function u = u(:Il, t) that is representative of the wave motion generated by a point source located at {O,0, O} and with source signature f = f( t) satisfies the three-dimensional scalar wave equation

(2.1)
where c is the wave speed in the medium and 6(Xl, X2, X3) denotes the threedimensional unit impulse (Dirac distribution) operative at the point {O,0, O}.

24

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

The exact solution to this equation is the spherical wave [49] u = J(t -lzl/c}/47rlzl, (2.2)

where Izl = (ZmZm)1/2. Its wave shape is the same as the source signature. Equation (2.2) can also be written as u(z,t)
= 8t

Jtl=lzl/c

J(t - t')G(z,t')dt',

(2.3)

where the source is assumed to start acting at the instant t = 0, and in which G(z, t) = H(t - Izl/c)/47rlzl, (2.4)

where H(t) = {0,1/2,1} for {t < O,t = O,t > O} denotes the Heaviside unit step function. Figure 2.1 shows a snapshot in space of this Green's function. How the well-known result of Eq.(2.2) can be obtained with the aid of the modified Cagniard method is shown in [23]. For this, as well as for the Thiele approximations of this chapter's subsequent sections, the spectral domain representation of u is needed. The latter follows by applying to Eq.(2.1) the one-sided Laplace transformation

u(z,s)

roo i.; exp(-st)u(z,t)dt,

(2.5)

followed by the spatial Fourier transformation with respect to the "horizontal" coordinates Z1 and Z2,

where Greek subscripts are used to indicate the horizontal components of vectors and tensors; to them the values {1,2} are to be assigned. Further, s is chosen to be real and positive, and al and a2 are real. The factor s in Eq.(2.6) is introduced for convenience in the transformation back to the space-time domain. The transformation inverse to Eq.(2.6) is given by

EXACT

SOLUTION

25

c = 750m/s ; t = 0.450s
SSP
source location 0.0

500.0

Figure 2.1: Snapshot of the exact Green's function in space for horizontal offsets are numbered sequentially.)

:1:3 ~

O. (The

26

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

The spectral-domain counterpart equation

u of u then

satisfies the ordinary differential (2.8)

where 8(:1:3) the one-dimensional unit impulse (Dirac distribution) operative is at :1:3 0 and =

(2.9)
is the spectral-domain by "vertical" slowness. The solution of Eq.(2.8) is given (2.10) To bring Eq.(2.10) into accordance with the form of Eq.(2.3), it is rewritten as u = 8j(8)G(ia,., :1:3, 8). (2.11) Comparing Eq.(2.11) with Eq.(2.10), it is found that

G = exp( -8')'1:1:31)/282,),.

(2.12)

Upon substituting Eq.(2.12) in Eq.(2.7), the algebraic factors of 8 cancel and only the 8 in the exponential function remains, which property serves as a criterion for selecting a Green's function in the Cagniard method. By differentiating the right-hand side of Eq.(2.12) twice with respect to :1:3, we can reconstruct the differential equation that G satisfies. In view of the non-analyticity at :1:3 0, this procedure leads to the required volume source = density in its right-hand side. Specifically, we obtain

8iG - 82')'2G =

-8-18(:1:3).

(2.13)

Now, under the transformations of Eqs.(2.5) and (2.6), we have 8t --+ sand 8,. --+ -isa,.. Applying these rules, we reconstruct from Eq.(2.13) the partial differential equation that G satisfies as (2.14) Note that the wave operator on the left-hand side does not admit, in threedimensional space, a factorization into space-time differential operators. (Such a factorization is restricted to only one space dimension.)

FIRST-ORDER

THIELE

APPROXIMATION

27

2.3.

The first-order Thiele approximation

(see Appendix A), I in Eq.(2.9) is (2.15)

In the first-order Thiele approximation replaced by

II

= lie + (e/2)0.,,0.,..

It is noted that this approximation is identical to the first-order Taylor expansion about e20."0.,, = o. The resulting spectral-domain wave function is given by -I -I. (2.16) u = sf(s)G (lo.",:l:3,S),
A

in which, in view of Eq.(2.12), (2.17) To reconstruct the latter's space-time equivalent, we substitute Eq.(2.17) into Eq.(2.7) and apply the modified Cagniard method. This method consists of transforming the variables of integration {all 0.2} with {0.1 E 'R,0.2 E 'R} into {p, q} with {p E I, q E 'R} (here, 'R denotes the real axis and I the imaginary axis) through

0.1 = -ip cos(8) - q sine8),

0.2

(2.18) (2.19)

-ip sine8) + q cos(8),

where {T,8} is related to {:l:ll:l:2} by

:1:1

= Tcos(8),

:1:2

= Tsin(8),

(2.20)

with 0 l' < 00,0 (J < 271", keeping q real, continuing the integrand analytically into the complex p-plane, and carrying out the integration along the modified Cagniard path that follows from a continuous deformation of the imaginary p-axis and satisfies the equation (2.21) where T is real and positive. In view of the condition that the joining paths at infinity must yield a vanishing contribution, only contours in the right half

s:

s:

28

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

of the complex p-plane are considered. Solving for pin Eq.(2.21), we see that the modified Cagniard path consists of p = pI(T,q) in the first quadrant of the complex p-plane, together with its complex conjugate p = pI*( T, q) in the fourth quadrant of the complex p-plane, where

in which (2.23) with (2.24) Equation (2.22) represents a straight line parallel to the imaginary p-axis that intersects the real p-axis at p = r / clz31. Along this path, T strictly increases when going from the point of intersection with the real p-axis to infinity. In the process of contour deformation we might possibly encounter the simple pole p = p~ of {jI at the simple zero of 'YI in the right half of the p-plane (cf. Eq.(2.17)). From Eqs.(2.15), (2.18) and (2.19) it follows that (2.25) This pole contributes to the Green's function when it lies to the left of the point of intersection of the modified Cagniard path with the real p-axis. In the case where r < IZ310 ("small" horizontal offset), the simple pole lies, for all values of q, to the right of the point of intersection of the modified Cagniard path with the real p-axis and is not passed in the process of contour deformation (see Figure 2.2).

The body wave

In the integral along the modified Cagniard path, T is introduced as the variable of integration. The corresponding one-dimensional Jacobian follows from Eq.(2.22) as (2.26)

FIRST-ORDER

THIELE

APPROXIMATION

29

c = lOOOm/s ; r = lOOOm ; X3 = lOOOm 4~-r. __ ~~_,-------.------~:-------'------, .......................................................[ "1"

q=0

~ 3:::::--::=::::::-f:::=:-:r-=I-::--~
.....-.t

2
•..•••••.••.. ••.....•...•....

········:······························t ..·············

1"""

'/;:
--

body ~ave
L _--. ---"'-r-"
!

1 __ ·········1
_

j .•....••.•. __

............................ ................................

\ t!; ······················.·T·· ..··. '

I !!

....--....-t.. ............................ 1" I

i

__.

-1 ····························· ..1····..·················_ ..+ ·

···············;·································1

-2 ..····························1················
............................ · 1 ··

·· ..

·.. f

·· ·..·!····

···

·,

·
·

·+ · · · +
·+ · · ·..·+·

+...................,
· · ·.. ;·. ·· · ·.. · .· 1
+ ,

-3 .. · ·

................................ l·· · ·.. +

-4~~~--~-4·------~------4-------~----~ I

2
-+

Re{pc}

Figure 2.2: The modified Cagniard contour and the pole location for a "small" horizontal offset in the first-order Thiele approximation.

30

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

Taking the contributions from p = pI and p = pl. together (note that the integrand satisfies Schwarz's reflection principle), the contribution from the Cagniard path is arrived at through Eq.(2.7) as

Uc{Zm's)
in which

•I

= sf(s)Gc{zm,s),

••

(2.27)

G~ = 7r-2Im G~ =

Jq=o

(00

dq

00

T=Tl(q)

exp[

_sr)[2-yI(pI,

q)tl(8pI

/8r)dr}.

(2.28)

Interchanging the order of integration yields

00

T=Tl(O)

exp(-sT)dr

(Ql(T) Jq=o

7r-2Im{[2-yI(pI,q)tl(8pI/8r)}dq,
is given by

(2.29)

in which q

= QI(r),

the (unique) inverse of r = TI(q),

QI(T)

= (2/clz31)1/2[r _ TI(0)]1/2.

(2.30)

Application of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function when taken at positive real values of s [54, pages 61-63] leads to the time-domain results G~ = and

[r= 7r-2Im{[2-yI(pI,q)tl(8pIj8r)}dq
Jq=o

] H[t-TI(O)]
t') dt'j

(2.31)

U~(Zm' t)

= 8t

it

I'=TI(O)

J(t - t')G~(Zm'
,

(2.32)

from which it is evident that TI(O), as given by Eq.(2.24), is the arrival time of the wave. The equation for the corresponding wave front follows from Eq.(2.24) as (2.33) which is a double oblate spheroid centered on r = 0 and with semiminor axis cTI(0)j2 and semimajor axis cTI(O)jV2. substitution

IZ31

cTI(0)j2

Through the

q = QI(t)sin(,p),
the integral in Eq.(2.31) can be evaluated in closed form. The result is G~ = 47rC(t2 _ ~r2 /C2)1/2 sgn(t -

(2.34)

r2/clz31)H[t

- TI(O)].

(2.35)

Figure 2.3 shows a snapshot in space of this Green's function for

Z3 ~

O.

FIRST-ORDER

THIELE

APPROXIMATION

31

c = 750m/s ; t = 0.450s

Figure 2.3: Snapshot ofthe Cagniard-path contribution to the Green's function in space for :1:3 ~ 0 in the first-order Thiele approximation. (The horizontal offsets are numbered sequentially.)

32

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

The precursor
When r > IZ31v'2 ("large" horizontal offset), the simple pole P6 lies to the left of the point of intersection of the modified Cagniard path with the real p-axis if q2 < (qJ)2, where qJ = [(r/clz31)2 - 2/c2Jl/2. In this case, the pole's contribution has to be taken into account in the process of contour deformation, as illustrated in Figure 2.4. The relevant contribution is found to be (2.36) which, upon replacing r(q2

+ 2/C2)1/2 by r, takes the form

G"
in which

I __ 0-

1_lTt.2
r=Ttl

21l"c

exp( -sr) d I -r, (r2 - 2r2/c2)I 2

(2.37)

(2.38) On account of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function when taken at positive real values of s, the time domain counterpart G~ of G~then follows as (2.39) This contribution (see Figure 2.5), which is present only if r > 1:l:31v'2, is to be added to the contribution from the modified Cagniard path Eq.(2.35). The corresponding wave front has a cylindrical shape and touches the spheroid of G~. It can be interpreted as a head-wave-like precursor. Figure 2.6 shows a snapshot in space of the total Green's function G~ + G~ for :1:3 ~ o.

The case

x3

=0

In the case :1:3 = 0, the expression in Eq.(2.22) breaks down; from Eq.(2.21) it follows that the modified Cagniard contour coincides with the real, positive p-axis. In view of the analyticity of the integrand away from the pole, the

-----

FIRST-ORDER

THIELE

APPROXIMATION

33

-_ ..............
c

lOOOm/s ;

lOOOm ;

X3

SOOm
, i
....................

....

q=O

./body wave

........ "head wave"

-1

i/

r0-

-2

-3

-4

_"' .. ~'
2
3 4 5

~ Re{pc}
Figure 2.4: The modified Cagniard contour and the pole location for a "large" horizontal offset in the first-order Thiele approximation.

34

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

c = 750m/s ; t = 0.450s
SSP
->

horizontal offset (rn)

source location

Figure 2.5: Snapshot of the pole contribution to the Green's function in space for :1:3 2:: 0 in the first-order Thiele approximation. (The horizontal offsets are numbered sequentially.)

FIRST-ORDER

THIELE

APPROXIMATION

35

c
SSP

= 750m/s

;t

= 0.450s

_, horizontal offset (m) source location 0.0 -448.0

448.0

500.0

Figure 2.6: Snapshot of the total Green's function in space for :1:3 ~ 0 in the first-order Thiele approximation. (The horizontal offsets are numbered sequentially. )

36

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

expression for Eq.(2.36))

Ol

reduces to the contribution from the pole only, i.e., (cf. All GO=211"c

00

q=O

exp[-sr(q2 + 2/C2)1/2J (q2+2/C2)l/2 dq.

(2.40)

The space-time counterpart to this expression is found to be (cf. Eq.(2.39)) (2.41) in which Ttl is given by Eq.(2.38) .

Partial differential equation for G1

By differentiating the right-hand side of Eq.(2.17) twice with respect to X3, and using that under the transformations of Eqs.(2.5) and (2.6) we have 8t --+ S and 8" --+ -isa", we can reconstruct the partial differential equation that Gl satisfies, viz.,

This equation factorizes as

[c-18t(83

+ c- 8t)
l

(1/2)8,,8,,][c-l8t(83

c-l8t)

+ (1/2)8v8vlGl
(2.43)

= -c-28;S(xt,

X2,

x3)H(t).

The first partial differential operator In brackets on the left-hand side of Eq.(2.43) is the parabolic approximation to the wave operator for a wave propagating in the direction of increasing X3; the second partial differential operator in brackets on the left-hand side of Eq.(2.43} is the parabolic approximation to the wave operator for a wave propagating in the direction of decreasing X3

[lJ.

SECOND-ORDER

THIELE

APPROXIMATION

37

2.4.

The second-order Thiele approximation

In the second-order Thiele approximation (see Appendix A), , in Eq.(2.9) is replaced by II 3 4/3c2 + a,..aIA , =. (2.44) c 4/c2+avav The resulting spectral-domain wave function is given by

u
in which

-II

sf(s)G
A

-II'

(zalA,:Z:3,s),

(2.45) (2.46)

To reconstruct the latter's space-time equivalent, we substitute Eq.(2.46) into Eq.(2.7) and apply the modified Cagniard method. First, however, it should be noted that upon replacing, in Eq.(2.9) by ,II according to Eq.(2.44), the inversion integral Eq.(2.7) becomes divergent, since ,II ~ 3/e as alAa,.. ~ 00. A physically meaningful interpretation of the integral must therefore be found. The asymptotic structure of the integrand Eq.(2.45) as a,..a,.. ~ 00 leads to the conjecture that the divergent integral contains a contribution that could be interpreted as a Dirac delta distribution operative at {:Z:I = 0,:Z:2 = O}. In view of this conjecture and of the contour deformation to follow, it is assumed that an interpretation in the sense of a Cauchy principal value around infinity will lead to a physically acceptable result. Adopting this interpretation, the integral on the right-hand side of Eq.(2. 7) is replaced by the limit ofthe integral over the rectangle {al E R, a2 E Rj - Llal < al < ~al' - Lla2 < a2 < Ll0<2} with Llal ~ 00 and Lla2 ~ 00. The corresponding evaluation will be carried out later in this chapter. The modified Cagniard method consists of applying the transformation given by Eqs.(2.18) and (2.19), in which q is kept real. The resulting integral with respect to p is extended over {p E T; -ill < p < ill} with Ll ~ 00. The integrand is continued analytically into the complex p-plane and an integration is carried out along the modified Cagniard path that follows from a continuous deformation of the imaginary p-axis and satisfies the equation (2.47)

38

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

where r is real and positive. The expression for '"'(II is written as (2.48) with
A2

8/3c2, 4/C2+q2.

(2.49) (2.50)

B2(q)

The non-real-axis part of the solution of the cubic equation (2.47) along which strictly increases, is the desired modified Cagniard path. In view of the condition that the joining circular arcs at infinity must, after compensation for the principal value around infinity, yield a vanishing contribution, only contours in the right half of the complex p-plane are considered. The desired contour consists of p = pII(r,q) in the first quadrant of the complex p-plane, together with its complex conjugate p = pII* (T, q) in the fourth quadrant of the complex p-plane. Its explicit form is obtained from Cardano's formula for the solution of the algebraic cubic equation; the relevant expression is given in Appendix B. For the moment we only need a number of properties that can directly be obtained from Eq.(2.47).
T

The first is that p = pII (r, q) and p = pll*( r, q) are finite arcs that leave the real p-axis at pF(q), at which point the value of r is denoted by TfI(q), and return to the real p-axis again at pf(q), at which point the value of r is denoted by TfI(q). Since both arcs leave the real p-axis at the same point and return again to the real p-axis at the same point, these points are double roots of the cubic equation (2.47). Consequently, or/op = 0 at these points. By differentiating Eq.(2.47) with respect to p, this condition can be expressed as (2.51) Inspection of the behavior of that pV(q) < B(q) and p~l(q)
t:

as a function of p along the real p-axis reveals

> B(q), which implies that (c£. Eq.(2.47))

(2.52)

SECOND-ORDER THIELE APPROXIMATION

39

c
6

lOOOm/s ;

""-..... -,

..

.....

__T[I(q) ...........
4

...

-.
..
\.
-5

lOOOm ;

X3

250m

<;

J/ »:
//
f

// ,

,/

// ...

~3

-T{l(q)

1\\

'.

'\

l
\

/ I

'. \

.../
o

-15

-10

10

15

~qc Figure 2.7: The of integration. and

Tf~ vs. ,

q-relationsj the area between the curves is the domain

(2.53)
The equality in the latter equation occurs in the exceptional case that r = 0, which will be discussed below. It follows that TfI(q) < TP(q) for all q, as shown in Figure 2.7. The second property is that along the arcs p = pH ( T, q) and p = pI! * ( T, q), strictly increases. This is proved by a reductio ad absurdum. Differentiation

40

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

of Eq.(2.51) with respect to p leads to

82T 8p2
from which it follows that

---

61x31
C

A2[B2(q) + 3p2] [B2(q) _ p2j3 ,

(2.54)

(2.55) and (2.56) Considering the local Taylor expansions about p = p{l(q) and p = pf(q), it then follows that along p = pII (T, q) and p = pII*( T, q), T increases away from pfI(q) and towards p~I(q). Now, if somewhere on the arc p = pII(T,q) there were other points at which 8T J 8p = 0, these points must be even in number with a minimum of two, while on the arc p = pII*( T, q) there would have to be the same number of such points. This would imply that the quartic equation (2.51) in p would have at least six solutions, which is impossible. Hence, along the modified Cagniard paths no zeros of 8T J 8p occur and T strictly increases. In the limit IX31 1 0 and r =I 0, Eq.(2.51) implies that both points of intersection, p{l(q) and p~I(q), approach the value B(q). The modified Cagniard contour then contracts to the point p( T, q) = B( q). Since, according to Eq.{2.47), in this limit T = pr is real, we have

Tf.~(q)

B(q)r,

(2.57)

which explains the equality sign in Eq.(2.53). Finally, differentiation of Eq.(2.47) with respect to q yields

8T 8q 8T!I(q)J8q>

31x31
C

2qA2 [B2(q)_p2J2'
and p = pf(q) 0 we have for all q

(2.58)

from which it follows that at p = pF(q) 0 and

8TiI(q)J8q>

> 0,

(2.59)

while q = 0 is a minimum of both T = TfI(q) and T = TP(q). To arrive at an equation for T = TII(q), we eliminate p from Eqs.(2.47) and (2.51). To this end, Eq.(2.48) is substituted into Eq.(2.47), which leads to (2.60)

SECOND-ORDER

THIELE

APPROXIMATION

41

where ( = 3Iz31/c. Differentiating the latter equation with respect to p and substituting the condition for double roots of Eq.(2.47), i.e., apT = 0, yields

2 T P - 2p~

B2(q) - -3= 0,

(2.61)

which equation is equivalent to Eq.(2.51). follows that

2 [ B2(q)r P - 2p 2(T _ ()

p2 _ B2(q) is obtained, which is substituted into Eq.(2.60).

From Eq.(2.61) an expression for In this way, it

2

+~

T-

(]

+B

(q)

+ 2(T

3(A2 _ () = O.

(2.62)

It is observed that Eqs.(2.61) and (2.62) are equivalent to Eqs.(2.60) and (2.61).
Using the fact that the roots of Eqs.(2.61) and (2.62) must coincide, a systematic elimination of the square-root expressions in the resulting equality leads to the final equation: 42 ["3B (q)

+ 2(TII(q)
_ ()

3(A2]

_ ()

- ()

B4(q)r 4 [ 6(TII(q)
= 0,

B2(q)(TII(q) 2r

(A2] [TII(q) - ( 2r 6r

+ 2(TIl(q)

B2(q)r] - ()
(2.63)

where we have set T = TIl (q). This equation is equivalent to Eq.(B.7) in Appendix B. (It is noted that the theory of resultants of two algebraic equations [55] can also be used to arrive at Eq.(2.63).) By setting q = 0, and substituting (, a quartic equation for the arrival times, i.e., the minima TB(O) of Tf.~(q), is obtained:

TIl(0)4 _1O(':31)TIl(0)3

-6 ('Zc31) [9(':3I
+27

C:

+ 4[9('z;l)

2_

2(~)2]TII(0)2

+ 4(~r]TIl(0)

3I

+ 36 ('ZC31

r(~r

+ 16(~f

= O.

(2.64)

Now that the properties of the modified Cagniard path {p = pIl (T, q) Up = pIl*( T, q)} have been investigated, the contour deformation leading from the integration along the imaginary p-axis to the one along the modified Cagniard path is further analyzed.

42

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

The body waves

In view of Cauchy's theorem the integral along the imaginary p-axis is equal to the sum of the integrals along the semicircle = {p E Cj Ipl = .6., Re(p) > O} and the modified Cagniard path (Figure 2.8). In addition, we may encounter the simple pole p = p~I of aU at the simple zero of "(II (Figure 2.9). From Eq.(2.48) it follows that

ct

(2.65)

(cf. Eq.(2.50)). The pole contributes to the Green's function when p~I ~ pfI j its contribution will be determined later on. (Since pV (q) > B( q), the
situation p~I > p~I can never arise.) On all is asymptotically equal to 2) (c/6s exp( -3slz31/c). Using the corresponding asymptotic integrand, and carrying out the transformation inverse to Eqs.(2.18) and (2.19), it follows that the contribution of results into

ct,

ct

(2.66) Furthermore, by taking the contributions from p = pll(T, q) and p = pII*(T,q) together (note that the integrand satisfies Schwarz's reflection principle), the contour part G{! of the expression for GIl is obtained:

where also the property has been used that the integrand is an even function of q. Interchanging the order of integration yields (see Figure 2.7)

SECOND-ORDER

THIELE

APPROXIMATION

43

4~--~,,----~;-----'r-----~----,
.................... ; ·T·········..··....····..·..· !
1

= lOOOm/s

;r

= lOOOm

X3

= lOOOm

L····..·......
i

···vertical singularity/
·····i··········..··..···..·.. ·····J

··.. .. · t··············. .

:g
~

!
.j

·",,····+..······.·········· ..·. ..
:..

2 .. ······. ····· · .. ··· ..i····::::::::·::::::::::tr:::::::: .

,.........

:::::I~" ; 1\
:.......... ..,

i
...........................

body wave

1...

'",I

\
"

,/

[\

\\.
"
!

-1~···......····················;·~·············.~·..........·..······························· .. · · · · ·· ·········i .. · · ·· · · , ; · .. · .. ..1

,~~
j

V···.....··· ···.. .J / /
+

-2 ~

j .............• ,

1;
, q=0
4

-3 ~

>

,./0'-:

·1

~
o

-4~--~~----~------~----+-----~
2 3
--t

-V

"

_ ~17."
Re{pc}

. 5

Figure 2.8: The modified Cagniard contour and the pole location for a "small" horizontal offset in the second-order Thiele approximation.

---

-------------

44

THIELE-APPROXIMATED

SCALAR

SPACE-TIME

WAVES

4~--~~------~!------r------r-----.

c = lOOOm/s; r

= lOOOm X3 = 250m ;

-..- - +- - i - j..! i

! i··-_··-·vertical singularity ~
!

·····t·········-

..-

-..-1.-- ! '
---t---- -

t.
.

-+-·-----········i·-········-······--·········

!
'

- -...1.......... .
......... -_ ! i_.__.._...·__···
! !

······················-···· ..r··-··-···············r·········-·--··-r··· .. ..·

................. -

..

·····;.·······-····--·--······-4·----··-.-

..;.-

-----

................. wave ,,_·+·-··-·-··_·····-··_.1··_···············~ bOd~ . j ..·· ··..············-·····-----·-··f

o
-1

· ..... tt -._. -........... ,

!
;

"head wave"
--·--·i···-········-·· ..

+_.-.-_ _.. +.. _-._ ..._._+-_ _ . _

;

;.·.·····_···_··---·.·.·-rj! ----1--·
,I

- -..--.- ..- ..

.. ··········· __ __ · .. ....._. .... .. _ ...

i- -_ _._-..

i:
!

··-·····--!·····~ ······· --.~.-.--.--i :

-2
.............. _ __

··I·····················_······
:

..·_············l -.. .... t······

I
I

!
!

-..~- --I

-..- -

--···_-·t·······-·······_ ..-·_··
_..-e+.•._..•__ ..•...__ .-.I

I!
I

, i ··········----···-······t··· .. ·-····--··-··········i·
!

_ .._

-3 .. ::~::::::·::::::::.:::::::t::~::::.:::::~:::~ ....

! I ! J -4+---~~-----+----~------~----4 3 4 o 2 -+ Re{pc}

·r~~~::~~-.~~::t:::::~:~:.::~.~J·· ..~-:~-··-·
5

!
!

Figure 2.9: The modified Cagniard contour and the pole location for a "large" horizontal offset in the second-order Thiele approximation.

SECOND-ORDER

THIELE

APPROXIMATION

45

IT=T{I(O)

Tfl(O)

exp( -sr)dr exp( -sr )dr

Q[I(T)

q=O

7r-2Im{[2-yll(pll, q)tl(8pII 7r-2Im{[2-yll (pI!, q)tl(

j8r)dq} 8pI! j ar )dq},

(2.68)

00

lQ[I(T) q=Q~I(T)

T=Tfl(O)

in which q = Q{l(r) and q = Q~I(r) are the (unique) inverses of r = TfI(q) and r = TfI(q) respectively. As Eq.(2.63) shows, these inverses follow from a cubic equation in q2. Application of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function, taken at positive real values of s, leads to the time-domain results

H[t - TfI(O)]
(2.69) and

ug (xm, t) = o, f,t=o f(t - t')Gg (xm, t - t') dt',

(2.70)

from which it is evident that t = T[I(O) and t = TfI(O) of the two wave fronts. Through the substitution

are the arrival times (2.71)

= QiI(t)sin('IjJ)

with

0::; 'IjJ::; 7rj2,

the inverse-square-root singularity at the upper limit of the first integral on the right-hand side of Eq.(2.69) is removed, while through the substitution

= {[Q~I(tWcOS2('IjJ)

+ [QiI(tWsin2('IjJ)}1/2

with

0::; 'IjJ::; 7rj2,

(2.72)

the inverse-square-root singularities at both limits of the second integral in Eq.(2.69) are removed. After these substitutions, the integrals are evaluated numerically with the aid of the trapezoidal rule. Figure 2.10 shows the result in the form of a snapshot in space for X3 ~ O. The wave with the arrival time T[I(O) will be denoted as the fast body wave; the wave with the arrival time TfI(O) will be denoted as the slow body wave.

46

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

c = 750m/s ; t = 0.4508
SSP
-+ horizontal offset (m)

source location 0.0 -448.0

448.0

500.0

Figure 2.10: Snapshot of the Cagniard-path contribution to the Green's function in space for :1:3 ~ 0 in the second-order Thiele approximation. (The horizontal offsets are numbered sequentially.)

SECOND-ORDER

THIELE

APPROXIMATION

47

c = 750m/s ; t = 0.450s
SSP
--+

horizontal offset (m) source location

-448.0

448.0

Figure 2.11: Snapshot of the pole contribution to the Green's function in space for 2:3 ~ 0 in the second-order Thiele approximation. (The horizontal offsets are numbered sequentially.)

48

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

The precursor
The contribution from the simple pole p&! to the Green's function is analyzed next. The location of the leftmost point of intersection p{l of the modified Cagniard path with the real axis satisfies the relation (d. Eq.(2.51)) (2.73) Upon substituting p&! (cf. Eq.(2.65)) into Eq.(2.73), the value q{/ for q at which p&! and p{1 coincide is found as (2.74) Note that ql! is real only if r ~ to q yields

~J3IZ31. Differentiating
q _ p2])

Eq.(2.73) with respect

8qp=
Since pf!(q)

+ ([4Iz31/rc3]/[B2(q)
8qpfI(q)~

forp=pf!(q)·

(2.75)

< B(q), Eq.(2.75) leads to

Pi

//()

(2.76)

when

s > O. Further,

from Eq.(2.65) it is clear that

q 8qp~! (q) = Il ().

Po q

(2.77)

Hence, when q = qlI, 8qp&!(q) > 8qp{1(q), so that the pole p&! lies to the left of prJ if q2 < (qJ!)2. In view of Eq.(2.74), this inequality cannot be satisfied when
r

< ~J3IZ31. Hence, the pole contributes only when r the residue, the relevant contribution is found to be
GIl _ ~
o97rc

> ~J3lz31. Evaluating

d q,

lqlI exp[-sr(q2+
q=O

(q2

4/3c2)1/2

+ 4/3c2)1/2]

(2.78)

which, upon replacing r(q2 Go A

+ 4/3c2)1/2 by T, takes the form

II _

2_ lTl.~
97rC

..-=Tl,{ (T2 - 4r2 /3c2)I / 2

exp( -ST)

d
T,

(2.79)

SECOND-ORDER

THIELE

APPROXIMATION

49

in which

T1,1 = (2/v'a) Of

(ric)

and

Tl,~ = (4r191x3D

(ric).

(2.80)

On account of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function, taken at positive real values of s, the time-domain counterpart G&I of then follows as (2.81) This contribution (see Figure 2.11), which is only present if r > IX31~v'3, is to be added to the contribution from the modified Cagniard path Eq.(2.69). It can be interpreted as a head-wave-like precursor to the fast and slow body waves.

The case r = 0
For an observation position straight above or below the source, we need the expression for the Green's function at r = o. In view of the occurrence of the Dirac distribution in Eq.(2.68), the limit as r 1 0 has to be taken with care, since the substitution of the value r = 0 in the Fourier representation for the Green's function leads to a divergent integral whose interpretation as a Dirac distribution is not obvious. To circumvent this difficulty, we leave the term in Eq.(2.68) containing the Dirac distribution intact and take in the subsequent terms in the expression the limit r 10. The corresponding contour then tends, in a non-uniform manner, to the imaginary p-axis (the points at infinity in the complex p-plane have to be identified with one another). When r = 0, the cubic equation for the modified Cagniard path reduces to a quadratic one (cf. Eq.(2.47))j the latter is solved by pII(r,q)=iB(q) for r r - TII(q)] [ TP~r
1/2

(2.82)

= 0,

TfI(q)

< r < TP, together with its complex conjugate, in which

(2.83)

50

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

and (2.84) For this case, the contour coincides with the imaginary axis and closes at infinity as opposed to the contours for r f:. 0 on which pIl remains bounded. The path of integration intersects the real axis at the origin where T = TfI(q) and tends to infinity as T I T.JI. The relation inverse to Eq.(2.83) is q = QF (T), with (2.85) for IX31/c < T < 3Ix31/c, while Eq.(2.85) implies that only the first integral of Eq.(2.69) contributes. This integral can be evaluated, and the Green's function reduces to

(2.86)

The case

x3

In the case X3 = 0, the modified Cagniard contour coincides with the positive real p-axis. In view of the analyticity of the integrand away from the pole we have, when X3 = 0, a contribution from the arc at infinity and a contribution from the pole. The latter is given by (cf. Eq.(2.78))

all = 2_ 1
o 97rC

00

q=O

exp[-sr(q2 + 4/3c2)1/2] dq. (q2 + 4/3c2)1/2

(2.87)

The space-time counterpart is found to be (cf. Eq.(2.81))

G~l
in which Eq.(2.69))

97rc[t2

~T6.D2P/2H(t

- Tt,{),

(2.88)

Tl,{ is given in Eq.(2.80).

The total Green's function is then (cf. (2.89)

ARRIVAL

TIMES

51

Figure 2.12 shows a snapshot in space of the total Green's function C{f for 2:3 ~ O.

+ C~[

Partial differential equation for GIl

By differentiating the right-hand side of Eq.(2.46) twice with respect to 2:3, and using that under the transformations of Eqs.(2.5) and (2.6) we have 8t ~ s and 8" ~ -isa", we can reconstruct the partial differential equation that ClI satisfies, viz.,

{c-28;(83 + c-18t) {c- 8;(83

2 2

[(1/4)83 + {3/4c)8t]8"o,,} [(1/4)83

-

= -[c- 8;

c- 8t)
1

(3/4c)8t]8v8v}ClI
(2.90)

(1/4) 8,xO,x]26(2:1' 2:2, 2:3)H(t).

The first partial differential operator in braces on the left-hand side of Eq.(2.90) is the second-order approximation to the wave operator for a wave propagating in the direction of increasing 2:3; the second partial differential operator in braces on the left-hand side of Eq.(2.90) is the second-order approximation to the wave operator for a wave propagating in the direction of decreasing 2:3 [33].

2.5.

Arrival times as functions of horizontal offset

A Thiele approximation of the vertical slowness is accurate for small values of 0.1 and 0.2. In this sense, it aims to approximate (cf. Eq.{2.7)) wave propagation in the vertical direction. How the results of this approximation compare with the ones of the asymptotic ray solution along a vertical ray is the topic of this and the next sections. In the present section, relations between horizontal offset and arrival time are derived for the different approximations. We keep 12:31 fixed and let r vary. The vertical travel, or intercept, time is denoted by Tmin = 12:31/c, and the arrival time by T.

52

THIELE-APPROXIMATED

SCALAR

SPACE-TIME

WAVES

c = 750m/s ; t = 0.450s
SSP
--+

horizontal offset (m) source location 448.0

0.0

-448.0

500.0

Figure 2.12: Snapshot of the total Green's function in space for :1:3 ~ 0 in the second-order Thiele approximation. (The singular part at r = 0 is not shown. The horizontal offsets are numbered sequentially.)

ARRIVAL

TIMES

53

For the exact wave motion only a single body wave is present and its arrival time satisfies the equation (d. Eq.(2.4)) (2.91) which leads to (2.92) This is also the relation predicted by ray theory (where straight rays emanate from the source to the point of observation). For the first-order Thiele approximation there is a single bodywave whose arrival time satisfies the equation (cf. Eq.(2.33)) (2.93) the solution of which is (2.94) In addition, there is a head-wave-like precursor whose arrival time is given by (cf. Eq.(2.38» (2.95) T6,1(T) = cT/V2 for T ~ 2Tmin; this precursor is only present if r 2.13).

> R6, where R6

CTmin v'2 (see Figure

For the second-order Thiele approximation we have Eq.(2.64),

T4 - 10TminT3 -6Tmin[9T.!in

+ 4 [9T.!in + 4(~rJT

2(D 2JT2
+ 27T,!in + 36T.!in (~r + 16(~r = 0,
(2.96)

the two real-valued positive solutions of which are (2.97)

54

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

lOOOm/s

X3

lOOOm

,---....

54000
-._;> <l)

Cd

~ o 3000 ~ o o
H

.~ 2000
,..q

1000

2
-t

time (8)

Figure 2.13: Travel time vs. horizontal offset (diffraction) curves associated with the arrivals of the exact and first-order Thiele-approximated Green's functions.

ARRIVAL

TIMES

55

c = lOOOm/s ;

X3

= lOOOm

5000r---------------------------------------------~~--_.

~ ~

(l)

~ ~

o 3000

N .i:: 2000

...c:

1000

slow wave

~ time (s)
Figure 2.14: Travel time vs. horizontal offset (diffraction) curves associated with the arrivals of the exact and second-order Thiele-approximated Green's functions.

56

THIELE-APPROXIMATED

SCALAR

SPACE-TIME

WAVES

for the fast body wave, and (2.98) for the slow body wave, where (2.99) and (2.100) At T = 3Tmin there is a spatial 5-function singularity on the vertical axis, i.e., right at the cusps of the wave front of the slow wave. In addition, there is a head-wave-like precursor whose arrival time is given by (d. Eq.(2.80)) (2.101) which is present only if r

> R{/ with R~I

CTmin ~J3(see Figure 2.14).

It is noted that the vertical travel times associated with the (fast) body wave in the first- and second-order Thiele approximations are exact. Furthermore, it follows that R~I > R~.

2.6.

The local behavior of the waves near the wave fronts in the vertical direction

In this section we compare the local behaviors of the waves near the wave fronts in the vertical direction for the exact Green's function and its first- and second-order Thiele approximations.
The exact Green's function at {O,0, Z3} reduces to (d. Eq.(2.4)) (2.102)

LOCAL

BEHAVIOR

NEAR THE WAVE FRONTS

57

The first-order

Thiele approximation at {O,0, X3} yields (cf. Eq.(2.35))

G (0,0, X3, t)

-4 -H(t
x d:

- IX31/c).

(2.103)

Replacing t in the denominator of Eq.(2.103) by t = IX31/c + (t -IX31/c), and using a Taylor expansion of the resulting expression about the arrival time t = IX31/c, yields

G1(0, 0, X3) t) = [-4 11 I 7r X3

~(t

47rX3

- IX31/c) + O[(t - IX31/c)2J] H(t -IX31/c).

(2.104) The right-hand side is the ray series expansion along a vertical ray. (The ray series expansion of the first-order Thiele approximation in an arbitrary direction is discussed by Holden and Gorman [56J.) The leading term in Eq.(2.104) coincides with the exact result (cf. Eq.(2.102)), whereas the next-order term already deviates from it. This result could also have been predicted from the parabolic equation in ray-centered coordinates occurring in the paraxial ray theory (see, e.g., [57]). In the second-order Thiele approximation at {O, 0, X3} the Green's function reduces to (cf. Eq.(2.86)) lim 1

(XI,X2)--->(O,O)

[G{j(Xl,X2,X3,t)
x~ ci

- :'(t - 3Ix31/c)H(t
6

- 3IX31/C)cS'(X1!X2)] (2.105)

= -( 7rct

II 3 X3 -

)2 [H(t -IX31/c)

- H(t - 3Ix31/c)J.

Replacing t in the denominator of Eq.(2.105) by t = IX31/c + (t - IX31/c), and using a Taylor expansion of the resulting expression about the arrival time t = IX31/c of the fast body wave, yields

II

(0,0, X3,t)

= 47rlx31

[1

3c:.! + 167rIX313 (t-IX31/c)

+0 [ (t-lx31/c)

3]

H(t-lx31/c),
(2.106)

which shows that now the first two orders are exact (cf. Eq.(2.102)). The arrival-time curves and the snapshots (see Figures 2.1, 2.6 and 2.12) show that the accuracy of asymptotic ray theory is achieved along the vertical

58

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

ray only; in oblique directions of propagation there is a discrepancy. In oblique directions a better approximation is found if the vertical slowness is expanded about the actual direction of the ray trajectory connecting the source and receiver. Such an approach requires an adapted form of Thiele's expansion about an oblique direction of propagation, which is discussed in Appendix A.

2.7.

Concluding

remarks

The foregoing analysis has revealed the general consequences that spectraldomain Thiele approximations for the vertical wave slowness have for the corresponding space-time domain wave motion. The most striking feature is that Thiele approximations lead to a number of non-physical artefacts that do not occur in the exact wave motion. First of all, for even-order approximants, a singularity in the wave motion is observed right above and below the source, which is absent in both the exact wave motion and in odd-order approximants. Further, for horizontal offsets that are large compared to the vertical offset, non-physical head-wave-like precursors occur. In addition, any Thiele approximant of order higher than one introduces "slow" waves, which again are absent in the exact wave motion (see also Collins [58]). In the odd-order approximations, however, they never occur along a direction of propagation which lies interior to a certain cone centered on the vertical axis. On the other hand, Thiele approximations reproduce amplitudes along the "vertical ray" as accurately as the asymptotic ray theory does. However, in contrast to asymptotic ray theory, the Thiele-approximated wave theory automatically initializes the ray amplitudes in the immediate neighborhood of the source. In conclusion, one can state that Thiele approximations are sufficiently accurate inside a cone for not too large horizontal-versus-vertical offset ratios and in a restricted time window, i.e., in the neighborhood of the arrival of the physical body wave. Further, in view of the associated non-physical singu-

CONCLUDING

REMARKS

59

c
SSP

= 750m
-->

8t

= 0.4508

horizontal offset (rn) source location

Figure 2.15: Snapshot of the total Green's function in space for :1:3 ~ 0 in the third-order Thiele approximation. (The horizontal offsets are numbered sequentially. )

60

THIELE-APPROXIMATED

SCALAR SPACE-TIME

WAVES

larities straight above and below the source, the even-order approximants are to be avoided. As an example of this, we show a snapshot of the third-order Thiele-approximated Green's function in Figure 2.15 (taken from Serafini and De Hoop [24]).

61

Chapter 3 Propagation of acoustic waves in fluid media in the first- and second-order Thiele approximations
The purpose of this chapter is to investigate for different source-receiver pairs the accuracy and the artifacts of the solutions to the acoustic wave equations in a fluid medium associated with the first- and second-order Thiele approximated slowness surface. Aspects of two canonical configurations are analyzed, viz., the structure of the body waves in an unbounded homogeneous fluid, and the wave reflection against the boundary with a medium with a higher wave speed. By evaluating, with the aid of the modified Cagniard method, the space-time Green's functions in the relevant two-medium configuration it is found that, instead of generating the physical head wave, the Thiele approximation generates artificial head-wave-like arrivals. In the analysis, the reflection and transmission coefficients are approximated in accordance with the Thiele expansions of the vertical slownesses. Numerical results are presented.

62

THIELE-APPROXIMATED

WAVES IN FLUIDS

3.1.

Introd uction

This chapter builds on the previous study [59] on the Thiele approximations of the scalar wave equation to derive a full Thiele-approximated theory of acoustic waves in fluids. In particular, we investigate for different sourcereveiver pairs the accuracy and the artifacts of the solutions to the acoustic wave equations in a fluid medium associated with the first- and second-order Thiele approximated vertical slownesses. Since our previous study revealed that the Thiele approximation generates non-physical artificial "head waves", we presently investigate the approximation's capability to describe the physical head waves that occur in the wave reflection against the boundary with a medium with a higher wave speed. By evaluating the space-time Green's functions in the relevant two-medium configuration it is found that, rather than generating the physical head wave, the Thiele approximation creates additional artificial head-wave-like arrivals. In the analysis, the reflection and transmission coefficients are approximated in accordance with the Thiele expansions of the vertical slownesses. Similar approximations have been introduced by Richards and Frasier [60] and Shuey [61], in combination with the linearization of the reflection coefficients in the medium parameter contrasts across the interface. The first-order Thiele, or parabolic, approximation has already been extended in different ways by Corones [62] and McCoy et al. [63] to media with interfaces across which the medium parameters jump by finite amounts. Further, the modification of Snell's law entailed by the first-order Thiele approximation for a two-medium configuration (represented by a transition zone) has been analyzed by Hatton et al. [64]. First, the wave-matrix formalism for acoustic waves in an arbitrarily vertically inhomogeneous fluid is developed, after which specifically the reflection against and transmission across a planar boundary between two different semi-infinite homogeneous, isotropic fluids are analyzed in detail. Both the exact solution to the problem and the solutions arising from the application of the first- and second-order Thiele approximations to the corresponding vertical slownesses are obtained with the aid of the modified Cagniard method

WAVE-MATRIX

FORMALISM

63

[23,49,50,51].

3.2.

The wave-matrix formalism

Consider linearized acoustic waves in an isotropic fluid occupying some subdomain of three-dimensional space 'R}. The point at which the wave motion is observed is specified by the coordinates {XI, X2, X3} with respect to a Cartesian reference frame with the origin 0 and three mutually perpendicular base vectors {i1' i2, i3} each of unit length. In the indicated order, the base vectors form a right-handed system. In accordance with geophysical convention, X3 is taken as the vertical, or depth, coordinate (which increases in the downward direction), leaving XI, X2 as the horizontal coordinates. The subscript notation for Cartesian vectors and tensors is employed. Lowercase Latin subscripts are used for this purpose, and the summation convention applies; the subscripts range through the values {I, 2, 3}. The position vector will be denoted by a: = xmim. The time coordinate is denoted by t. Differentiation with respect to Xm is denoted by am j at is a reserved symbol for differentiation with respect to t. Given these notational preliminaries, the acoustic wave motion is characterized by the system of first-order partial differential equations akP where p Vr p
K,

+ patVk K, atp + arvr

(3.1)
q,

(3.2)

acoustic pressure [Pal, particle velocity [m/s], volume density of mass [kg/m3], compressibility [Pa-1], = volume source density of force [N/m3]. volume source density of injection rate [s-1],

Jk

64

THIELE-APPROXIMATED

WAVES IN FLUIDS

It is assumed that

and It are both time independent and translationally invariant in the :1:1- and :l:2-directions; they depend on :1:3 only. First, the case where the changes with :1:3 are continuous is considered; the case with a discontinuous change in medium properties will be handled separately. Both continuous and discrete reflection and transmission require the decomposition of waves. The homogeneous fluid is a special case of the continuum. However, since the decomposition is no more difficult for the continuum than for the one for the homogeneous medium, we discuss in the general analysis the former.
p

Causality of the wave motion is enforced. Assuming that the sources generating the wave field are switched on at the instant t = 0, causality implies that the wave field quantities satisfy the initial conditions

p{:l:m,t) Vr{:l:m,t)

0 0

for for

t<O

and all and all

(3.3) (3.4)

t<O

In view of the time invariance of the medium, the causality of the wave motion can be taken into account by carrying out a one-sided Laplace transformation with respect to time, taking the time Laplace-transform parameter 8, which is in general complex-valued, to be in the right half {Re{s} > 0) of the complex s-plane and requiring that the transform-domain wave quantities are bounded functions of position in all space. To show the notation, we give the expression for the Laplace transform of the acoustic pressure

(3.5)
Employing the translational invariance of the configuration in the horizontal directions, also a Fourier transformation is applied in the horizontal plane with real transform parameters 8Ql and SQ2. Lowercase Greek subscripts - ranging over the values {I,2} - are used to indicate the horizontal components of a vector in three dimensions; to them, too, the summation convention applies. For the acoustic pressure, the relevant transformation is

WAVE-MATRIX

FORMALISM

65

The variables ia,.. are the horizontal components of the slowness vector of the wave motion. The transformation inverse to Eq.(3.6) is given by

The Fourier components exp( -isa6:t6) are the horizontal eigenmodes of the system of differential equations (3.1) and (3.2); they propagate independently. Under the Laplace transformation, assuming zero initial conditions, we have at ~ s. The Fourier-domain counterpart of oiJ is given by -isaw With these rules, the first-order acoustic wave equations (3.1) and (3.2) transform into
03P -isaK.p
SK

+ SPV3
+
sp VK.

(3.8) (3.9) q. (3.10)

P-

isapvp

+ 03V3

Equation (3.9) is rewritten as (3.11) which result is used to eliminate the horizontal components of the particle velocity from Eq.(3.10). The derivatives of the remaining components of the wave field in space in the :t3-direction can be expressed in terms of the Fourierdomain counterparts of the derivatives of the wave field with respect to the horizontal coordinates. The remaining equations are arranged as a single matrix differential equation. In this procedure, a subscript notation for the matrices involved will be employed. Uppercase Latin subscripts are used for this purpose; they range through the values {1,2}. Following De Hoop [65], we write 83Pr

+ sAI,JF]

= NI,

(3.12)

in which the elements of the acoustic field matrix

PI

are given by (3.13) (3.14)

FI =p,

= V3,

66

THIELE-APPROXIMATED

WAVES IN FLUIDS

the elements of the acoustic system's matrix A[,J by

AI,1
A
2,}

A2,2
-1

= 0,

(3.15) (3.16)

A},2 = p,

=P
=

Ot.pOt.p

+ "',
fh
by

(3.17)

and the elements of the notional source matrix

NI N2
Note that

j3,

= iOt./tp-I

t. + ij.

(3.18) (3.19)

A is independent

of s,

Through an appropriate linear transformation to be carried out on the acoustic field matrix, a wave-matrix formalism will be arrived at from which, at each depth level, a decomposition into up- and downgoing waves (with respect to the :l:3-direction) will be made manifest. To this end, we write (3.20) to transform Eq.(3.12) into

.L[,J (83WJ

+ SAJ,MWM)
A[,J.LJ,M

-(83.L[,J)

WJ

+ NI,

(3.21)

such that AJ,M, satisfying the equation

= il,JAJ,M,

(3.22)

is a diagonal matrix. Obviously, AJ,M is the diagonal matrix of the eigenvalues of A[,], and .L[,J consists of the corresponding eigencolumns. We denote W[ as the wave matrix and .L[,J as the composition matrix. To determine the solutions of Eq.(3.22), we introduce the column matrices i}±), according to
L[
-(+)

- (-) L[

= L[,I ,

= L[,2
f

(3.23) (3.24)

Let the diagonal elements of AJ,M be written as - (+) A-1,1 - ......

,

(3.25) (3.26)

A2,2

= "/(-),

WAVE-MATRIX

FORMALISM

67

then Eq.(3.22) decomposes into the two systems of equations (3.27) (3.28) The 1'(±) are known as the vertical slownesses. Note that the elements .i~±) compose the acoustic pressure and that the elements .i~±) compose the vertical particle velocity. The eigenvalue problem represented by Eqs.(3.27) and (3.28) (or Eq.(3.22)) implies the characteristic equation

det(A - A) = 0,
I.e.,

(3.29)

(3.30) The expression A2,1Al,2 can in turn be considered as the eigenvalue of the elliptic partial differential operator ..4.2,1..4.1,2, given by (d. Eqs.(3.16)-(3.17)) (3.31) corresponding to the eigenmode exp( -isQ"z,,). The square roots of b(±)]2 are taken such that Re{-y(+)} ~ 0 and Re{-yH} ::; o. We then have (d. Eq.(3.30)) (3.32) (3.33) where
l'

=(

QpQp

-2)1/2

(3.34) (3.35)

with Re{ 1'}

> 0 and

which is the local acoustic wave speed. Now that the eigenvalues of the system's matrix have been calculated, the columns of the composition matrix can be constructed, the normalization still being a degree of freedom. Three different choices will be discussed, each having its own specific advantages. First, we discuss the acoustic-power-flux normalization. This choice is the most appropriate one for mathematical considerations of, for example, the proof of

--------

-------------------------------------------------------

68

THIELE-APPROXIMATED

WAVES IN FLUIDS

convergence [65] of the generalized ray series. Second, the amplitude of the acoustic pressure in each eigencolumn is normalized to unity. This choice is the most appropriate one if the acoustic pressure is being measured. Third, the amplitude of the vertical particle velocity in each eigencolumn is normalized to unity. This choice is the most appropriate one if the vertical particle velocity is being measured.

Normalization to the vertical acoustic power flux

Imposing vertical transform-domain acoustic-power-flux normalization to the eigen-columns [66J , we have the condition

L-(-) 1

.i~+) .i~+»)
L-(-)
2

(0 1) (.i~+) .i(+)
10
2

L-(-)
2

.i~-») _
-

(1 0)
0 -1

(3.36)

Further, since Eq.(3.22) leads to

pL~±) = -y(±) Ll±),

, ,

we obtain (3.37) (3.38)

.i~+) = (2Yt1/2 .i~-) = (2y)-1/2

where

.i~+) = (Y/2)1/2, .i~-) = _(Y/2)1/2,

(3.39)
is the (transform-domain) vertical acoustic wave admittance. Hence,

_
L= The decomposition matrix L

(2Y)-1/2 (Y/2)1/2

(2Y)-1/2) _(Y/2)1/2

(3.40)

L:'
-

follows from Eq.(3.40) as

- -1 _

(Y/2)1/2 (Y/2)1/2

(2Yt1/2) _(2y)-1/2

(3.41) is then found to (3.42)

From Eqs.(3.40) and (3.41), the interaction matrix be

-.i-183.i

WAVE-MATRIX

FORMALISM

69

in which (3.43) is the local (transform-domain) reflection coefficient density. Equation (3.42) shows that the interaction matrix is, for the normalization under consideration, off-diagonal. Setting

Xl = X2 =

(1,-1 h.MNM = (Y/2)1/2 (1,-1 h.MNM

Nl + (2y)-1/2 N2,

= (Y/2)1/2

s. - (2Yt1/2

(3.44) (3.45)

N2,

Eq.(3.21) finally transforms into the coupled system of "one-way" wave equations for the normalized wave amplitudes 83W1

+ 8,Wl = ;W2 + Xl,

-

(3.46) (3.47)

83W2

8,W2

= ;W1

+ X2,

from which the interpretation of ; as a reflection coeffient density is evident. In a sub domain of space in which the medium is homogeneous, we have r = O. In this case, the interaction terms on the right-hand sides of Eqs.(3.46) and (3.47) vanish and the upgoing waves (characterized by W2) and downgoing waves (characterized by WI) travel independently.

Normalization to the acoustic pressure

The solution of Eqs.(3.27) and (3.28), normalized so that the acoustic pressures are set equal to unity, yields the composition matrix L

- = (1 1)
Y _Y

(3.48)

Note that under this condition p = matrix represent acoustic pressures. sponding to Eq.(3.48) is given by

WI + W2,

i.e., the elements of the wave The decomposition matrix L-1 corre-

(3.49)

70

THIELE-APPROXIMATED

WAVES IN FLUIDS

hence the interaction matrix

-£-183£

becomes (3.50)

with (3.51) Note that the interaction matrix can be related to the scattering matrix in a discontinuous medium (see Appendix C, Eqs.(C.16)-(C.19)) upon identifying its diagonal elements with a transmission coefficient density i according to -r = i-I - 1. Furthermore, the notional sources transform into (d. Eqs.(3.44)-(3.45))

Xl X2

= (2Y

(Y Nl

+ N2),

(3.52) (3.53)

= -(2Yt

(-YNI

+ N2).

The coupled system of "one-way" wave equations for the acoustic-pressure wave amplitudes now becomes

83Wt 83W2

+ s/WI
-

= (2Y

tI(8 Y)
3

(-WI

+ W2) + (2y)-1
- W2)
-

(YNt

+ N2)'
(3.54)

s/W2 = (2y)-I(83Y)(Wl

(2y)-1 (-YNl

+ N2).
(3.55)

In the Thiele approximations to be discussed later on, Y becomes a rational function of the horizontal slownesses and Eqs.(3.54 )-(3.55) are transformed into space-time partial differential equations upon replacing 8 ~ 8t and -isa
jJ ~

8jJ"

Normalization to the vertical particle velocity

The solution of Eqs.(3.27) and (3.28), normalized so that the vertical particle velocity is set equal to unity, yields the composition matrix L-

__ (Y-11 _ y-1 ) 1

(3.56)

WAVE-MATRIX

FORMALISM

71

Note that under this condition V3 = WI + W2, i.e., the elements of the wave matrix represent vertical particle velocities. The decomposition matrix L-1 corresponding to Eq.{3.56) is given by

L-1 = (Y/2)
hence the interaction matrix

(~1 ~=:)
j

(3.57)

-L-183L

becomes (3.58)

with (3.59) Furthermore, the notional sources transform into (cf. Eqs.{3.44)-{3.45))

Xl

= =

(1/2) (Y N1 + -(1/2) (Y Nl

X2

N2), - N2).

(3.60) (3.61)

The coupled system of "one-way" wave equations for the vertical-particlevelocity wave amplitudes now becomes

83W1 83W2

+ 8')'W1 = _(2y)-1{8 Y)(-W1

3 -

+ W2) + (1/2)(YN1 + N2),

(3.62)

8')'W2 = _(2y)-1{83Y)

(WI -

W2)

(1/2) (Y N1 - N2).
(3.63)

In the Thiele approximations to be discussed later on, Y becomes a rational function in the horizontal slownesses and Eqs.{3.62)-{3.63) are transformed into space-time partial differential equations upon replacing 8 ~ 8t and -i8aiJ ~

8w
With this, the transform-domain formulation of the acoustic wave problem in an arbitrary, continuously layered isotropic fluid, and its reformulation as a coupled-wave problem of up- and downgoing constituents, have been completed. At interfaces across which the constitutive parameters jump by finite

72

THIELE-APPROXIMATED

WAVES IN FLUIDS

amounts, the acoustic wave equations have to be supplemented by boundary conditions which require the continuity of the acoustic pressure and the vertical particle velocity (see Appendix C). Apart from their adaptations to the physics, the different normalizations also have consequences for the corresponding first- and second-order Thiele approximations. In both cases the vertical wave slownesses and, consequently, the vertical wave admittance, which contain square-root expressions in the horizontal slownesses, are approximated by rational functions. Now, in the acousticpressure and vertical-particle-velocity normalizations of the eigencolumns, the elements of the decomposition and composition matrices become rational expressions as well. In the acoustic-power-flux normalized eigencolumns, however, a quartic root occurs which leads to square-root expressions for the elements of the decomposition and composition matrices after the introduction of the Thiele expansion. To avoid having to deal with another Thiele expansion of these, purely as a result of the normalization chosen, we employ in our analysis only the acoustic-pressure normalization; the results for the vertical-particle-velocity normalization follow along a similar line of analysis. The space-time counterparts of the decomposition and composition matrices can for the latter two normalizations subsequently be identified with partial differential operators. This implies that the initialization of the wave matrix at the source position and the transformation back to observables at the receiver position can be represented by partial differential equations.

3.3.

The transform-domain solution for an unbounded homogeneous fluid

In this section we investigate the propagation of acoustic waves in an unbounded homogeneous fluid. As far as the excitation of waves is concerned, we investigate two cases of practical interest: the explosive point source (model for an airgun in marine seismics and for a monopole transducer in ultrasonics)

TRANSFORM-DOMAIN

SOLUTION

73

and the point force source (model for a vibrator in land seismics, where the equivalent fluid model for seismic waves in the earth is often employed, and for a dipole transducer in ultrasonics). As far as the observables are concerned, we consider the acoustic pressure (measured with a hydrophone) and the vertical particle velocity (measured with a geophone). In view of the spatial invariance of the medium, it can be assumed that the source is operative at the origin without loss of generality. Further, the composition matrix is taken to be normalized with respect to the acoustic pressure. For an explosive point source, we have (3.64) in which et( s) denotes the s-domain source signature. For a point force source, we have (3.65) in which j2(s) denotes the (vectorial) s-domain source signature. Eqs.(3.52)-(3.53) show, XI is of the general form Xl =
Wi

Then, as

5(:1:3),

(3.66) (3.67)

X2 = -W25(:l:3),
in which

(2y)-I et( s), (2Ytltt(s),

in the case of an explosive source, and

(3.68) (3.69)

WI = (2y)-l[p-Iia,J~(s) W2

+ yn(s)],
- yjg(s)],

(3.70)

(2YrI[p-lial<j~(s)

(3.71)

in the case of a force source. It follows from Eqs.(3.54) and (3.55), with T = 0, interpreted in distributional sense, that

Wi
W2

= wIH(:l:3) exp( -S,:l:3), = W2H( -:1:3) exp(s,Z3),

(3.72) (3.73)

74

THIELE-APPROXIMATED

WAVES IN FLUIDS

where H denotes the Heaviside unit step function. Employing Eqs.(3.20), (3.23) and (3.24), we then obtain an expression for the observables: (3.74) The acoustic-pressure and vertical-particle-velocity responses due to both the vertical force and injection space-time point sources can be grouped together in a Green's tensor G, that satisfies (3.75) Substituting Eqs.(3.68)-(3.71) into Eq.(3.74), and using the equalities (3.76) and (3.77) leads to (3.78) where

represents the direct source contribution,

1 -p S -183 )
_1

(3.79)

and

G = (2s2y)-1 G,
-(sp)-183G1 -(8p)-1[83G2

exp( -81Iz31).

(3.80)

We analyze the independent elements of

G through
(3.81)

= (2s t1sgn(z3)exp(-s1Iz31),
2 -

(3.82) (3.83)

s-26(Z3)]

= (Yj2s2)exp(

-811:2:31).

TRANSFORM-DOMAIN

SOLUTION

75

Then, we arrive at the following expressions for the explosive source:

pq

= S2(t(S)01'

(3.84) (3.85)

ii~ = S2qO(S)02' and the following expressions for the force source:

PJ
-J

-p-

'0 sf",(s)(-isa",)G1 .

'0 +s 2f3(S)G2-,

(3.86) (3.87)

V3 = -p

-1'0

sf",(s)(-~sa",)G2

+ s fJ(s)GJ.
2 '0

Using Eq.(3.11), the expression for ii!,J can be obtained from Eqs.(3.85) and (3.87). Through the necessary spatial and temporal differentiations and temporal convolutions, the space-time expressions for the acoustic pressure and the particle velocity can be obtained from the space-time Green's functions. The Green's function given by Eq.(3.81) has been discussed in detail in Chapter 2, so here we shall concentrate on the other two. Application of the inverse Fourier transformation Eq.(3.7) to Eqs.(3.82) and (3.83) yields 02(:V~":VJ's) = (8'/l'2t1sgn(:VJ) and OJ(:v,..,:V3,S)= (8'/l'2pt1
JOI'E'R.2

exp(-isa,..:v,..)exp(-s-yI:V31)dalda2 (3.88)

JOI'E'R2

-yexp(-isa,..:v,..)exp(-s-yI:VJl)da1da2'

(3.89)

It is observed that the inversion integrals in Eqs.(3.88) and (3.89) become

oscillatory integrals as :V3 0, which have to be interpreted in distributional -+ sense. The last integral is characteristic of the Dirac distribution singularity in Eq.(3.83). Strictly speaking, the contour deformation (in the exact and the first-order approximated OJ) applies to the difference

(1/ ps2)(-y /2) exp( -S-yI:V31) - (1/ ps3)5(:V3)

(in a weak sense), since lim (TJ/2)exp(-sTJI:V31) 1'11-+00 when Re{TJ}

= (1/s)5(:VJ)

> O. Note that Re{-y} > O.

76

THIELE-APPROXIMATED

WAVES IN FLUIDS

3.4.

The exact solutions in the space-time domain

The exact space-time-domain expressions for the Green's functions introduced in the previous section follow from Eqs.(3.7) and (3.5) upon using the modified Cagniard method. They have been discussed by several authors, e.g., [67]. In our notation they are given by G1 = 47rR H(t - R/c), G2 = -R3
47r t P

(3.90) (3.91)

sgn(z3)lz31 H(t - R/c)

(see Figure 3.1) and G3

87rpc

12R3[r2+c2t2(2-3r2/R2)]H(t-R/c)+21

5(Zl,Z2,Z3)t2

H(t) (3.92)

(see Figure 3.2, apart from the singular term), in which r

z,..Z,. )1/2 ,

(3.93)

are respectively the horizontal distance and the radial distance from the source to the receiver, and c = (pKt1/2 is the acoustic wave speed. The last term in Eq.(3.92) is interpreted as the direct source contribution.

3.5.

The first-order Thiele-approximated solutions in the space-time domain

I

The first-order Thiele approximation involves replacing 'Y by 71 and hence Y by yl, where (d. Eq.(2.15» 7 yI
1

c-1

+ (c/2)a,..a,.,
+ (c2/2)a,..a,.].

(3.94) (3.95)

(pctl[1

FIRST-ORDER

THIELE

APPROXIMATIONS

77

c
SSP

= 7 SOmis;

= 0.450s

_. horizontal offset (m) source location 0.0

500.0

Figure 3.1: Snapshot of the exact Green's function G2 in space for (The horizontal offsets are numbered sequentially.)

:1:3

> O.

78

THIELE-APPROXIMATED

WAVES IN FLUIDS

c
SSP

= 7 SOmis;
-+ horizontal

= 0.450s

offset (m) source location

0.0

-448.0

448.0

] ...
>

'" ~ o

....

.,

500.0

Figure 3.2: Snapshot of the exact Green's function G3 in space for (The horizontal offsets are numbered sequentially.)

:1:3

>

O.

FIRST-ORDER

THIELE

APPROXIMATIONS

79

To construct the corresponding space-time equivalents of Gf, G~ and G~, we apply the modified Cagniard method. This consists of transforming the variables of integration {all a2} with {a1 E 'R, a2 E 'R} into {p, q} with {p E I, q E 'R} (here, I denotes the imaginary axis) through
a1

= -ipcos(O) - qsin(O),

(3.96) (3.97)

a2

-ipsin(O)

+ qcos(O),
= r sine0),

where {r, O} is related to {:l:1' :l:2} by

:1:1

= r cos(0),

:1:2

(3.98)

with 0 ~ r < 00, 0 ~ 0 < 211', keeping q real, continuing the integrand analytically into the complex p-plane and carrying out the integration along the modified Cagniard path that follows from a continuous deformation of the imaginary p-axis and satisfies the equation (3.99) where T is real and positive. In view of the condition that the joining paths at infinity must yield a vanishing contribution, only contours in the right half of the complex p-plane are considered. Solving p from Eq.(3.99), it is found that the modified Cagniard path consists of p = pI (T, q) in the first quadrant of the complex p-plane, together with its complex conjugate p = pI*( T, q) in the fourth quadrant of the complex p-plane, where

in which (3.101) with (3.102) Equation (3.100) represents a straight line parallel to the imaginary p-axis that intersects the real p-axis at p = r/cJ:l:3J. Along this path, T strictly increases when going from the point of intersection with the real p-axis to infinity.

80

THIELE-APPROXIMATED

WAVES IN FLUIDS

Then, the integral along the imaginary p-axis is equal to the integral along the modified Cagniard path. In the resulting integral, r is introduced as the variable of integration. The corresponding one-dimensional Jacobian follows from Eq.(3.100) as (3.103) For of, given by Eq.(3.81), the procedure has been carried out in a Chapter 2, where it is found that

G{

47rc(t2_ ~Tl,l)2F/2
1

sgn(t - T1,2)H[t

- TI(O)] (3.104)
(3.105)

+ 27rC[t2 (Tl,l)2j1/2X[Tt,1.Tt,21(t), _
where in which

(3.106) Let us next consider the Green's function given by Eq.(3.82). Taking the contributions from p = pI and p = pI* together (note that the integrand satisfies Schwarz's reflection principle), Eq.(3.88) leads to

oL

G~=

(1/2)7r-2sgn(z3) 1m {

roo lq=Q

dq

00

T=TI(q)

exp( -sr)(8l

/8r)dr}.

(3.107)

Interchanging the order of integration yields (3.108) in which q = QI(r) is the (unique) inverse of r = T1(q) and is given by (3.109) Application of Lerch's theorem [54]on the uniqueness of the Laplace transform of a causal time function when taken at positive real values of s, leads to the time-domain result (3.110)

FIRST-ORDER

THIELE APPROXIMATIONS

81

c
SSP

= 750m/s;

= 0.450s

0.0

-448.0

500.0

Figure 3.3: Snapshot of the first-order Thiele approximated Green's function G~ in space for :1:3 2: O. (The horizontal offsets are numbered sequentially.) from which it is evident that TI(O) as given by Eq.(3.102) is the arrival time of the wave. Through the substitution (3.111) the integral in Eq.(3.1l0) can be evaluated in closed form. The result is

G2 =

41rC :1:3

I I sgn(:l:3)H[t

- T (0

(3.112)

(see Figure 3.3). An alternative way to derive Eq.(3.1l2), viz. with the aid of retarded potentials and the exact wave equation, is given in Appendix D. In a similar way, we evaluate the Green's function G~, given by Eq.(3.83). Taking the contributions from p = pI and p = »" together (note that, again,

82

THIELE-APPROXIMATED

WAVES IN FLUIDS

the integrand

satisfies Schwarz's reflection principle), Eq.(3.89) leads to p-11:VJI-t 1m { (ex> dq lex>

Jq:::o

G~,c = (1/2)7r-2

T:::TI(q)

exp( -sr)( r -lr)(

8pI / 8r )dr} ,

(3.113) where we have eliminated using Eq.(3.99). However, in addition to this we have a contribution from the arcs at infinity, leading to the direct source contribution. In the right-hand side of Eq.(3.113) we distinguish two terms, i, and given by il

t..

it and

= (1/2)7r-2p-ll:VJI-tlm{ i, = -(1/2)7r-2p-tl:VJI-lr

lq:::o

(ex> dqlex>

T:::TI(q)

exp(-sr)r(8pI/8r)dr}

(3.114)

1m { (ex> dq (00
lq:::o

IT:::TI(q)

exp( -sr)l(8pI

/or)dr}.
(3.115)

Interchanging it and

the order of integration

in both of them yields

lex>
IT:::TI(o)

exp(-sr)dr(1/2)7r-2p-tl:VJI-lr

Jq:::o

t:J (T)lm{8pI/8r}dq

(3.116)

i, = lex>
Application

T:::TI(O)

exp( -sr)dr

(-1/2)7r-2 p-tl:VJI-lr

Q1

q:::O

(T) Im{pI (8l /8r)}dq.

(3.117) results (3.118)

of Lerch's theorem once again leads to the time-domain

and

12

[-(1/2)7r-2p-ll:VJI-lr

.;:::(t) Im{l(8pI/8r)}dq]
it is found that

H[t - TI(O)]. (3.119)

With Eqs.(3.110)-(3.112)

11 =

47rpczJ

H[t - TI(o)].

(3.120)

SECOND-ORDER

THIELE

APPROXIMATIONS

83

Through the substitution Eq.(3.111), the integral in Eq.(3.119) can be evaluated in closed form as well. The result is (3.121) Taking II and 12 together, and adding the direct source term, we arrive at G~ = 47rpc :1:3
\

2 [et - (1'2/1:1:31)] H[t - TI(O)] + _!_5(:l:l,

2p

:1:2, :1:3) t2 H(t)

(3.122)

(see Figure 3.4, apart from the singularity). How the approximation of the vertical slowness affects AI,J, and hence the system of first-order differential equations, follows from Section 3.2 upon performing the steps in the reverse order. Equations (3.8) and (3.9) remain the same, but Eq.{3.10) changes upon replacing", by p-l [(-yI)2 - upup] into (3.123) and hence, in view of Eq.{3.94), into

s",[1

+ (1/4)c (upup)2]p
4

- isu/vp

+ 03V3 = q.

(3.124)

Upon replacing s by and -isup Eq.(3.124) is found as

4

at

by op, the space-time counterpart of

",[a: + (c /4)(opOp)2]

P + o:Orvr = a:q,

(3.125)

which replaces Eq.(3.2) in the first-order Thiele approximation.

3.6.

The second-order Thiele-approximated solutions in the space-time domain

The second-order Thiele approximation involves replacing "'I by "'Ill, and hence y by yll, where (cf. Eq.(2.44)) (3.126)

84

THIELE-APPROXIMATED

WAVES IN FLUIDS

c
SSP

= 7 SOmis;
-t

= 0.450s

horizontal offset (m) source location

0.0

-448.0

448.0

500.0

Figure 3.4: Snapshot of the first-order Thiele approximated Green's function G~ in space for :1:3 ~ O. (The horizontal offsets are numbered sequentially.)

SECOND-ORDER

THIELE

APPROXIMATIONS

85

and

yll = 3(pC)-1 [4/3~

A]«

+ a.... a....

+ aIJaIJ].

(3.127)

To construct the space-time domain counterparts of Gf, G~ and G~ we I I, again apply the modified Cagniard method. First, however, it is observed that upon replacing "I in Eq.(3.88) by "III according to Eq.(3.126), the inversion integral Eq.(3.88) becomes divergent since "III -t 3/c as laIJaIJI -t 00 and a physically meaningful interpretation of the integral must be found. The asymptotic structure of the integrand Eq.(3.88) as IaIJaIJ I -t 00 leads to the conjecture that the divergent integral contains a contribution that could be interpreted as a Dirac delta distribution operative at {Xl = 0,X2 = OJ. In view of this and of the contour deformation to follow, it is assumed that an interpretation in the sense of a Cauchy principal value around infinity will lead to a physically acceptable result. Adopting this interpretation, the integral on the right-hand side of Eq.(3.88) is replaced by the limit of the integral over the rectangle {al E 'R, a2 E 'Rj - Ll"'l Ll"'l -t 00 and Ll"'2 -t 00.

<

al

< Ll"'l' - Ll"'2 <

a2

< Ll"'2} with

The modified Cagniard method consists of applying the transformation of Eqs.(3.96) and (3.97), keeping q real. The resulting integral with respect to p is extended over {p E T; -iLl < p < iLl} with Ll -t 00. The integrand is continued analytically into the complex p-plane and the integration is carried out along the modified Cagniard path that follows from a continuous deformation of the imaginary p-axis and satisfies the equation (3.128) where T is real and positive. The expression for "III is written as (3.129) with

A2 = 8/3c2,
B2(q)

(3.130) (3.131)

q2+4/c2•

86

THIELE-APPROXIMATED

WAVES IN FLUIDS

The non real-axis part of the solution to the cubic equation (3.128) along which r strictly increases, is the desired modified Cagniard path. This part consists of p = plI(r, q) in the first quadrant of the complex p-plane, together with its complex conjugate p = pIl*(r,q) in the fourth quadrant of the complex p-plane. Its explicit form is obtained from Cardano's formula for the solution of the algebraic cubic equation; the relevant expression is given in Appendix B. In Section 2.4 the properties of the modified Cagniard path {p = pIl(r, q)Up = plI*(r, q)} have been investigated. In particular, it was found that p = pIl(r, q) and p = pIl*( r, q) are finite arcs that leave the real p-axis at pfI (q), at which point the value of r is denoted by T[I(q), and return to the real p-axis again at p~l(q), at which point the value of r is denoted by TfI(q). For GIl, as it follows from Eq.(3.81), the procedure has been carried out in Chapter 2, where it was found that

with

Tt,{ =

(2/V3)(r/c)

and

Tt,~ =

(4r/9Ix31)(r/c),

(3.133)

and in which q = QV(r) and q and r = T[l(q), respectively.

= Q~l( r)

are the (unique) inverses of r

= T[I(q)

Next, we consider G~l as it follows from Eq.(3.82). The contour deformation leading from the integration along the imaginary p-axis to the one along the modified Cagniard path is further analyzed. In view of Cauchy's theorem the integral along the imaginary p-axis is equal to the sum of the integrals along the semi-circle {p E Cj Ipl ~,Re(p) > O} and the modified Cagniard path. On the integrand is asymptotically equal to G~ = (1/2s2)sgn(x3)exp(-.!Ipr - 3.!1lx31/c).Using the asymptotic integrand

ct = ct

SECOND-ORDER

THIELE

APPROXIMATIONS

87

and carrying out the transformation inverse to Eqs.(3.96) and (3.97), it follows that the contribution of C! results in
Ll"l -~~"2 --+00

G~,

L~~~Ll"2

(s /27r)2

c:

dQl sgn(z3) exp( -3slz31/ c)da2

Z2)'

exp] -isQpzpJ(1/2s2)

= (1/2s2)

sgn(z3) exp( -3slz31/c)5{zl'

(3.134)

Taking, next, the contributions from p = pIl(r,q) and p = pIl*(r,q) together (note that the integrand satisfies Schwarz's reflection principle), the expression for of is obtained as

where also the property that the integrand is an even function of q has been used. Interchanging the order of integration yields

Application of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function when taken at positive real values of s, leads to the time-domain result

88

THIELE-APPROXIMATED

WAVES IN FLUIDS

c
SSP

= 750m/s

; t = 0.450s

0.0

-448.0

448.0

500.0

Figure 3.5: Snapshot of the second-order Thiele approximated Green's function in space for :1:3 ~ o. (The horizontal offsets are numbered sequentially.)

Gf

from which it is evident that t = Tfl{O) and of two wave fronts. Through the substitution
q=

TP{O) are the arrival times

(3.138)

Q~I(t)sin(1/1)

with

05:.1/1 5. 7r/2,

the inverse square root singularity at the upper limit of the first integration on the right-hand side of Eq.(3.137) is removed, while through the substitution

the inverse square root singularities at both limits of the second integration in Eq.{3.137) are removed. After substitution, the integrals are evaluated numerically with the aid of the trapezoidal rule. The result is shown in Figure 3.5.

SECOND-ORDER

THIELE

APPROXIMATIONS

89

c
SSP

= 750mjs;
-t

= 0.4508

horizontal offset (m) source location

0.0

-448.0

448.0

500.0

Figure 3.6: Snapshot of the second-order Thiele approximated Green's function G~I in space for :1:3 ~ o. (The horizontal offsets are numbered sequentially.) Next, we consider of as it follows from Eq.(3.83). Repeating the steps in Eqs.(3.134)-(3.137), we arrive at

from which it is evident that t = TfI(O) and t = TP(O) are the arrival times of two wave fronts. Through the substitutions Eqs.(3.138)-(3.139), the integrals are evaluated numerically with the aid of the trapezoidal rule. The result is

90

THIELE-APPROXIMATED

WAVES IN FLUIDS

shown in Figure 3.6. In addition, in the process of contour deformation we may encounter the simple pole p = pf of G{l at the simple zero of [-ylltl. From Eqs.{3.129) and (3.131) it follows that (3.141) Since the points of intersection of the modified Cagniard path with the real paxis satisfy the inequalities (see Section 2.4) pF (q) < B( q) < pF (q), the pole, as long as :1:3 =I 0, always lies inside the closed modified Cagniard contour and hence does not contribute to the result in the process of contour deformation. However, in the limit 1:1:31 --+ 0, the contour contracts to the point p = B(q) and G~I can be evaluated as the residue of the pole at B(q). This leads to

(3.142)

G3 (:l:J,:l:2,0,S)

°11

(3/2pc)(1/s 2 )6(:1:1,:1:2) __ 2_ roo exp( -sr) dr 7rpc3},r=2r/c [r2 - 4r2/c2J1/2 .

=

(3.143)

On account of Lerch's theorem on the uniqueness of the Laplace transform of a causal time function taken at positive real values of s, the time-domain counterpart of G~I then follows as

Gf

(3.144) in which (3.145) denotes the arrival time. How the approximation of the vertical slowness affects AI,]' and hence the system of first-order differential equations, follows from Section 3.2 upon

REFLECTION

AND TRANSMISSION

COEFFICIENTS

91

performing the steps in the reverse order. Equations (3.8) and (3.9) remain the same, but Eq.(3.10) changes upon replacing It by p-l [(-y1I)2 - ll:pll:p] into (3.146) and hence in view of Eq.(3.126) into (3.147) Upon replacing s by Ot and -isll:p Eq.(3.147) is found as
It

by op, the space-time counterpart of

[o?(o; - (c2/4)OpOp)2 + (c6/16)(ol'ol')3]p + Ot(o; - (c2/4)opOp)20rvr = Ot(o; - (c2/4)opOp)2q,

(3.148)

which replaces Eq.(3.2) in the second-order Thiele approximation.

3.7.

Exact, and first- and second-order Thiele-approximated reflection and transmission coefficients

In this section we investigate what influence the first- and second-order Thiele approximations of the spectral-domain vertical slowness pertaining to the acoustic wave motion in a fluid have on the reflection and transmission properties of a planar interface between two different fluids. To this end, we consider the wave motion in a configuration consisting of two adjacent half-spaces. In each of the two half-spaces a homogeneous fluid is present. In one of the halfspaces a source generates an acoustic wave motion. At the interface this wave is partially reflected and partially transmitted. The configuration is typical of a seismic exploration arrangement. To enable the wave matrix formalism of Section 3.2 to be employed, the source level is also treated as an interface.

92

THIELE-APPROXIMATED

WAVES IN FLUIDS

I Pl,Cl I

+
S

>

Xl.2

1 R+

T-

+
R

Figure 3.7: The configuration of the reflection problem.

REFLECTION

AND TRANSMISSION

COEFFICIENTS

93

The configuration is shown in Figure 3.7 and its nomenclature is given in Table 3.1. The source is located at :1:3= OJ the interface coincides with :1:3 h. = For the wave quantities we shall employ the acoustic-pressure normalization. In the two half-spaces with medium parameters {Pm, II':m}, = 1,2, we set (C£. m Eq.(3.20)) (m) ~ 12 F-[(m) - L-(m)W- J (3.149) [,J lor m = , , in which

w (m)
1

WJm)

= =

w~m) w~m)

exp[-S"Ym(:l:3 h)], exp[s"Ym(:l:3 h)]. -

(3.150) (3.151 )

Here, the vertical slownesses associated with the two half-spaces are indicated as ±"Ymfor m = 1,2. They imply the vertical acoustic wave admittances

Table 3.1: Nomenclature of the configuration. layer vertical range material properties vertical slownesses admittances

-00

<:1:3 < 0 0<:1:3 < h h <:1:3 < 00

PI, PI,

11':1, Cl 11':1, Cl 11':2, C2

P2,

-"Yl ±"Yl "Y2

Ym = p;;.I"Ym

for

(see Table 3.1). The boundary conditions of the continuity type and V3 lead to the condition (3.152)

By introducing reflection and transmission coefficients R±, T± according to (3.153)

94

THIELE-APPROXIMATED

WAVES IN FLUIDS

it is found in Appendix C that

R+ - - R- __

Vi - 1'1 Vi+1'1'
T- = 2

(3.154)

+_

- 2y;

1'1

2+

y;' 1

y.:2+ y;I

Vi

(3.155)

Note that T± = 1 + R±, which holds for all values of alA away from singularities of Ym• So far, we have not substituted any particular approximation of the vertical slownesses in the vertical acoustic wave admittance Y j hence, the relations are generally valid. In this respect, it is noted that an approximation of the vertical slownesses is to be applied to the delay factors and the reflection and transmission coefficients simultaneously, to avoid the introduction of artificial surface source distributions at the interface (see also Appendix C).

The exact case

In the exact case, we have to substitute (3.156) in Eqs.(3.154)-(3.155). The corresponding exact transmission and reflection coefficients are analytic away from the branch points that arise from the vertical slownesses in the different media.

The first-order Thiele approximation

In the first-order Thiele approximation we have (3.157) For this case the reflection and transmission coefficients become
R/'+ _ -R/'__

(P2C2)-1 (P2C2)-1

- (PlCltl + (PI Ct}-l

+ (1/2)[(c2/ P2) - (cd

+ (1/2)[( C2/ P2) + (cd

pdlalAalA pdlallall '

(3.158)

REFLECTION

AND TRANSMISSION

COEFFICIENTS

95

and

= =

2 (P2C2)-1 2 (P2C2)-l

(PICt}-l[l + (CU2)lllJlllJl + (PICt}-l + (1/2)[(c2/ P2) + {ct/ pdlll..,ll..,' (P2C2)-l [1 + (C~/2)lllJlllJl . + (PlCt}-l + (1/2)[(cdP2) + (ct/pdlll..,ll..,
(3.159)

The first-order Thiele-approximated reflection and transmission coefficients contain simple poles at the simple zeros of their common denominators. In the right half of the complex p-plane these are given by (note that lllJlllJ = q2 _ p2, cf. Eqs.(3.158)-(3.159) and Eqs.(3.96)-(3.97)) (3.160) in which (3.161) The following situations occur:

Ct/C2 < 1 then Ct/C2 = 1 then Ct/C2 > 1 then

(C)-2 < ci-2 , I (CI)-2 =Cl,-2 (CI >Cl·-2 )-2

(3.162) (3.163) (3.164)

In the case Cl = C2, the pole actually vanishes and the reflection coefficient reduces to the exact expression R1,+ = -(P2l - PIl )/(PII + P2l), while the same holds for the transmission coefficient T1,+ = (2PIl )/(PIl + P2l). (When, in addition, PI = P2, the reflection coefficient vanishes and the transmission coefficient becomes unity.)

The second-order Thiele approximation

In the second-order Thiele approximation we have

II _ II/ Ym - 1m

pm -

3(

pmCm

)-1 [4/3C!. 4/ 2

cm

+ lllJlllJ] + ll..,ll..,

for m = 1,2.

(3.165)