T R U N C A T E D A S Y M P T O T I C S O L U T I O N O F T H E O N E - D I M E N S I O N A L I N H O M O G E N E O U S W A V E E Q U A T I O N by B A R R Y C U R T I S Z E L J B.Se., The University of Vic to r ia , 1984 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A February 1987 © B a r r y Curt is Zelt, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Geophysics a n d A s t r o n o m y The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date February 25 1987 DE-6(3/81) ii A B S T R A C T I present a new time-domain method for solving for the stress and particle velocity of normally incident plane waves propagating in a smoothly varying one-dimensional medium. Both the Young's modulus E and the density p are allowed to vary smoothly with depth. The restriction of geometrical optics, that the wavelength be much less than the material stratification length, is not required in this new method. The infinite geometrical optics expansion is truncated after n terms, imposing a condition on the acoustic impedance I for exact solutions to exist. For the case n = 2 there are three general classes of impedance functions for which the resultant expansion is uniform and exact. To check the numerical! validity of the "truncated' asymptotic" (TA) solution results are calculated for the case of a medium with a linear velocity gradient for which there is an exact solution in the frequency domain. Since a linear velocity gradient is not one of the foregoing classes of impedance functions, a curve-fitting approach is necessary. The TA method compares favourably to the exact solution and is accurate to within the error of the required curve fit. Two classes of synthetic seismograms are calculated for smooth velocity and density variations. The same impedance as a function of traveltime is used for both classes. In the first class the principal variation in impedance is due to velocity, while in the second it is mainly due to density. The amplitudes in both classes of synthetic seismograms are very similar, but, as expected, the traveltime curves for each class are widely separated. For the case p = constant the TA solution is used as a bench-mark to analyse a two-term WKBJ approximation for three classes of velocity functions. The velocity functions are such that the TA solutions are exact. For two of the classes the WKBJ i i i solution performs well when the length of the transition zone is of the same order, or larger, than the length of the applied wavelet. For steeper velocity gradients the W K B J solution begins to differ significantly from the exact T A solution. The W K B J solution for the th i rd class performs extremely well even for steep gradients. Equations governing the validity of the W K B J solution are examined to explain the above results. Equations are derived to describe the distortion of a stress pulse propagating through a transit ion zone. For small velocity gradients (relative to the length of the applied pulse) the wavelet changes amplitude but its phase is not effected. As the gradient increases and the velocity function becomes a discontinuity at z — 0 the wavelet travels through undistorted. Only when the transition zone width is of the order of the length of the wavelet is there any visible phase distortion. Reflection and transmission coefficients as functions of time are calculated for low, intermediate and high gradient transition zones. The transmission coefficient is a delta function i n each case. The reflection coefficient has the shape of a Hilber t transform for low gradients. For higher gradients the reflection coefficient approaches the shape of a delta function. iv T A B L E O F C O N T E N T S A B S T R A C T . i i L I S T O F F I G U R E S v i A C K N O W L E D G E M E N T S v i i i C H A P T E R I I N T R O D U C T I O N 1 C H A P T E R II Q U A L I T A T I V E D E S C R I P T I O N O F T H E E Q U A T I O N S . 7 2.1 The Homogeneous M e d i u m 7 2.2 The Smoothly Varying Inhomogeneous M e d i u m 12 2.3 Exac t Solutions 16 C H A P T E R III R E S U L T S 21 3.1 Val id i ty of the Truncated Asymptot ic Solution 21 3.2 The T A Solution for Variable Velocity and Variable Density M e d i a . 26 3.3 Analysis of the W K B J Solution Relative to the T A Solution . . . . 31 3.4 Val id i ty of the Second Order W K B J Solution 44 3.5 Distor t ion of Stress Pulse Based on T A Solution 53 3.6 Reflection and Transmission Coefficients Calculated from the T A Solution 55 C H A P T E R I V C O N C L U S I O N S 64 R E F E R E N C E S 67 A P P E N D I X A Derivation of Equations (20) 70 V A P P E N D I X B Derivation of the Rayleigh Solution for Stress 74 A P P E N D I X C Derivation of the First Four Terms in the W K B J Solution for Stress 78 A P P E N D I X D Derivation of the Coefficients in the Exponential Series Form of the W K B J Solution for Stress 82 V I L I S T O F F I G U R E S 1. Stress pulse in a homogeneous medium 11 2. Al lowed forms of impedance for exact solutions 19 3. App l i ed stress pulse 23 4. T A solution for stress in an inhomogeneous medium 24 5. T A solution for particle velocity in an inhomogeneous medium 25 6. Absolute difference between Rayleigh, T A and W K B J solutions . . . . 27 7. Velocity-depth and density-depth curves for (i) p = u 1 / 4 , (ii) p = 3v1/4 . 29 8. T A solution for stress for case (i) p = v 1 ' / 4 ' 30 9. T A solution for stress for case (ii) p= Sv.1/4' 31! 10. T A and W K B J solutions for low gradient tanh velocity model! . . . . 33 11. T A and W K B J solutions for intermediate gradient tanh velocity model 35 12. T A and W K B J solutions for high gradient t anh 2 velocity model . . . 36 13. T A and W K B J solutions for low gradient co th 2 velocity model . . . . 37 14. T A and W K B J solutions for intermediate gradient co th 2 velocity model 38 15. T A and W K B J solutions for high gradient co th 2 velocity model . . . . 39 16. T A and W K B J solutions for low gradient t a n 2 velocity model 41 17. T A and W K B J solutions for intermediate gradient t a n 2 velocity model . 42 18. T A and W K B J solutions for high gradient t a n 2 velocity model . . . . 43 19. | 5 n + i / ( S ' r l w m a x ) | , n = 0 ,1 ,2 and \S3/LJ^^\ versus traveltime for low gradient t anh 2 velocity model 46 20. | < S , n + i / ( 5 „ w m a x ) | , n = 0 ,1 ,2 and | S , 3 / o ; 2 n a x | versus traveltime for intermediate gradient t anh 2 velocity model 47 v i i 21. \Sn+i/(Snumax)\, n = 0 ,1 ,2 and | 5 3 / w m a x | versus traveltime for high gradient t a n h 2 velocity model 48 22. iSn+iKSnUma.*)], n = 0 ,1 ,2 and | 5 3 / c j m a x | versus traveltime for high gradient t a n 2 velocity model 49 23. T A and W K B J solutions for high gradient t a n 2 velocity model and Rayleigh solution for best linear fit 50 24. T A and W K B J solutions for different t a n 2 velocity model and Rayleigh solution for best linear fit . 52 25. Scattering model ; . . . . 56 26. Incident, t ransmitted and reflected waves, transmission; andi reflection coefficients for low gradient t a n h 2 velocity model! . 59 27. Incident, transmitted and reflected waves, transmission and reflection coefficients for intermediate gradient t a n h 2 velocity model 61 28. Incident, transmitted and reflected waves, transmission and reflection coefficients for high gradient t a n h 2 velocity model 62 vi i i A C K N O W L E D G E M E N T S I thank my supervisor M a t t Yedl in for his constant guidance and friendship through-out the past three years. Its been a genuine pleasure working wi th M a t t and I must say that I feel honoured to have been his first graduate student. Thanks also to Dr . Br ian Seymour of the Department of Mathematics for making many helpful suggestions. In fact, he conceived the idea of the truncated asymptotic solution method. I also thank K e n W h i t t a l l for his help in the latter stages of this thesis and Sonya Dehler for allow-ing me to use some of her T g X files. To al l of the people here at the Department of Geophysics and Astronomy: Thanks for making this department the C l u b M e d of U B C . This thesis would not be complete without an extra special thank-you to my parents for their support over the years and to R o n and Lorraine M a c L e o d for making my stay in Vancouver even more enjoyable. This work was supported by a grant from Imperial O i l L t d . and a student scholarship from the Natura l Sciences and Engineering Research Counci l of Canada. O h yeah, T H I N G S R U L E ! 1 C H A P T E R I I N T R O D U C T I O N The purpose of this research is to study the propagation of waves in vertically inhomogeneous acoustic media. The results of this study can eventually be used to generate synthetic seismograms for transition zones found in continuous velocity and density logs. Currently, most techniques used to calculate one-dimensional ( l - D ) synthetic seis-mograms in a vertically inhomogeneous acoustic medium operate in the frequency do-main . These methods are the basis for the synthetic seismograms generated from well data. In 1955, Peterson et a l . , using a layered-earth approximation, investigated the synthesis of seismograms from well-log data. Such synthetic seismograms are then matched to the zero-offset seismic data, to extrapolate lithologic information. Berry-man et al . (1958) assumed that a continuous velocity log could be decomposed into a number of layers wi th linear gradients in velocity. Seismograms from multiple transi-t ion layers were calculated in the frequency domain, under the assumption of constant density. Wuensehel (1960) applied Thomson's layer matr ix method to calculate l - D synthetic seismograms, and he included multiples and transmission coefficients. He as-sumed constant velocity and density layers. Bortfeld (1960) analysed acoustic waves travelling in a transition layer, w i th constant density and the velocity changing linearly wi th depth. These results were extended by Scholte (1961), who allowed for arbitrary variat ion in velocity and density, and obtained a solution using successive approximation techniques. Brekhovskikh (1960) obtained exact solutions for the Helmholtz equation, using hypergeometric functions, while Gupta (1965) obtained similar solutions involv-ing Bessel functions of imaginary order, parameterised by the frequency. In addit ion, 2 G u p t a assumed a linear velocity variation as a function of depth. Foster (1975) inves-tigated transmission effects in a continuous, 1-D seismic model, in which he concluded that "transmission effects in a continuous 1-D seismic model are unsubstantial". Gray (1984) analysed Foster's approximation of the transmission coefficient by unity. He con-cluded that "transmission effects are not negligible in a case where impedance, though continuous, has large variations". Ganley (1981) used a modification of the commu-nication theory approach to include the effects of absorption and dispersion. He used a fiat-layered model, wi th each layer having constant velocity, constant density, and constant spatial attenuation factor. Pa i (1985) developed a modular solution process for the 1-D inhomogeneous wave equation in which the medium is divided into layers and vertically averaged within each layer. Bis method, which is a generalisation of the Haskell matrix method, requires solving an eigenvalue equation in the horizontal! coordinate and a wave equation i n the vert ical coordinate for each layer. The solution throughout the medium is obtained by matching layer solutions at layer interfaces. Work has also been done on synthetic vertical seismic profiles (VSPs ) . Dietr ich and Bouchon (1985) presented a numerical simulation of V S P s using the discrete horizontal wavenumber representation of seismic wave fields. They used the Thomson-Haskel l mat r ix formulation to propagate the discretised source wave field. Seymour and Varley (1970) and Varley et al . (1971) used the concept of the modulated simple wave to investigate both small and finite amplitude waves in weakly inhomogeneous and rate-dependent materials. This concept, which is a generalisation of the classical linear theory of geometrical acoustics, was extended by Seymour and Varley (1984), (1987) to derive a representation for the exact solution of the wave equation for a general class of functions, thus defining a depth variable velocity. There are a l imited number of exact solutions to the l - D inhomogeneous wave equation dz2 v2(z) These can be categorised in terms of the velocity function v(z). These include: 1. v2(z) = l + cz where c and d are constants. This case was first studied by Cans (1915). The solution of the wave equation in this case is expressed in terms of cyl indrical functions of the 1/3 order. 2. v2{z} = (1 + cz)m This case was studied by Wallot (1919) for arbitrary m. v(z) = cz + d. This linear velocity case was first studied by Rayleigh (1894, §148). The Rayleigh solution is derived explicit ly in Appendix B . 4 Th i s case was investigated by Rytow and Iudkevich (1946). are in terms of Hankel functions. Solutions in this case 5. A*) 1 — CiZ2 — C3Z3 — c4z4 — This case was studied by M u l l i n (1953). He found solutions in terms of parabolic cylinder functions. 6. c2{z) = d2 e e 1 - N 4Af 1 + e c (1 i e c z ) 2 - 1 This ease, studied by Epstein (1930), has solutions in terms of hypergeometric functions. 7. v{z) = de-cz. This case was investigated by Elias (1931) and Heller (1953). Solutions of the wave equation for this ease are Hankel functions. A much more thorough and detailed review of the existing solutions is given in Brekhovskikh (1960, §35). The work described in this thesis wi l l serve to extend this list. I present a new method for solving the 1-D inhomogeneous wave equation. The method is analytically exact for a large class of transition zones. Furthermore, the calculations are done directly in the time domain. The main advantage of this new method is that it does not require that the length scale denned by the material inhomogeneity be much greater than the applied wavelength. In addition, internal multiple reflections 5 and transmissions are summed as time evolves. Thus the divergence of the solution, as calculated by the standard ray-sum techniques of geometrical optics, is avoided. In chapter II I begin with the basic equations governing 1-D elastic wave propaga-tion. I then solve for stress o and particle velocity u for the case of a homogeneous medium. This analysis is used as a building block for the case when both the den-sity p and Young's modulus E are allowed to be smoothly varying functions of depth. In this case I begin with the initial conditions and boundary condition for the prob-lem being considered. I then proceed to outline the conventional approach to solving the problem and point out the weaknesses and inadequacies of this method. A new approach is suggested which solves the 1-D inhomogeneous wave equation exactly for specific classes of transition zones. This approach is based on truncating the infinite geometrical optics expansion after n terms. For the case n = 2 there are three general classes of impedance functions for which the resultant expansion is uniform and exact. These classes are stated explicitly and are general enough to allow curve fitting to other forms of transition zones. The "truncated asymptotic" (TA) solutions for stress and particle velocity are presented to end the second chapter. Chapter III begins with an analysis which shows that the TA solution is indeed accurate to within the error of the required curve-fit for transition zones which are not of the required form for an exact solution. I consider the case p = constant and velocity v to be a linear function of depth. For this case an exact solution exists in the frequency domain. The TA method is shown to perform well relative to the exact solution and a two term WKBJ solution. Results are also presented for the case of a medium with variable velocity and density. This is followed by an analysis of the second order WKBJ solution for stress relative to the TA solution. The purpose here 6 is to determine conditions under which the WKBJ solution is a valid approximation. Results are presented for three different velocity gradients for each of three classes of velocity functions. These classes are such that the TA solution is exact. Equations are presented which govern the validity of the second order WKBJ solution and the previous results are explained in terms of these equations. Equations describing the distortion of a stress pulse propagating through a transition zone are derived. Amplitude and phase distortion is discussed in terms of the velocity gradient. Finally, reflection and transmission coefficients as functions of time are calculated for transition zones with low, intermediate and high velocity gradients. 7 C H A P T E R II Q U A L I T A T I V E D E S C R I P T I O N OF T H E E Q U A T I O N S 2.1 The Homogeneous Medium The basic equations that govern waves propagating in a l-D elastic medium are Newton's Law, which can be written as du do dt dz' (1) and Hooke's Law do „ du — = E —. dt dz (2) In equations (l) and (2) o(z,t) represents the stress, u(z, t) the particle velocity, p the density and E Young's modulus. In general, both p and E are assumed to vary smoothly with depth z. Equations (l) and (2) can be combined as a first-order system of equations as follows d_ f u l = I" 0 1//>1 d_ U l . » dt [o\ [E 0 J dz [o\' Equation (3) describes all the basic physics relevant to l-D linear elastic wave propaga-tion. As a starting point I consider the case when E and p are constant so that the medium is homogeneous. Let u o 8 Then there exists a transformation matrix T which determines the trajectories along which linear combinations of u and a are constant. To find T, let u = T u' with u' « i u 2 ' (4) Then the original system (3) becomes dt T- 1 0 l/p E 0 du' = 0. (5) The system (5) decouples if T is chosen to be the matrix which diagonalises the matrix 0 1/p E 0 With T chosen in this way, u[ and u'2 are constant on the curves dz ~dt (6) where v is the velocity of wave propagation defined by v = y/E/p. The positive sigh corresponds to waves travelling to the right (or downward), while the negative sign corresponds to waves travelling to the left (or upward). Thus, this transformation matrix T decomposes the fundamental physical variables into their upgoing and downgoing wave components. The differential equations (6), which define the ray trajectories can 9 be integrated and parameterised by the constant of integration. Let a be the time intercept of the downgoing wave and /? the time intercept of the upgoing wave; a = - t 2 v / and P '-) v ) The curves defined by constant a and /? are known as the characteristic curves. These curves, which are straight lines in a homogeneous medium, can be determined by numer-ical integration in an inhomogeneous medium. I use either a or /? to define travelling waves since these parameters are constant for rays travelling in a smoothly varying medium. To demonstrate the usefulness of the foregoing parameterisation, I first explicitly calculate the transformation matrix T and its inverse T - 1 . These are given by \/2 j-l/2 j-1/2 _jl/2 j l / 2 (7) and T 1 = V2 / 1 / 2 _ / - l / 2 j r l / 2 J - l / 2 (8) where / is the acoustic impedance given by pv (also / = y/Ep). It then follows from equation (4) that r /1 1 r rl/2 r-l/2 1 (9) 1 • jl - J - l / " U jl/2 j - 1 / 2 a Combining equations (5), (7) and (8) gives + v-^r1 = 0 dt dz 10 where du\ du\ dz du\ dt dt dt dz by the chain rule. Thus, on alpha characteristic curves where dz/dt = +v, and which in this case are straight lines du\ dt Hence u\ (given in equation (9)) is an arbitrary function of a; ±ll^u-I-^o} = -F(«). (10) Similarly, u'2 is an arbitrary function of /?•; ±[i*'*u + rv**]=s{0}. (ii) Solving equations (10) and (11) for u and a in terms of F and S gives a = I^2[S{f3)+F(a)} and u — J - 1 / 2 [S((3) — F(a)] (12) with the \/2 absorbed in F and S. The solutions of equations (12) can be interpreted as follows. First, assume that S(/3) is identically zero. Then o = I1/2 F(ot). This is the source pulse applied to the 1-D medium. For any fixed a, a0, this source pulse amplitude is constant for all values of z and t satisfying the equation (13) 11 Figure 1. Solution for stress a as a function of time (seconds) and depth (kilometres). Plot shows a stress pulse propagating vertically in a homogeneous medium with constant velocity 3.5 km/s. As z increases, time must also increase for ao to remain constant. Therefore, the source pulse travels to the right without distortion through the medium. This model of l-D wave propagation is illustrated in Figure 1. A similar analysis may be applied to waves travelling to the left, by considering the curves /? equal to a constant. Waves travelling to the right are known as a-simple waves, while those travelling to the left are /9-simple waves. 12 2.2 T h e S m o o t h l y V a r y i n g I n h o m o g e n e o u s M e d i u m In what follows, both p and E are assumed to vary smoothly with depth z. The initial conditions are given by u(z,0) =a(z,0) =0, (14) with z ranging from zero to infinity. At z = 0, the boundary condition is given by the applied stress pulse a(0,t) b(t). (15) The a and (i characteristics are now given by * \(t T{z)) and 0=±(t + T(z)) (16) with To solve the system (3) with initial conditions (14) and boundary condition (15) I start by defining two new dependent variables: F(a,/3,I) = Fdown(«,I) + Fup((3,I) = l-rl/2{o-Iu), (18a) and 5(a,/?,/) = Sd own(a,/) + S u p 0 M ) = ^ / " 1 / 2 ( a +/u), (18b) 13 where I(z) = p(z)v(z) is the acoustic impedance. The subscripts in equations (18a) and (18b) correspond to the physical situation for which the applied stress pulse is principally downgoing or upgoing. For general initial conditions and boundary conditions, all four quantities Fd O Wn(«,^), Fup((3,I), Sdowi^ a,/), Snp(/3,I) must be calculated. For the boundary condition given by equation (15), I need only consider i*down arid Sdown-Henceforth, I drop the subscripts and write F(a,I) =• -/•• 1/2(<7 - lu) and S{a, I) = -I~1/2(o• + Iu). (19) 2 2 /'' and S in equation (19) correspond to the "fast" and "slow" Riemann invariants for the a-wave and are regarded as functions of a and l(z), with dl jdz normalised1 with respect to a small parameter e, the ratio of the applied wavelength to the material stratification length. (Thus, for example e = \dl/dz\/u> for a periodic input.) Note that the classical geometrical acoustics expansion is valid for e <C 1 (Brekhovskikh, 1960, §16). For a downgoing applied stress pulse at z — 0, travelling in a homogeneous medium, S is identically zero, while F is a time-shifted version of the applied stress pulse, and does not vary with depth. Therefore, S can be interpreted as a measure of the inhomogeneity of the medium, or the accumulated partial reflections for a given a and I. Similarly, F represents the accumulation of all the transmitted waves. With the change of independent variables from (t, z) to («,/), and the dependent variables from u and o to F and S, equation (3) reduces to the following first order 14 system of equations (see Appendix A) dF S + 21— = 0, (20a) dS d In/ / dS 1 \ , , , T« = ^ {'T, + 2F) • <20b> Equation (20b) shows clearly the dependence of S on the reflectivity |dln //dT. These equations may be solved by using asymptotic expansions (Seymour and Varley, 1970) of the following form: oo oo F(I,a) = £ V hn(a)fn(I) and S(I,a) = £ £B M a K ( J ) . (21) n=0, n=l These are substituted into equations (20) and coefficients of powers of e are matched to obtain the following hierarchy of equations valid for n > 1: + - A l i = 0, (22) dhn(a) da = hn. ,(a> (23): and Equations (22)-(24) are solved recursively with I(z) [and hence I(T)] and /io(£*) given. The initial equations which begin this hierachy are fo(I) = 1 and SQ(I) = 0. Then, for 15 example, hi(a) = / ho(x) dx, .,(«) = 5 S and /,(,) = - - / As the a-wave travels, energy is continuously transmitted and reflected. This pro-cess occurs continuously as the direct a-wave propagates through the medium. In the conventional mathematical analysis of the above situation, an infinite series is used to represent al l accumulated partial reflections and transmissions. In practice, however, only the first one or two terms are calculated, and presented as the solution to equation (3). Unfortunately, the foregoing expansions are invariably secular, being valid only for times or distances of order l / e . The secular behavior arises from the accumulation of internal reflections caused by stratification of the medium. The zeroth order solution contains only the dominant transmitted wave propagating in the direction of increas-ing 2, while the first order terms represent the accumulation of a l l double reflections. Eventually, as the signal is dispersed, the multiple reflections wi l l become of the same order as the primary signal, and the expansion w i l l no longer be val id. A multiple scale technique (Bender and Orszag, 1978, p. 556) can be used to extend the range of the validity of the expansion. However, this two-term series does not reproduce the physics of propagation correctly. This can be inferred by application of the following analysis. A s the zeroth order transmitted wave travels through the medium, it continually sheds reflected energy. Therefore, after it has propagated through the material for some time, 16 its amplitude will be of the same order as the partial reflections which have been ig-nored. Since it is not possible in practice to sum an infinite series of partial reflections and transmissions, another approach must be developed. A new method of solving the system (3) which avoids infinite series may be achieved via the following algorithm: a) Higher order transmitted waves are included in the calculation of the total en-ergy moving to the right. This modification to the conventional infinite series approach does not in itself alleviate the problem discussed earlier. b) Insist that the infinite series expansion terminate after n terms {n is usually 1 or 2). This imposes a condition on the acoustic impedance for an exact solution of the system (3) to exist. For the case n=2 there are three classes of acoustic impedance which satisfy the condition for the termination of the infinite se-ries. In practice, these three classes are sufficiently general to fit many material inhomogeneities. 2.3 Exact Solutions It has been shown in the previous section that the usual approach of taking a given, but slowly varying, impedance function and assuming an asymptotic solution in powers of e leads to seeularities and a limited range of validity. In this section I show that a different approach yields exact (and hence uniformly valid) solutions for a general class of impedance functions. I make no assumption on the form of I[z) but instead insist that the expansions (21) terminate after n terms. This yields a condition on I(z) for such an exact solution to exist. For the case n — 2 there are three classes of I(z) which 17 satisfy this condit ion, and these are sufficiently general to be able to curve-fit many forms of inhomogeneity. Fi rs t , I consider the a-wave, and again write o and u in terms of F(a, I) and S(a, I) using equations (19). I make no assumptions on the size of | / ' (^ ) | but assume that jP and S can be wri t ten as finite series of the form N N F = J2hn(a)fn(I) and S = >T hn{a)sn{I), (25) n=0 n=l where the hn, fn and sn again satisfy equations (22), (23) and (24). (This is equivalent to put t ing e = 1 in equations (21).) Since the series (25) terminate, equation (22) becomes SN{I} + 2 / ^ f = 0 (26); at n = N, and equation (24) becomes (,dsN 1 r \ at n = N + 1. The partial differentiation symbols may be dispensed wi th since both /pf. and sjy depend only on J . Equations (26) and (27) may then be combined into the following two Euler equations for fx and sjy: (27) IN Al^r + 4/2 SN ~ 4 / dl ds pj 41 d2fN dP ,d2sN dl2 = 0 18 The solutions of these equations determine and syy as ' N Ar'^ + Bl'/2 and sN = Ar1'2 - BI1'2 (28) where A and B are arbitrary constants. For the case of N = 1, a combination of equations (24) and (28) yields the following ordinary differential equation for /: ^ = 2 / ( A / - 1 / 2 - B / 1 / 2 \ . (29) dT v / Equation (29) is the condition imposed on the inhomogeneity for exact solutions to exist. Since the method of solution is exact, the usual restrictions of geometrical optics are not applicable. There are three classes of impedance functions satisfying equation (29) depending on the sign and magnitude of A and B. As a function of the traveltime T defined in equation (17), these are given by A A J(T) = - t a n h 2 [ V A B ( r + C)] for AB > 0 and I"< — (30a) I{T) = ~- coth 2 [\/AB{T + C)j for AB > 0 and / > ^ (30b) I(T) = - — tan 2 W-AB{T + C)] for AB < 0 (30c) B where C is a constant of integration. Figure 2 shows the allowed forms of the impedance functions for C > 0 and C < 0. In practice of course, the impedance of a medium will not exactly match one of the allowed forms. In this situation the impedance can be fitted to one of the forms given in equations (30), with the resulting solution accurate to within the error of the fit. 19 F i g u r e 2. Allowed forms of impedance J as given by equations (30). Plots on left are for C > 0. Plots on right are for C < 0. 20 The solutions for o and u can be writ ten explicitly for the case N = 1 by combining equations (19) and (25) to get a{a,I) = Ix'2 (h0[a) + (Sl{I) + h(I))hx{a)) , (31a) u(a,I) = r1'2 {-h0{a) + ( 5 l ( J ) - / 1 ( / ) ) / l l ( a ) ) (31b) wi th fx and Sx given in equation (28) and ho and hx related by equation (23). hx{a\ is obtained by substituting equation (23) into equation (31a), giving rise to the ordinary differential equation a(a,J) = I1/2(~^ + + M ' ) ) M " ) ) - (32) Equat ion (32) is solved by insisting that the surface (z = 0) be an absorbing boundary from which no backscattered energy can reflect. This is equivalent to requiring that o at z = 0 be equal to the applied stress pulse b(t) for all t > 0. This gives rise to the solution A»(•«)'•= ^TTT / ^xb{x)dx (33) 1 / 2 where 7 = fx{I'o) + ^ 1 ( ^ 0 ) = 2 A / / 0 is a constant and I0 = I[z — 0). /^o(° !) from equation (23) is then ho{«) = bM-lhl(a). (34) The foregoing calculations can be repeated for any order TV with I(T), sn, and / „ for n = l , 2 , . . . ,N — 1 determined from a system of 2N first order ordinary differential equations (22) and (24). 21 C H A P T E R III R E S U L T S 3.1 Validity of the Truncated Asymptotic Solution I present synthetic seismograms calculated with the T A series method. A s demon-strated in the previous sections, the solution obtained is val id for smoothly variable density p(z) and velocity v(z). Initially, however, I consider the case of a constant den-sity medium, so that all of the spatial variation in impedance is associated with velocity. In this si tuation the condition imposed on the inhomogeneity (29) becomes ~=2v{Av-l'2-Bv1'2) (35). where A = p - 1 / 2 A and B = plf2B are new constants. The three classes of velocity functions for an exact solution to exist are then given by the three equations defined in equation (30), wi th / replaced by v and A and B replaced by A and B. In addition to constant density, I consider the particular case of a medium with a variable velocity of the form v(z) = cz + d, or v(T) — decT wi th T the traveltime given in equation (17). For this combination of density and velocity, an exact solution (Rayleigh, 1894, §148) exists in the frequency domain (see Appendix B ) . Comparing the T A solution wi th this "Rayleigh" solution wi l l give a check of the numerical validity of the T A solution. The Rayleigh solution was chosen as the simplest exact canonical solution for normally-incident plane waves propagating in a smoothly varying l - D elastic medium. It should be stressed, however, that while the Rayleigh solution requires the assumption of constant density and a linearly variable velocity function, the T A method is valid for a broader class of inhomogeneities, (see equations (30)) and does not require the 22 assumption of constant density. For completeness the synthetic seismograms will also be calculated using a two-term W K B J solution (see Appendix C) and a comparison will be made of all three methods. For comparison purposes I calculate synthetic seismograms using a velocity function given by v{z) = {z + 3.5) km/s or v{T) = 3.5eT km/s. (36) There are several reasons for choosing this particular function. First, the surface veloc-ity, at z = 0, of 3.5 km/s is fairly typical of velocities encountered in seismic exploration. Second, the gradient chosen (1 s _ 1 ) is fairly steep (gradients larger than this are usually modelled by steps)1. In the following examples I use a shifted Ricker wavelet (scaled by the square root of the velocity at z = 0) of peak frequency / m a x = 10 Hz as the ap-plied stress pulse (15)'. Figure 3 shows this wavelet and its frequency spectrum. Figure 3b shows that this wavelet has no energy at 0 Hz. Having chosen the velocity model! given by equation (36), so that the exact Rayleigh solution may be used, I curve-fit dv(T)/dT to the form given in equation (35). Carrying out this procedure using a linear least squares fit over 100 points for 0 < T < 0.5 s (or 0 < z < 2.2 km) resulted in an I2 error of approximately 0.0085 between v(T) — 3.5er and the allowed form. For the velocity function given in equation (36) AB < 0 so that the allowed form for this case was the tangent squared form (30c). Figures 4 and 5 show the result of propagating the wavelet of Figure 3a vertically through an inhomogeneous medium with constant density whose velocity varies as in equation (36). Figure 4 shows the T A solution for stress as a function of time and depth. Figure 5 shows the T A solution for particle velocity. These figures show that the initial 23 1 L I L 0.20 0.30 Time (seconds) 0 10 20 30 40 Frequency (Hz) F i g u r e 3. (a) Applied stress pulse b(t) at z = 0 (Ricker wavelet), (b) Fourier transform of Ricker wavelet. 24 F i g u r e 4. T A solution for stress a as a function of time (seconds) and depth (kilome-tres). Plot shows the solution for a stress pulse propagating vertically in an inhomoge-neous medium with v(z) = z + 3.5 km/s. The density of the medium p is constant. stress and particle velocity pulses propagate essentially undisturbed except for changes in amplitude. At the scale of these plots, it is impossible to detect any differences between the Rayleigh, W K B J and T A solutions for either stress or particle velocity. This is to be expected because in the T A solution the single largest source of error arises in the curve-fitting of dv{T)jdT to the allowed form for an exact solution (35). The error in this case was very small (0.0085) and thus the T A solution should be expected to be extremely close to the exact solution. 25 Figure 5. TA solution for particle velocity u as a function of time (seconds) and depth (kilometres). Plot shows the solution for a particle velocity pulse propagating vertically in an inhomogeneous medium with v{z), = z + 3.5 km/s. The density of the medium p is constant. The second order WKBJ solution is very close to the exact solution because the dominant frequency of the input wavelet is relatively high (10 Hz). At this frequency, the ratio of the applied wavelength to the material stratification length is 1:10. For much lower frequencies, when the applied wavelength is of the same magnitude or greater than the material stratification length, the WKBJ solution will fail. The accuracy of the TA solution, however, is independent of the frequency of the input pulse. To estimate more accurately the difference between the exact Rayleigh solution and the two other solutions, plots were obtained showing the point by point absolute 26 difference between the methods at a depth of z = 2.0 km. Figure 6a shows the difference between the Rayleigh solution and the T A solution (/2 error over 200 points is 4.7 x 10~ 4), while Figure 6b shows the difference between the Rayleigh solution and the W K B J solution (I2 error over 200 points is 1.1 x 10~ 5). It is clear from these plots that the second order W K B J solution is the more accurate of the two approximation methods. The I2 error is much smaller and the spikey difference plot indicates that the error between the solutions is of the order of the round-off errors accumulated in the calculations. Figure 6a shows a much smoother difference curve indicating a larger systematic error in the T A solution. This error is due to the error in the curve-fit. As the velocity gradient increases, the error in the curve-fit is expected to increase. Since these results were calculated using a fairly high velocity gradient the results shown for the TA solution are essentially a worst case. Examining the magnitude of the differences, however, indicates that the error between the exact solution and the T A solution is still very small and well within the range of acceptability. 3.2 The T A Solution for Variable Velocity and Variable Density Media In the preceding examples the T A solution was compared with a canonical solution for the case of constant density, with the variation in impedance solely due to the variation in velocity. This served as a check of the numerical validity of the TA solution. I now consider waves travelling in a medium where both density and velocity vary. For an impedance model I use equation (30b) so that the TA solution will be exact (no curve fitting is required). This function represents a simple monotonic transition zone. 27 CO 1 O 80 -/ \ I2 error = 4.7 X 1 0 " 4 X 60 Diff. 4 0 Abs. 2 0 0 i | 1 i i i i 0 .60 0.70 0.80 Time (s) F i g u r e 6. Point by point absolute difference between the exact Rayleigh solution and (a) the T A solution, (b) the second order W K B J solution for stress a at a depth of z = 2.0 k m . 28 Also, I assume a density-velocity relation of the form (Sheriff and Geldart, 1983) p = avl/A, (37) where a is a constant. For the same impedance, given by equation (30b), I calculate two different velocity-depth functions corresponding to two different density-velocity functions. In Figure 7a I present the velocity-depth functions for the two density-velocity relations: (i) p = v 1/ 4; (ii) p = 3u 1 / / 4. Figure 7b shows the two corresponding density-depth functions. The impedance function (equation 30b) for these two cases are identical in the traveltime do-mairu, the difference being that in case (ii), more of the spatial variation in the impedance is associated with the density. Figures 7c and 7d show the velocity and density gradients for; the two cases. Figure 7c shows a larger velocity gradient for velocity function (i)>. Figure 7d shows that the density gradient for the density function (ii) is much larger at shallow depths, while at a depth of approximately 0.67 km the two gradient curves cross. This intersection occurs because density function (ii) approaches an asymptotic value faster than (i). Otherwise, no special significance is attached to this. In Figures 8 and 9, I present the results for a for the two cases described above. Figure 8 corresponds to case (i) while Figure 9 corresponds to case (ii). A comparison of these two plots shows the major difference to lie in the time of the arrivals of the stress pulses. In case (ii), the velocity is much slower, delaying the arrivals as seen in Figure 9. 29 0:5 0 .0 0.4 0.8i Depth- (km) 4.0 E 4> 3 . 0 w 2 . 5 2? m 2 . 0 c «) Q 1.5 1.0 \ (i) (b) . i 0.0 014! 0.8 Depth (km) 0.4 0 .8 Depth (km) o.o 0.4 0 .8 Depth (km) Figure 7. (a) Velocity-depth curves for two different density-velocity relations: (i) p = v 1 / / 4, (ii) p = 3V1/4. The impedance as a function of traveltime is the same for both cases, (b) Corresponding density-depth curves, (c) Corresponding velocity gradient curves, (d) Corresponding density gradient curves. 30 F i g u r e 8. T A solution for stress ( r a sa function of time (seconds) and depth (kilometres) for case (i) p = v1/4. This is an exact solution. Al though the amplitude of the pulses appear to be identical 1 on both plots they are different. This is because the impedance functions (which control the amplitude) corresponding to the two cases are identical at equal traveltimes only. For two different velocity functions at a given depth the traveltimes are different, and hence, so are the amplitudes. Only at z = 0, where T = 0 do the amplitudes on the two plots match. 31 F i g u r e 9. T A solution for stress a as a function of time (seconds) and depth (kilometres); for case (ii) p = ZvllA. This is an exact solution. The variation in impedance is due mainly, to density variations relative to case (i)!. 3.3 Comparison of T A and W K B J Solutions In this section I compare solutions for stress as calculated by the T A method and the second order W K B J approximation method. I take p = constant for simplicity and v to be of the form given in equations (30) (with / replaced by v and A and B replaced by p-ll2A and pl/2B). Thus, the T A solution is exact and I will be able to determine conditions under which the W K B J solution of (3) fails significantly. The second order W K B J solution consists of the first three terms in the WKBJ expansion for the solution of (3). The solution for stress (derived in Appendix C) is 32 given by a{u,T) ^a{u,0) [F1{T) + (38) where 1 / 2 F0{T) = r , (39a) and F2(T) =AFl(T) (v-ll2{T) - 1 / 2 S V0 } (39b) and uo = v{T = 0). I now consider three examples for each of the three forms of velocity functions. The three examples vary only in terms of the velocity gradient. Equivalently, I could have varied the frequency of the applied stress pulse and achieved the same results^ In the first example I consider a velocity function which starts at V'Q = 3 km/s when z = T = 0 and asymptotically approaches a value of v x = 4 km/s (Figure 10a)'. The length of this transition zone is ~3.8 km or, in terms of traveltime, 1 second, i.e., in this example the transition zone width is about 1 order of magnitude larger than the width of the applied stress pulse l / / m a x - Figure 10b shows the TA and second order WKBJ solutions for stress at T = 1.0 seconds, while Figure 10c shows the absolute difference between the two solutions. In this case the applied stress pulse has travelled through the medium essentially undisturbed. The peak amplitude of the wave has increased from vlJ2 to vlJ2. The £ 2 error (over 200 points) between the two solutions is (i) v ~ tanh 2 33 4.0 \ 3.8 E 3.6 3.4 o > 3.2 3.0 2.0 1.5 1.0 0:5 0 0 -0.5 -1.0 A = 4.89 B = 1.22 C - 0.538 (<0 0 1 2 3 4 Z (km) TA r I \ ' WKBJ 2 i i i r \ 1 ! Jl \ / (b) • v i' ' i i 1.10 1.20 1.30 Time (s) 0.0010 a) g 0.0008 3= 0.0006 Q O. 0.0004 8 0.0002 X> < 0.0 1. 10 V I?' error = 0.0057 1.20 1.30 Time (s) Figure 1 0 . (a) Hyperbolic tangent squared velocity model, (b) T A and W K B J solu-tions for stress a as a function of time at traveltime T = 1.0 seconds. In this example the two solutions overlay one another, (c) Absolute difference as a function of time between T A and W K B J solutions of (b). 34 ~0.0057 indicating very good agreement between the second order W K B J solution and T A (exact) solution. In the second example I increase the velocity gradient by a factor of 10 so that the transition zone width is now of the same order as the length of the applied stress pulse (Figure 11a). This is equivalent to decreasing the peak frequency of the applied stress pulse by a factor of 10 in the previous example. Figure l i b shows the T A and WKBJ solutions for o at T — 0.1 seconds. Examination of the sidelobes shows that the wavelet has been slightly distorted after travelling through the inhomogeneity. Comparing the two solutions shows a small but noticable difference. The I 2 error (see Figure 11c) in this ease is still reasonably small (~0.39). In the third example I again increase the velocity gradient by a factor of 10 so that the transition zone width is about 0.01 seconds, i.e., one order of magnitude smaller than the width of the applied stress pulse (Figure 12a). The two solutions (Figure 12b) in this case are markedly different. The T A solution shows a wavelet which has propagated through the medium without undergoing phase or amplitude distortion. The W K B J solution, on the other hand, shows severe phase distortion and in no way resembles the correct (TA) solution. The £ 2 error (see Figure 12c) in this case is very large (~6.9) indicating major failure in the W K B J solution for this high gradient case. The validity of the W K B J solution will be discussed further in Section 3.4. (ii) v ~ coth 2 This function is similar to the hyperbolic tangent squared function, and can be used to model zones of decreasing velocity. In the examples to follow, the velocity will start at v0 = 4 km/s and drop to vt ~ 3 km/s. Figures 13-15 show the results for three different 35 0.30 0.40 Time (s) Figure 11. (a) Hyperbolic tangent squared velocity model. Transition zone width is 10 times smaller than model in Figure 10a. (b) TA and WKBJ solutions for stress a at T = 0.1 seconds, (c) Absolute difference between TA and WKBJ solutions of (b). 36 I* i i f t 0.20 0.30 Time (s) 3 1.0 0.8 i / V / V I 2 error = 1 V I V 6 - 9 0.6 l 1 \ 0.4 - / 0.2 0.0 1 J _L 0.20 0.30 Time (s) F i g u r e 12. (a) Hyperbolic tangent squared velocity model. Transition zone width 10 times smaller than model in Figure 11a. (b) T A and W K B J solutions for stress a T = 0.01 seconds, (c) Absolute difference between T A and W K B J solutions of (b). 37 4.0 ^ 3.8 E - * 3.6 .t? 3.4 o _0 » 3.2 3.0 A = 4.23 B = 1.41 - \ C = 0.538 : i i i (o) < i 0.0 1.0 2.0 3.0 Z (km) 1.5 j / V TA 1.0 / V WKBJ2 0.5 1 \ 0.0 f | It ' | 1 f 0.5 [ v i i \ / W 1. 10 1.20 1.30 Time (s) 0.0006 o 0.0005 c (D 0.0004 H— 5 0.0003 •2 0.0002 » 0.0001 < 0.0 1.10 I 2' error = 0.0033 1.20 1.30 Time (s) Figure 13. (a) Hyperbolic cotangent squared velocity model, (b) TA and WKBJ solutions for stress a at T = 1.0 seconds. In this example the two solutions overlay one another, (c) Absolute difference between TA and WKBJ solutions of (b). 38 o 4.0 3.8 3.6 . A ^ B = 42.3 = 14.1 3.4 \ C = 0.0538 3.2 3.0 (a) i i i i i i i 0.0 0.10 0.20 Z (km) 0.30 1.5 L fx 1 X It V TA T.O 1 \ I WKBJlj 0.5 1 ! V V 0.0 . — \ h \ \Ji \ i i: i 0.5 { i • / (b) i i | ; 0.030 » 0.025 S 0.020 a> 0.015 -- § 0.010 v> „ , 0.30 0.40 Time (s) - / \ I? error = -If ' li 0.24 -It V (c) r > i i . . . . i I 0.30 0.40 Time (s) F i g u r e 14. (a) Hyperbolic cotangent squared velocity model. Transition zone width 10 times smaller than model in Figure 13a. (b) T A and W K B J solutions for stress a T = 0.1 seconds, (c) Absolute difference between T A and W K B J solutions of (b). 39 F i g u r e 15. (a) Hyperbolic cotangent squared velocity model. Transition zone width is 10 times smaller than model in Figure 14a. (b) TA and WKBJ solutions for stress a at T = 0.01 seconds, (c) Absolute difference between TA and W K B J solutions of (b). 40 velocity gradients. The transition zone width for these three examples are 10// m a x, l//max a n a " l/(l0/max) respectively as in the previous set of examples. The results obtained here are also similar to the previous examples—the W K B J solution performs very well for low gradients and very poorly for high velocity gradients. Further, for low 1 / 2 1 / 2 velocity gradients the amplitude of the stress pulse decreases from v0 to vx while for 1 / 2 high gradients the amplitude does not vary from v0 . The principal difference is seen in the distortion of the stress pulse relative to the applied pulse for the intermediate gradient case (cf. Figures 14b and l i b ) . ( i i i ) v ~ t a n 2 I now present results for stress using the third allowed form of velocity. Figures 16a, 17a and 18a show the velocity models used. These functions are essentially linear over the range displayed, increasing from VQ = 3 km/s to vt — 4 km/s over three different widths. The widths are the same as those used in the previous examples (i.e. 10//jnax) l/Zmax, 1 / ( l 0 / m a x ) respectively). In all three eases there has been virturally no distortion of the applied stress pulse (see Figures 16b, 17b and 18b). The W K B J solution is extremely accurate in all three examples, with the largest i2 error being ~ 1.1 X 1Q - 7 for the high gradient case. The spikey difference curves (Figures 16c, 17c and 18c) indicate that the difference between the two solutions is exclusively due to round-off error. To summarise, the second order W K B J solution performs well relative to the exact T A solution when the velocity function is of the form tanh 2 and coth 2 and when the length of the transition zone is of the order of l / / m a x or greater. For higher velocity gradients the W K B J solution begins to break down. For the case v ~ tan 2 the W K B J 41 CO o CD 4.2 A = 3.72 X 10"5 4.0 - B = -0.0773 3.8 - C = 919. 3.6 -3.4 -3.2 - (o) 3.0 " i • i i • 0 1 2 3 4 Z (km) 2.0 1.5 1.0 N : 0.5 0.0 -0;5 -1.0 1.20 1.30 Time (s) 0.10 03 II 0.04 (0 0.02 < 0.0 -I I? error = -; 1.1 X 10~9 (c) 1 1 1 1 1 1.10 1.20 1.30 Time (s) F i g u r e 16. (a) Tangent squared velocity model, (b) TA and W K B J solutions for stress a at T = 1.0 seconds. In this example the two solutions overlay one another, (c) Absolute difference between T A and W K B J solutions of (b). 42 4.2 A = 3.72 X 10"4 4.0 - B = -0.773 y 3.8 C - 91.9 / 3.6 3.4 -3.2 - (o) 3.0 - f i i i i 0.0 0.1 0.2 0.3 Z (km) 0.4 -1 .0b I L L L L 0.30 0.40 Time (s) 0.10 00 if 0.08 ; I 2 error = O ! 1.4 X 10"fl X 0.06 v — ' vi—' 0.04 Ci CO 0.02 X ) ( c ) < 0.0 i i i i i 0.30 0.40 Time (s) F i g u r e 17. (a) Tangent squared velocity model. Velocity gradient of this model is 10 times greater than model in Figure 16a. (b) TA and W K B J solutions for stress o at T = 0.1 seconds. In this example the two solutions overlay one another, (c) Absolute difference between T A and W K B J solutions of (b). 43 10 O 8 >S 6" CO JO; < 2 -I 2 error, = 1.1 X 10~7 u s u m . 0.20 0.30 Time (s) (c) F i g u r e 18. (a) Tangent squared velocity model. Velocity gradient of this model is 10 times greater than model in Figure 17a. (b) T A and W K B J solutions for stress o at T = 0.01 seconds. In this example the two solutions overlay one another, (c) Absolute difference between T A and W K B J solutions of (b). 44 solution appears to perform well, even for steep velocity gradients (though it will even-tually fail for extremely steep gradients). It should be noted, however, that while the tanh 2 and coth 2 velocity functions can be considered as transition zones, the tan 2 func-tion is not a "true" transition zone in that it does not level off to some constant value. This feature of the tan 2 velocity model will be discussed further in the following section. 3.4 V a l i d i t y o f t he S e c o n d O r d e r W K B J S o l u t i o n In an effort to explain and understand the performance of the second order W K B J solution I now consider the conditions under which it is a valid approximation. First I rewrite the W K B J approximation for stress (38) as an exponential series OJV(U>,T) ~ a(w,0)exp iuy^(iu)"nSn(T) n=0 (40) This series will be a good approximation of a if the condition (Bender and Orszag, 1978, p. 494) Rr « i M . o;-»oo (41) is satisfied for all T in the interval of approximation. Here R is the relative error between the asymptotic approximation ON(T) and the exact solution. For the case N = 2 equation (41) requires R = « 1 (42) '''max where c j m a x = 27r / m a x and / m a x is the peak frequency component of the input wavelet. In addition, cr2{T) will be a valid approximation on an interval if each of the functions 45 Sn+i(T)/Sn(T) (n = 0,1,2) are bounded functions of T on the interval. Equivalently, this requires that the asymptotic relations | S n + 1 ( T ) | < \wSn{T)\ n = 0,1,2 as w oo (43) hold uniformly in T. To find So,..., S3 I equate the third order W K B J expansion (given in equation (C-2) of Appendix C) to 03(T) (see Appendix D), obtaining S 0 T, Si - In F\, F-2 F3 1 / F 7 \ 2 S 2 = and S 3 = ' - - { -± i where F:i(T), the first term neglected in the two term expansion (38), is derived in Appendix C and is given by M n = -A % 1 / 2 Fa{n (45), Figures 19a,, 19b and 19c show |S n+i/(S nu; m a x)| versus traveltime for n = 0,1,2. These curves correspond to the first example discussed previously (i.e., v ~ tanh 2, low gradient) in which the W K B J solution very accurately approximated the exact T A solution. On the interval of interest (0 < T < 1.2) \Sn+i/(Sn0Jma,y.) \ <C 1 for n = 0,1,2 so that the W K B J approximation (40) will be valid on the interval. Figure 19d shows that R (defined in equation (41)) as a function of traveltime is bounded by 0.00025 indicating that a very accurate two-term W K B J representation is to be expected. Recalling that R gives the relative error between the W K B J and T A solutions (44) 46 — 0.010 0.008 (/^ 0.006 <^ 0.004 0.002 A - 4.89 B - 1.22 C - 0.538 • • (o) 1 0.0 0.4 0.8 1.2 Travel Time (s) 0.4 0.8 Travel Time (s) 0.4 0.8' Travel Time (s) 0.6 0.4 o.8; Travel Time (s) (<0 T.2! Figure 19-. (a-e); | 5 n + i / ( 5 n w m a x ) | versus traveltime T for n = 0,1,2'. (d) R = l^s/Umaxl versus traveltime T. These plots correspond to the low gradient, hyper-bolic tangent squared velocity model example of Figure 10. The terms Sn, n = 0,1,2 are related to the first four terms in the WKBJ expansion for a by equations (44). means that the maximum relative error for the range of traveltimes indicated is about 0.025%. This analysis agrees very well with the results shown in Figure 10. Figures 20 and 21 show the same curves for the intermediate and high gradient examples of Figures 11 and 12. These figures show that decreasing the transition zone width by a factor of 10 results in a corresponding tenfold increase in the amplitude of |5 n +i/(5 nu; n- i a x)|, (n = 0,1,2). The amplitude of ISS/W^^I increases by a factor of 102. These effects follow directly from equations (44), (45) and (39). For the intermediate gradient case 47 — 0.10 J 0.08 t/f 0.06 CO 0.04 0.02 A - 48.9 B - 12.2 C - 0.0538 (o) Li i_ i i 1 1 u 0.0 0.04 0.08 Travel Time (s) 0.12 0.0 0.04 6.08) Travel! Time (s) 0.12: a „ 0.435 to ^0.430 CO 0.425 0.420 -0.04 0.08 0.12 Travel Time (s) 0.0 0.04 0.08 0.12 Traveli Time' (s) F i g u r e 20. (a-e) iS„ +i/(S nw max)|; versus T for n = 0,1,2. (d) R = |53/wmax|: versus T. These plots correspond to the intermediate gradient, hyperbolic tangent squared velocity model example of Figure 11. (see Figure 20) |S" n + 1 /(5„u; n K l x)| (n = 0,1,2) are not <C 1, but only < 1. As well, l^/^maxl ' s bounded by 0.025 indicating a relative error of 2.5%. It is not surprising then that a slight degradation in the accuracy of the W K B J solution is seen in Figure l i b . For the high gradient ease (see Figure 21) | 5 n + 1 / ( 5 n w m a x ) | > 1 (n = 1,2) while R is bounded by 2.5 indicating a relative error of 250% between the W K B J and TA solutions. This extremely large discrepancy is seen clearly in Figure 12b. 48 — 1.0 J^O.8 (/? 0.6 CO 0.4 0.2 A - 489.05 8 - 122. C - 0.00538 -(o) "l l_ i i i 0.0 0.004 0.008 0.012 Travel Time (s) 0.0 0.004 0.008 0.012* Travel: Time (s) 0.0 0.004 0.008 0.012 Travel Time (s) 0.004 0:00ffi 0.012 Travel Time (s) F i g u r e 21 . (a-c) | 5 n + 1 / ( 5 n w m a x ) | versus T for n = 0 ,1 ,2 . (d) R = ^ / w ^ J versus-T. These plots correspond to the high gradient, hyperbolic tangent squared velocity model example of Figure 12. The results for the case t; ~ c o t h 2 are analogous to the above results. The relative errors were 0.02%, 2% and 200% for the low, intermediate and high velocity gradient cases respectively. Figure 22 shows \ S n + i / ( S n u > m 3 X ) \ , (n = 0,1,2) and | « S 3 / W m a x | versus traveltime for the high gradient, tangent squared velocity model example of Figure 18. For this case | S i / ( 5 0 u ; m a x ) | is < 1 but not <C 1, while | 5 „ + i / ( 5 n u ; m a x ) | <C 1 for n = 1,2. Also , l ^ 3 / w m a x l * s bounded by 1.8 x 1 0 ~ 1 0 so that the maximum relative error over the range 49 F i g u r e 22 . (a-e) |.S'n, i / ( 5 n w m a x ) j versus T for n — 0,1,2. (d) i? = \S3/u)fu:iX\ versus T. These plots correspond to the high gradient, tangent squared velocity model example of Figure 18. of interest is 1.8 x 1 0 _ 8 % , agreeing with the error curve of Figure 18c. In Figure 23a I increase the velocity gradient by another factor of 10. Figure 23b shows the TA and W K B J solutions for this velocity model, while Figure 23c shows the absolute difference between these two solutions. Again, the error between the two solutions is extremely small. This tenfold increase in the gradient increases the values of |S„4-i/(5 riu;max) |, (n = 0,1,2) by a factor of 10, and i S s / w ^ J by a factor of 100. Thus, in this case \Si/(S0u>ma.*)\ > 1 over the range of interest (0 < T < 0.0012) while \Sn+i / (Snojm^) \ <C 50 2.0 1.5 y — N 1.0 N i 0.5 0 o \J • w -0:5 -1.0 1.0 1 1 O 0.8 T ~ X 0.6 w H— H— 0.4 Q CO 0.2 < 0.0 0.20 0.30 Time (s) l2 error = 5.2 X 10r 7 0.20 0.30 Time (s) F i g u r e 23 . (a) Tangent squared velocity model. Velocity gradient of this model is 10 times greater than model in Figure 18a. (b) T A and W K B J solutions for stress o at T = 0.001 seconds. Also shown is the Rayleigh solution for a velocity model which is the best linear fit to the model shown in (a), (c) Absolute difference between T A and W K B J solutions of (b). £ 2 error between T A and Rayleigh soutions is 1.5. 51 1 and |5'3/u;2nax| is bounded by 1.8 x 10~ 8. It is evident that the condition for n = 0 in equation (43) is not necessary for an accurate approximation. This condition is required to be satisfied only by the highest order terms included in the expansion. The velocity function in Figure 23a appears to be nearly linear over the range displayed. To further investigate the tan 2 velocity model I calculate the solution for stress using a velocity function which is the best linear fit to the function in Figure 23a. Carrying out this fit over the range 0 < z < 0.0035 resulted in the linear velocity model v(z) =• (290.9z + 2.943) km/s. Using this v{z) I computed the Rayleigh solution (see Appendix B) for stress shown in Figure 23b. This solution (which is exact for linear velocity functions) differs substantially with the T A and W K B J solutions also shown in Figure 23b {£2' error between the Rayleigh and T A solutions is ~ 1.5). In Figure 24a I present another velocity curve of the tan 2 variety which,, over the range displayed, is very similar to Figure 23a. Both curves are nearly linear functions increasing from 3 km/s at z = 0 to 4 km/s at z = 0.0035 but the curve of Figure 24a increases more slowly for larger z. Figure 24b shows the T A and W K B J solutions for such a velocity model. This figure also shows the Rayleigh solution for the linear velocity model mentioned previously. In this case the T A and Rayleigh solutions agree fairly well (P error ~ 0.17), while the W K B J solution differs substantially (E2 error between W K B J and TA solutions is ~ 3.4, see Figure 24e). So the T A method, which is exact for the tan 2 velocity model, gives two very different solutions for two velocity functions (Figures 23a and 24a) which are very similar over the range of interest (i.e., 0 < z < 0.0035). It is clear that one must take into account the entire length of the velocity function being considered when comparing results for v ~ tan 2. This is necessary because the T A solution is exact and includes all 52 o 4.2 A = 248. 4.0 " B = -5.77 3.8 C = 0.00683 y 3.6 3.4 3.2 («) 3.0 " i i i i 0.0 0.002 0.004 0.5 F <D O c 0.4 t-a> £ 0.3 Q a> 0.2 8 o.i JO < o.o Z (km) 2.0 A ' A 1.5 ' / V — 'A V, '/ V, -1.0 II: V I f i — .'/ V N l 0.5 III J, ll ft V — • ' I ; Hi if; ill 0.0 // { ^ \ 1 | -0.5 V / \ -1.0 \ \ — TA — WKBJj — RAYLEIGHl 0.20 0.30} Time (s) -/ \ I 2 error = -/ V / \ 3.35 -1 V i i ) 0.20 0.30 Time (s) F i g u r e 24. (a) Tangent squared velocity model. Over the range displayed this model is very similar to the curve shown in Figure 23a but remains more linear over a larger range of depths, (b) TA and W K B J solutions for stress a at T = 0.001 seconds. Also shown is the Rayleigh solution seen in Figure 23b. (c) Absolute difference between T A and W K B J solutions of (b). I2 error between T A and Rayleigh soutions is 0.17. 53 contributions of all partial reflections and transmissions throughout the entire medium. For the tan 2 velocity function, v —• oo as traveltime T —> To = ir/(2y/— A B ) — C (while depth z —* oo) so that the larger To is, the more linear the function appears over a larger range of depths. The velocity function of Figure 23a blows up at To ^ 0.0075 while the velocity function in Figure 24a blows up at To ^ 0.035. Thus, the latter velocity curve approximates the linear model over a much larger range of traveltimes and depths. 3.5 Distortion of Stress Pulse In this section1 I look at the equations for stress, hi order to better understand the behavior of a stress pulse which has propagated through a transition zone. A similar analysis could be applied to a particle velocity pulse. The distortion of a stress pulse which has travelled through an inhomogeneous medium for a length of time given by the traveltime T is* A o ( a , v) — o ( a , v ) i — a(a(t — T, 0), V Q ) ' . where VQ = v(z = 0) and a(t, T) is defined in equation (16). ACT is given explicitly (from (31a) and (34)) as A a ( a , v) = f - 1 b(a) + 2 A U - J hx (a) (46) where 6(a) is the applied stress pulse and hi(a) is given by equation (33). The first term in equation (46) only effects the amplitude of the stress pulse. As expected, at a traveltime T, the amplitude distortion is proportional to v1//2 (T) /1 ; ^ 2 — 1 (Gray, 1986). This effect is easily seen in Figure 10b where (recalling that the applied 54 1 / 2 stress pulse has been scaled by v0 ) the peak amplitude of the stress pulse has increased from \/3 at T = 0 to \f\ at T = 1.0 s. Phase distortion is controlled by the second term of equation (46). The distortion arises because h i ( a ) is not symmetric wi th respect to the peak of the applied stress pulse (i.e. Ricker wavelet). I now consider the effect of increasing the velocity gradient. This is equivalent to increasing the constants A and B (which occur in the velocity functions (30)), and 1 / 2 hence, the constant 7 = 2A/v0 . Firs t I write A c as A C T = Da + Dp where fv1/2 Da = I . — 1 j b(a) ~ amplitude distortion term, w / v'l2\ Dp = 2A I 1 - T ^ ' h i ( a ) = phase distortion term. V vo J 1 / 2 As the gradient increases (i.e. as 7 —> 00), /ii(ft) —> b ( a ) / ( - y v 0 ) so the phase distort ion term, / v^2\ b{a) Thus , in this case the amplitude distortion cancels the phase distortion. So, in the l imit , as the velocity function becomes essentially a velocity discontinuity at z = 0, Aa —• 0 so that the stress pulse is not distorted at a l l , i.e., the velocity discontinuity is "invisible". 1 / 2 Thus , the wavelet, which starts wi th a peak amplitude of vQ , is not distorted at all after travelling through the inhomogeneity. This result is clearly seen in Figures 12b and 15b. For extremely small gradients A , B —* 0 so that Dp —> 0. Thus the distortion in the stress pulse is due solely to the amplitude distortion term Da. The peak amplitude 55 of the wavelet at traveltime T is then given by (i;1/2(r)/«J/2)6(a). Th i s result can be seen in the low velocity gradient examples (Figures 10b and 13b). It is clear then that there is very little phase distortion for either very high or very low velocity gradients. I find that the wavelet undergoes max imum phase distortion when the length of the transition zone is of the order of 1 / / m a x where, again, / m a x is the peak frequency component of the applied stress pulse. This can be seen in Figures l i b and 14b where the phase distortion (though slight) is plainly visible. 3.6 Reflection a n d Transmission Coefficients In this section I calculate reflection and transmission coefficients as a function of time for the hyperbolic tangent squared velocity model. Results for the hyperbolic cotangent squared velocity function are similar. The equation used to generate reflection and transmission coefficients as functions of time and traveltime is (Bruckstein et al . , 1985) WR(t,T)~ TL(t,T) RR(t,T)' ^(*,o)" WL(t,0) RL(t,T) TR(t,T) WL(t,T) (47) where * denotes the convolution operator. In equation (47) WR and WL are the right and left propagating waves, Xf, and RL are the transmission and reflection coefficients for right travelling waves, TR and RR are the transmission and reflection coefficients for left travelling waves (see Figure 25). WR(t,T) and Wi(t,T) are simply the Riemann invariants, F(a,v) and S(a,v), writ ten as functions of traveltime and time, i.e., WR(t,T) = F(a(t,T),v(T)) WL(t,T) = S(a(t,T)MT)). 56 W R(0.t) wR(T.t) Inhomogeneous w L(o,t) Medium WL(T,t) ^ 0 T F i g u r e 25 . Situation modelled by hyperbolic tangent squared velocity function.; WR and WL are right and left travelling waves related by equation (40). For large enough traveltime WL(t,T) = 0. If To is a large enough traveltime so that the hyperbolic tangent squared velocity function has leveled off to its asymptotic value (i.e., v{T) —> A/B as T —> oo) then equations (19), (25) and1 (28) combine to give WL(t,T0) = h^t) (^A ( | y / 2 - B ( | ) 1 / 2 ^ J ! = Thus , fbr traveltimes T > TQ there wi l l be no backscattered energy. I w i l l determine TL and at such a traveltime. (Computing reflection and transmission coefficients when W^{t,T) ^ 0 is more difficult and w i l l not be considered here.) Equat ion (47) then gives WR{t,T0) =TL{t,T0)*WR{t,0) WL{t,0)=.RL{t,To)*WR{t,0) (48a) (48b) 57 so that TL and Ri can be obtained by deconvolution. WR and WL can be written expl ici t ly from equations (25), (28) and (34) as ^ ( t , r ) = ^ f » - M « ( ( , r ) ) Bv1/2(T) 1 / 2 V J and WL{t,T)=h1(a{t,T)) y / 2 ( r ) - B v1 / 2 ^ ) where h\{a) is defined in equation (33). Then W « ( . , 0 ) = ^ ^ + A i ( « ( « . 0 ) ) flu 1 / 2 0 1 / 2 V J WR(t,T0) = d ( a ( < ; / f o ) ) + 2/i1(a(t,ro)) 1 / 2 _ 1 / 2 V J (49a) (49b) and W L M ) =/n(a(i,0)) T / i ~ B V (49c) 1 / 2 where 6(a) is the applied stress pulse (Ricker wavelet of peak amplitude v0' ). Note that WR(t,0)+WL(t,0) = b{a)/vl/2, the applied source pulse scaled by v^2. For extremely low velocity gradients, as A and B (and hence 7) —> 0, hi(a) —* u 0 1/f2 J" b(x) dx so that WR(t,0) WR{t,T0) b{a{t,0)) 1 / 2 b(a(t,T0)) , 1 / 2 as A , B - > 0 , as A,B-^0 (50a) (50b) 58 and WL{t,0) l / 2 v, 0 a ( t , 0 ) 6(x) dx as A , £ —> 0. (50c) For very high velocity gradients hi(a) —> b(a)/(2A) so that W R ( * , 0 ) 5 1 / 2 , 1 6(a(*,0)) as y i , /? -> oo 1 / 2 (51a) WR{i,T0) -> ( ^ ) 6 ( a ( i ,T 0 ) ) as -> oo (51b) and WL(t,0) B 1 / 2 ,4 I / 2 fc(a(t,0)) as oo (51c> Figure 26 shows WR{t,T0), WL{t,0), WR{t,0), TL{t,T0) and RL{t,T0) for the low gradient, t a n h 2 velocity model of Figure 10a. For this ease equations (50a) and (50b) indicate that both the incident and transmitted waves should look like the applied stress 1 / 2 l / 2 pulse 6(a) scaled by v0 . Recall ing that the peak amplitude of 6(a) is v0 means that both WR(t,0) and WR(t,To) should be Ricker wavelets of peak amplitude near 1 as seen in Figures 26a and 26b. It is not surprising then to find that the transmission coefficient in this case is a time-shifted delta function. For very low velocity gradients equation (50c) indicates that the reflected wave Wi(t,0) wi l l appear as a scaled and integrated version of the applied pulse as seen in Figure 26d. The peak amplitude of W^(t,0) is several orders of magnitude smaller than the peak amplitude of WR(t,Q) indicating very little reflected energy. The reflection coefficient (Figure 26e) for this 59 Incident Wave 0.15 0.30 0.45 Time (s) Transmitted^ Wave Transmission* Coefficient 1.20 1.35 1.50 Time (s) 1.0 0^ 8 0.6 0.4 0.2 0.0 (c) J . l _ 0.8 1.0 1.2 Time (s) Reflected: Wave; Reflection Coefficient W L(t,0) (d) j i ' ' ' ' i • ' 0.15 0.30 0.45 Time (s) 0.0005 F 0.0004 -0.0003 0.0002 h 0.0001 0.0 h -0.0001 -0.0002 -0.0003 -0.2 0.0 0.2 Time (s) F i g u r e 26. (a) Incident wave WR(t,0), (b) transmitted wave Wft(t,To), (c) trans-mission coefficient T£,(£, To), (d) reflected wave Wi(t,0) and (e) reflection coefficient RL{1,TO) for velocity model shown in Figure 10. 60 case has the characteristic shape of a Hilbert transform. The small oscillations seen here are numerical artifacts. Figure 27 shows the results corresponding to the intermediate velocity gradient example seen in Figure 11a. The incident and transmitted waves are again very similar (but not exact), both being slightly distorted Ricker wavelets with peak amplitudes near 0.97. These results show that the second terms on the right-hand sides of equations (49a) and (49b) are small compared to the dominant applied stress pulse term but do cause a visible amount of distortion (Figures 27a, 27b). The transmission coefficient is again a delta function, indicating the similarity of the two wavelets. The peak amplitude of the reflected wave (Figure 27d) is about 5% of the peak amplitude of the incident wave indicating a significant amount of backscattered energy. The reflection coefficient again has the form of phase-rotated delta function. This occurs because W/,(<,0) is approximately a phase-rotated: version of Wn(t,0). Figure 28 shows the results corresponding to the high velocity gradient example seen in Figure 12a. For this case equation (51a) indicates that the incident wave Wn(t,0) should be a scaled version of the applied stress pulse, with a peak amplitude of ap-proximately [BVQ/A + l]/2 = 0.875 as seen in Figure 28a. From equation (51b) the transmitted wave should also be a scaled version of the applied pulse with peak am-plitude near {voB/A)1/2 ~ 0.866 as seen in Figure 28b. The incident and transmitted waves are again Ricker wavelets with peak amplitudes of about 0.874. The transmission coefficient is a delta function as expected. The reflected wave W^{t,Q) is also a scaled version of b(a) with peak amplitude near 0.125. In this case a large proportion of inci-dent energy has been reflected. The reflection coefficient for this case is a band-limited 61 Incident Wave 0.15 0.30 0.45 Time (s) Transmitted Wave Transmission Coefficient 0.30 0.45 0.60 Time (s) Reflected Wave 0.15 0.30 0.45 Time (s) 1.0 0.8 0.6 0.4 0.2 0.0 T L(t, (c) a L 0.0 0.2 Time (s) Reflection Coefficient 0.005 | 0.004 - R L ( t J 0 ) 1 0.003 1. 0.002 \ ( e ) 0.001 0.0 — • i i i -0.2 0.0 0.2 Time fs) Figure 27 . (a) Incident wave Wi?(t,0), (b) transmitted wave Wn(t,To), (c) trans-mission coefficient Ti(t,To), (d) reflected wave WL(t,0) and (e) reflection coefficient RL{t,To) for velocity model shown in Figure 11. 62 Incident Wave 0.15 0.30 0.45 Time (s) Transmitted Wave Transmission Coefficient 0.15 0.30 0.45 Time (s) 1.0 0.8 0.6 0.4 0.2 0.0 I T L(t,T 0) 1 (c) 1 I I . -0 .2 0.0 0.2 Time (s) Reflected Wave Reflection Coefficient 0.15 0.30 0.45 Time (s) 0.015 R L ( t J 0 ) 0.010 i 0.005 _ (e) 0.0 • ' W \ / | A /w^ 1 1 -0.2 0.0 0.2 Time Cs^ F i g u r e 28. (a) Incident wave Wp(t,0), (b) transmitted wave WR(t,To), (c) trans-mission coefficient Ti(t,To), (d) reflected wave WL{1,Q) and (e) reflection coefficient RL{1,TO) for velocity model shown in Figure 12. 63 delta function. The band-limited nature arises in the numerical deconvolution used solve for Ri(t,To). 64 C H A P T E R I V C O N C L U S I O N S The truncated asymptotic solution is an accurate t ime-domain solution of the inho-mogeneous wave equation for smoothly varying media. Solutions for stress and particle velocity are accurate to wi th in the error of the required curve fit and fare well in com-parison wi th the second order W K B J solutions. The methods validity was confirmed using a linear velocity function wi th a high gradient—a fairly difficult function to fit. Other inhomogeneities wi l l allow much better fits and and hence more accurate solu-tions. The main advantage of the T A solution is that the assumption of geometrical' optics,, that the material! stratification length is much greater than the wavelength, is not required. The T A solution includes a l l accumulated partial reflections and trans-missions and, for the case N = 1, is exact for three classes of inhomogeneities given by equations (30). Another advantage of the method is that both p and E are allowed to vary smoothly wi th depth. The assumption of constant density is not required. The principal disadvantage of the method is that it is exact only for those three classes. However, reasonable results for other monotonic, smoothly varying impedance functions may be obtained using curve-fitting. As well , taking N > 1 w i l l increase the number of classes of inhomogeneities for exact solutions to exist. Increasing the number of classes w i l l also improve the accuracy of any required curve fits. Results have been presented for media wi th smoothly varying density and velocity. T w o solutions were calculated using the same impedance as a function of traveltime. The first solution was calculated for a medium in which the principal variation in impedance was due to velocity. In the second, the principal variation was mainly due to density. The 65 amplitudes of the seismograms calculated for each case were similar but the traveltime curves were widely separated. A n investigation into the validity of the second order W K B J solution lead to the following conclusions: For two of the three allowed forms of velocity ( tanh 2 and co th 2 ) the W K B J solution is an accurate approximation to the exact (TA) solution when the length of the transition zone is of the order of, or larger than, the length of the applied wavelet. For higher gradient transitions the W K B J solution fails. For the thi rd allowed form of velocity ( tan 2 ) the W K B J solution was accurate even for very steep gradients. This case, however, must be considered separately from the previous two because the t a n 2 function does not represent a true transition zone. The entire length of the t a n 2 function must be considered when comparing results for this case. These results agree wi th equations (42) and (43) which govern the validity of the W K B J solution. The condition (43) was found to be stricter than necessary for an accurate approximation. In part icular, only the highest order term of this equation need be satisfied for the approximation to be good. For a stress pulse propagating through a t anh 2 t ransit ion zone phase distortion is apparent only when the length of the transition zone is of the order of the length of the wavelet. Even in this case the phase distortion is slight. For longer transition zones the wavelet is distorted in amplitude only. For extremely short transition zones, as the velocity function becomes a discontinuity at z = 0, the wavelet is unaffected. Reflection and transmission coefficients as a function of time have been calculated for the t a n h 2 transition zone. Three cases were studied—low, intermediate and high velocity gradients. In all three cases the transmission coefficient was a time-shifted delta function. The reflection coefficient for the low gradient case has the shape of a Hilbert 66 transform. For steeper gradients the reflection coefficient appears as a phase-rotated delta function. For extremely steep velocity gradients the reflection coefficient is a delta function. The amount of reflected energy is largest for the steep gradient case and is controlled by the magnitude of the difference between vo and the asymptotic value of the velocity. The truncated asymptotic method presented here has been used to solve the l - D inhomogeneous wave equation. Possible extensions include applying the method to two and three-dimensional wave propagation and solving other equations using the truncated asymptotic approach. 67 R E F E R E N C E S Bender, C M . , and Orszag, S .A. , 1978, Advanced mathematical methods for scientists and engineers: M c G r a w - H i l l Book Co . Ber ryman , L . H . , Goupi l laud, P . L . , and Waters, K . H . , 1958, Reflections from multiple transit ion layers. Part I. Theoretical results: Geophysics, 23 , 223-243. Bortfeld, R . , 1960, Seismic waves in transition layers: Geophys. Prosp., 8, 178-217. Brekhovskikh , L . M . , 1960, Waves in layered media: Academic Press Inc. Die t r ich , M . , and Bouchon, M . , 1985, Synthetic vertical seismic profiles in elastic media: Geophysics, 50, 931-949. El las , G.I'., 1931, Das Verhalten elektromagnetischen Wellen bei raurnlich veranderlichen eleetrisehen Eigenschaften. Elektr . Nach. Tecknik. 8, 4. Epste in , P., 1930, Reflection of waves in an inhomogeneous absorbing medium. Proc. Nat . A c a d . Sci. U . S . A . 16, 627-637. Foster, M . , 1975, Transmission effects in the continuous one-dimensional seismic model: Geophys. J . Roy. As t r . S o c , 42 , 519-527. Ganley, D . C . , 1981, A method for calculating synthetic seismograms which include the effects of absorption and dispersion: Geophysics, 46 , 1100-1107. Gans , R . , 1915, Fortpflanzung des Lichtes durch ein inhomogenes M e d i u m . A n n . Physik . 47 , 709. G u p t a , R . N . , 1965, Reflection of plane waves from a linear transition layer in l iquid media: Geophysics, 30, 122-132. 68 Gray, S .H . , 1984, A problem of discrete approximations to an arbitrari ly varying one-dimensional seismic model: Geophys. J . Roy. As t r . S o c , 78, 431-437. Heller, G . S., 1953, Reflection of acoustic waves from an inhomogeneous fluid medium. J . Acoust . Soc. A m . 25, 1104. M u l l i n , C . J . , 1953, Solution of the wave equation near an extremum of the potential. Phys. Rev. 92 , 1323. P a i , D . M . , 1985, A new solution method for the wave equation in inhomogeneous media: Geophysics, 50, 1541-1547. Peterson, R . A . , Fi l l ippone, W . R . , and Coker, F . B . , 1955, The synthesis of seismograms from well log data: Geophysics, 20, 516-538. Rayleigh, L o r d , 1894, Theory of Sound, volume I: M a c M i l l a n P u b l . C o . Ry tov , S . M . , and Yudkevich, F . S . , 1946, Electromagnetic wave reflection from a layer wi th a negative dielectric constant. J . E x p t l . Theoret. Phys. (U.S.S.R.) 10 , 285. Scholte, J . G . J . , 1961, Propagation of waves in inhomogeneous media: Geophys. Prosp., 9, 86-115. Seymour, B . R . , and Varley, E . , 1970, High frequency periodic disturbances in dissipative systems. I. Small amplitude finite rate theory: Proc. Roy. Soc. London, A 3 1 4 , 387-415. Seymour, B . R . , and Varley, E . , 1984, Exact representation for acoustical waves in strat-ified media: Institute of Appl ied Mathematics , Un iv . B r . Columbia , Tech. Report 84-11. Seymour, B . R . , and Varley, E . , 1987, Exact representation for acoustical waves when the sound speed varies in space and time: Studies A p p l . M a t h . , 76 (to appear). 69 Sheriff, R . E . , and Geldart , L . P . , 1983, Explorat ion seismology, 2: Cambridge Univ . Press. Varley, E . , Venkataraman, R . , and Cumberbatch, E . , 1971, The propagation of large ampli tude tsunamis across a basin of changing depth. Part 1. Offshore behavior: J . F l u i d . Mech . , 49 , 775-801. Wal lo t , J . , 1919, Der senkrechte Durchgahg elektromagnetiseher Wellen dureh eine Schicht raumlich veranderlicher Dieleektrizitatsknostante. A n n . Physik 60 , 734. Wuenshcel, P . C . , 1960, Seismogram synthesis including multiples and transmission co-efficients: Geophysics, 25, 106-1291 70 A P P E N D I X A Derivation of Equations (20 ) In this appendix I rewrite equations (l) and (2) in terms of F and 5, the Riemann invariants defined in equation (19), to obtain equations (20). This requires changing variables from {z,t) to (J, a) in equations (l) and (2), which I rewrite here as oz = put 1 uz = —o E Now OL = |(< — T(z)) so that t can be considered a function of z and a. Thus, if £(z,a) = o{z,t(z,a)) U(z,a) = u(z,t{z,a)) then, by the chain rule (A - 1) {A - 2) T,z = oz + attz and Uz = uz + uttz. (A - 3) Substituting equation (A-3) into equations (A-l) and (A-2), and recalling that tz — l/v, gives rise to E z = put + — and Uz = -^at + — . (A - 4) t; E v 71 Multiplying through by I — pv in the second equation of (A-4) (remembering that v = y/E/p) gives IUZ = — + put v which is the right-hand side of the first equation of (A-4). Thus E z = IUZ. {A - 5) Using the chain rule once again gives E z = E j / Z and Uz — Ujlz so that equation (A-5); becomes Ej, = IVi. (A - 6) In terms of T, and U, the Riemann invariants are F(atI} = -rl^{i:-IU} and 5(a,7) = ^/ _ 1 / ' 2 ( E + IU). 2 2 Solving for E and U and differentiating with respect to J gives E = J 1 / 2 ( S + F) and E/ = J 1 / 2 (S/ + Fj) + \ r l ' 2 { S + F), {A - 7) U — 7 _ 1/ 2(5 - F) and Ut = T1'2 {Sr-Fj)-\rz'2{S - F). (A - 8) Substituting equations (A-7) and (A-8) into equation (A-6) gives S + 11Fi = 0 {A-9) 72 which is equation (20a). To obtain equation (20b) I start with the second equation in (A-4) and recall that Uz = Ujlz so that EIzUj = at + Iut = — {a + Iu). But JT = | f #- = \4~ so that at at oa I aa — - (E + IU) = 2ElzUt. {A- 10) Oft From the first equations of (A-7) and (A-8) Y,a = 7 1 / 2 {Sa Fa) and Ua r 1 / 2 ( S „ - Fa). These, together with the expression for Uj in equation (A-8) are substituted into equa-tion (A-10) to get I1/2(Sa + Fa) + I1'2 (Sa - Fa) = 2EIZ (j-1/2 (Sr -Fj)- \ir^{S - F}j .. (A-11). S is given in equation (A-9) as —2IFj so that equation (A-ll) reduces to I z ( „ 1 Sa = E j U S J + -F 73 But Iz = ITTZ = J / r and, since E/v = I If 7 Z is normalised with respect to a small parameter, e, this last equation becomes Sa = e*f (lST + l-F which is equation (20b). 74 A P P E N D I X B Derivation of the Rayleigh Solution for Stress In this appendix I derive the Rayleigh solution for stress o. I take p to be constant for s implici ty so that I can replace I by v. The wave equation for stress, which follows from equation (3) is dz2 E{z) dt2 v2{z)dt2' 1 J I want to solve this equation for velocities which vary linearly wi th depth. I take v(z) - cz -I- d where c and d are constants and assume a solution of the form a{z,t) = E{z,u)eiujt {B-2) for time-periodic inputs. Substituting equation (B-2) into equation ( B - l ) leads to the ordinary differential equation d 2 S w 2 , ' + 7 - ^ £ = 0, (B-S) dz2 (cz + d)' a case first studied by Rayleigh (1894, §148) for periodic boundary data. To solve equation (B-3) I set f the characteristic exponent gives 75 — cz + d and try as a solution S = $ k . Solving for so that Z{z,u) = P(u)t*±m (B-4) where m = —; {B ~ 5} For forward travelling progressive waves I take the minus sign (Rayleigh, 1894, §148) in the exponent of equation (B-4). The total frequency-domain solution is obtained by mult iplying the right-hand side of equation (B-4) by B(u>),. the Fourier transform of the applied stress pulse b(t). Thus the solution is now E(2 ,w) = P{u)B{uj)$*-im. (B - 6) To solve for the constant P(w) I demand that this solution match the truncated asymp-totic solution exactly at z = 0. A t z = 0 the T A time-domain solution is given by equation (15), i.e., <7 T A (0,0 = b(t). So the truncated asymptotic frequency-domain solution at z ~ 0 is ETA{0,u;) = B{u). 76 Equa t ing this solution with E(0 , UJ) from equation (B-6) gives P(w) = d~l(2+im. So the final frequency-domain Rayleigh solution is E(z ,w) = B(u>) • , \ 1 / 2 exp ( — iml 'n cz + d d (B-7) and the time-domain solution is a(z,t) = F - l 1 / 2 R , v « f£ + rf\ / • r (ez + d B(u>) | ——— J exp I; irn In - 8)! the inverse Fourier transform of equation (B-7). The Rayleigh solution for particle velocity u is obtained in a similar manner. The wave equation for particle velocity is d2u _ p d2u I dE du ~dz^ ~ E[z)'di2~ + E(z)~dzd~z 1 d2u 1 dEdu + v2{z)dt2 E(z) dz dz (B - 9) and its t ime-domain solution is found to be u{z,t) = F - l 1/2 n( \ i'cz + d \ ( fcz + d -B(w) [ — - — ) exp( —tnm where n — — for forward travelling progressive waves. 78 A P P E N D I X C Derivation of the First Four Terms in the W K B J Solution for Stress In this appendix I derive the first four terms in the W K B J solution for stress o. I again take p = constant for simplicity and replace / by v. The wave equation for stress, which follows from equation (3) is <7,, = v2(z) ( C - l ) In order to solve equation ( C - l ) I assume a solution of the form Then z) = c(u,0) {Fx(z) + •• £ & ) e *"'«.<«>. (C - 2) F\zz + — ^ + - iuFizF0z - iuFxF0zz - F2zF0l ILJ -F2FQZZ FZZFQZ FZF0ZZ. no iu> iuF\zFoz — F2ZF0. + [iuj^FiFol + iuF2F0l + F3F0 ( C - 3 ) Equations (C-2) and (C-3) are substituted into equation ( C - l ) and coefficients of powers of UJ are matched. The coefficients of the u2 term give rise to the ordinary differential equation 79 Choosing the positive sign for waves travelling to the right gives i.e., the term i*o is just the traveltime T(z). Coefficients of the u1 term result in the following ordinary differential equation F0zF2 - 2F0zFiz - F0zz = — (C - 5) where equation (C-4) gives FQZ — l/v, FQZZ — —vz/v2 so that equation (C-5) reduces to 2Flz - —Fj = 0. v Separating variables and integrating between 0 and z gives vl!2iz\ Coefficients of the u;0 term give F2F0zz + 2F0zF2z = Flzz {C - 7) with F\zz found from equation (C-6). Substituting for F0zz and F0z into equation (C-7) 80 and mul t ip ly ing both sides by vxl2(z) reduces equation (C-7) to Thus , z z z' z F2(z) = \v"\z) lZv"\x)Flzz{x)dx. ( C - 8 ) * Jo Coefficients of the u ; - 1 term give rise to the following ordinary differential equation: FiFozz + '^F0zF3z = F2zz which is identical to equation (C-7). The solution for Fs is then analogous to the solution for F 2 ; F3(z) = ^(z) j\^(x)F2zz{x) dx. (C - 9) w i t h F2zz found from equation (C-8). The second order W K B J solution for a as a function of depth z is comprised of the first three terms F0, F'i and F2 given in equations (C-4), (C-6) and (C-8). In Section 3.3, where I compare the T A solution to the W K B J solution, I am using velocity functions parameterised by traveltime T. Changing variables from z to T gives the W K B J terms as functions of traveltime. These are found to be FQ(T) = T, F1(T)= ' , 1 fT F2{T) = -v1'2 / (v~1/2FlTT - v-3/2vTF1T) dT ( C - 10) 2 Jo 81 and MT) = U 1 ' 2 I {v~1/2F2TT - v-3'2vTF2T) dT. [C - 11) 2 Jo For velocity functions of the form given in equation (30), equation (C-10) reduces to F2(T) = AFl(T)(v-l'2{T)-v-112) while equation (C-ll) becomes F3(T) = -Av-1/2F2(T} which is equation (45). The second order WKBJ solution for u is obtained in a similar manner to the solution for o described above. 82 A P P E N D I X D Derivation of the Coefficients in the Exponential Series Form of the W K B J Solution for Stress In this appendix I relate the W K B J coefficients Fo,...,Fs of equation (C-2) (see Appendix C) to the coefficients S o , . . . , S3 of the exponential series o3(T) (see equation (40)). I do this by equating (FX{T) + SF2(T) + 62F3(T]\eF^/6 = exp J £ ^ ( T ) (D-l) w here 6 = —1/iuj. The right-hand side of equation ( D - l ) can be writ ten as exp J exp (S i + 6S2 + 62S3) (D - 2) to give S o /S = FQ/S from the left hand side of equation ( D - l ) or S 0 ( T ) - F0{T) = T. The second term in expression (D-2) is exp (S1+6S2 + 62S3) = exp Si exp(tSS 2 + S2S3) = exp S x ( l + (SS2 + S2S3) + (6S2 + 62S3)2/2 + = exp S i ^1 + 6S2 + 62(S3 + ^ S 2 ) + O ^ 3 ) ^ 83 where the right-hand side is obtained by expanding the exponential and collecting terms in 6. Equating this last expression with Fx + SF2 + 62F3 and ignoring terms of order 63 and higher gives F2 F3 1/F2\2 S1 = lnF0l S2 = - and S a = - - - { - ) which are equations (44).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Truncated asymptotic solution of the one-dimensional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Truncated asymptotic solution of the one-dimensional inhomogeneous wave equation Zelt, Barry Curtis 1987
pdf
Page Metadata
Item Metadata
Title | Truncated asymptotic solution of the one-dimensional inhomogeneous wave equation |
Creator |
Zelt, Barry Curtis |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | I present a new time-domain method for solving for the stress and particle velocity of normally incident plane waves propagating in a smoothly varying one-dimensional medium. Both the Young's modulus E and the density ⍴ are allowed to vary smoothly with depth. The restriction of geometrical optics, that the wavelength be much less than the material stratification length, is not required in this new method. The infinite geometrical optics expansion is truncated after n terms, imposing a condition on the acoustic impedance I for exact solutions to exist. For the case ռ = 2 there are three general classes of impedance functions for which the resultant expansion is uniform and exact. To check the numerical validity of the "truncated asymptotic" (TA) solution results are calculated for the case of a medium with a linear velocity gradient for which there is an exact solution in the frequency domain. Since a linear velocity gradient is not one of the foregoing classes of impedance functions, a curve-fitting approach is necessary. The TA method compares favourably to the exact solution and is accurate to within the error of the required curve fit. Two classes of synthetic seismograms are calculated for smooth velocity and density variations. The same impedance as a function of traveltime is used for both classes. In the first class the principal variation in impedance is due to velocity, while in the second it is mainly due to density. The amplitudes in both classes of synthetic seismograms are very similar, but, as expected, the traveltime curves for each class are widely separated. For the case ⍴ = constant the TA solution is used as a bench-mark to analyse a two-term WKBJ approximation for three classes of velocity functions. The velocity functions are such that the TA solutions are exact. For two of the classes the WKBJ solution performs well when the length of the transition zone is of the same order, or larger, than the length of the applied wavelet. For steeper velocity gradients the WKBJ solution begins to differ significantly from the exact TA solution. The WKBJ solution for the third class performs extremely well even for steep gradients. Equations governing the validity of the WKBJ solution are examined to explain the above results. Equations are derived to describe the distortion of a stress pulse propagating through a transition zone. For small velocity gradients (relative to the length of the applied pulse) the wavelet changes amplitude but its phase is not effected. As the gradient increases and the velocity function becomes a discontinuity at z = 0 the wavelet travels through undistorted. Only when the transition zone width is of the order of the length of the wavelet is there any visible phase distortion. Reflection and transmission coefficients as functions of time are calculated for low, intermediate and high gradient transition zones. The transmission coefficient is a delta function in each case. The reflection coefficient has the shape of a Hilbert transform for low gradients. For higher gradients the reflection coefficient approaches the shape of a delta function. |
Subject |
Wave equation |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052971 |
URI | http://hdl.handle.net/2429/26677 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A6_7 Z44.pdf [ 4.12MB ]
- Metadata
- JSON: 831-1.0052971.json
- JSON-LD: 831-1.0052971-ld.json
- RDF/XML (Pretty): 831-1.0052971-rdf.xml
- RDF/JSON: 831-1.0052971-rdf.json
- Turtle: 831-1.0052971-turtle.txt
- N-Triples: 831-1.0052971-rdf-ntriples.txt
- Original Record: 831-1.0052971-source.json
- Full Text
- 831-1.0052971-fulltext.txt
- Citation
- 831-1.0052971.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052971/manifest