Teleseismic waveform modelling with a one-way wave equation

It is now widely accepted that elastic properties of the continental lithosphere and the underlying sublithospheric mantle are both anisotropic and laterally heterogeneous at a range of scales. To fully exploit modern three-component broad-band array data sets requires the use of comprehensive modelling tools. In this work, we investigate the use of a wide-angle, one-way wave equation to model variations in teleseismic 3-D waveforms due to 2-D elastic heterogeneity and anisotropy. The one-way operators are derived based on a high-frequency approximation of the square-root operator and include the effects of wave propagation as well as multiple scattering. Computational cost is reduced through a number of physically motivated approximations. We present synthetic results from simple 1-D (layer over a half-space) and 2-D (subduction zone) models that are compared with reference solutions. The algorithm is then used to model data from an array of broad-band seismograph stations deployed in northwestern Canada as part of the IRIS-PASSCAL/LITHOPROBE CANOE experiment. In this region radial-component receiver functions show a clear continental Moho and the presence of crustal material dipping into the mantle at the suture of two Palaeo-Proterozoic terranes. The geometry of the suture is better defined on the transverse component where subduction is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. The modelling reproduces the main features of the receiver functions, including the effects of anisotropy, heterogeneity and finite-frequency scattering.     

Numerical scheme for the inversion of acoustical impedance profile based on the Gelfand-Levitan method

Summary. An exact method for the solution of the inverse problem in plane wave propagation modelled after the Gelfand-Levitan technique is reviewed and refined. A numerical scheme for the solution of the integral equation that arises in the method is proposed. A discussion on the stability and an error analysis of the numerical approximation are presented. The applicability of the inversion algorithm is demonstrated in a numerical experiment.     

Surface wave tomography: global membrane waves and adjoint methods

We implement the wave equation on a spherical membrane, with a finite-difference algorithm that accounts for finite-frequency effects in the smooth-Earth approximation, and use the resulting 'membrane waves' as an analogue for surface wave propagation in the Earth. In this formulation, we derive fully numerical 2-D sensitivity kernels for phase anomaly measurements, and employ them in a preliminary tomographic application. To speed up the computation of kernels, so that it is practical to formulate the inverse problem also with respect to a laterally heterogeneous starting model, we calculate them via the adjoint method, based on backpropagation, and parallelize our software on a Linux cluster. Our method is a step forward from ray theory, as it surpasses the inherent infinite-frequency approximation. It differs from analytical Born theory in that it does not involve a far-field approximation, and accounts, in principle, for non-linear effects like multiple scattering and wave front healing. It is much cheaper than the more accurate, fully 3-D numerical solution of the Earth's equations of motion, which has not yet been applied to large-scale tomography. Our tomographic results and trade-off analysis are compatible with those found in the ray- and analytical-Born-theory approaches.     

Surface-wave propagation in a cracked poroelastic half-space lying under a uniform layer of fluid

The effect of cracks on the elastic properties of an isotropic elastic solid is studied when the cracks are saturated with a soft fluid. A polynomial equation in effective Poisson's ratio is obtained, whose coefficients are functions of Poisson's ratio of the uncracked solid, crack density and saturating fluid parameter. Elastic and dynamical constants used in Blot's theory of wave propagation in poroelastic solids are modified for the introduction of cracks. The effects of cracks on the velocities of three types of waves are observed numerically. The frequency equation is derived for the propagation of Rayleigh-type surface waves in a saturated poroelastic half-space lying under a uniform layer of liquid. Dispersion curves for a particular model of oceanic crust containing cracks are plotted. The effects of variations in crack density and saturation on the phase and group velocity are also analysed.     

Wave propagation characteristics in frozen saturated soil

Ultrasonic detection technology is of great significance in the detection and evaluation of physical and mechanical properties of frozen soil, but wave propagation characteristics in frozen soil are unclear. Based on the three-phase composition of frozen saturated soil and the mixture theory, considering Bishop's effective stress formula, the wave propagation equations are establish for frozen saturated soil. In wave propagation, an entropy inequality was introduced to describe the coupling of different phases. The analytic expressions of propagation velocity and attenuation law of waves in frozen soil are obtained, and wave propagation characteristics in frozen saturated soil are discussed. Results show that four types of waves(i.e., P1, P2, P3 and S) are found in frozen saturated soil and all four wave types are dissipative waves, in which the attenuation of P3 is the maximum. The velocity of four waves increases sharply at the excitation frequency range of 10~3–10~9 Hz,but the wave velocity at high-frequency and low-frequency is almost constant. When volume ice content increases, the wave propagation velocity of P1 and S decreases dramatically, and the velocity of P2 increases gradually, but P3 velocity increases first and then decreases to zero with increasing saturation. The attenuation coefficients of P1 and S waves begins to increase gradually when the volume ice content is about 0.4, P2 increases first and then decreases with an increase of volume ice content and P3 increases with the volume ice content and decreases rapidly from extreme to zero.     

Split-step complex Padé-Fourier depth migration

We present a split-step complex Padé-Fourier migration method based on the one-way wave equation. The downward-continuation operator is split into two downward-continuation operators: one operator is a phase-shift operator and the other operator is a finite-difference operator. A complex treatment of the propagation operator is applied to mitigate inaccuracies and instabilities due to evanescent waves. It produces high-quality images of complex structures with fewer numerical artefacts than those obtained using a real approximation of a square-root operator in the one-way wave equation. Tests on zero-offset data from the SEG/EAGE salt data show that the method improves the image quality at the cost of an additional 10 per cent computational time compared to the conventional Fourier finite-difference method.     

Quasi-isotropic approximation of ray theory for anisotropic media

Wave propagation in weakly anisotropic inhomogeneous media is studied by the quasi-isotropic approximation of ray theory. The approach is based on the ray-tracing and dynamic ray-tracing differential equations for an isotropic background medium. In addition, it requires the integration of a system of two complex coupled differential equations along the isotropic ray.
The interference of the qS waves is described by traveltime and polarization corrections of interacting isotropic S waves. For qP waves the approach leads to a correction of the traveltime of the P wave in the isotropic background medium.
Seismograms and particle-motion diagrams obtained from numerical computations are presented for models with different strengths of anisotropy.
The equivalence of the quasi-isotropic approximation and the quasi-shear-wave coupling theory is demonstrated. The quasi-isotropic approximation allows for a consideration of the limit from weak anisotropy to isotropy, especially in the case of qS waves, where the usual ray theory for anisotropic media fails.     

Interaction of a compressional impulse with a slot normal to the surface of an elastic half space

Summary. An improved finite difference scheme is applied to simulate wave propagation in the vicinity of a slot normal to the surface of an elastic half space. It provides visualization of the scattered wave pattern at a sequence of time steps, and also the components of displacement as functions of time at a series of observation points.
After being hit by a normally incident plane P pulse, the slot oscillates with two main cycles and two shear-compressional pairs of diffracted waves, and also Rayleigh pulses, are scattered from it. The resulting wavefronts are parallel to the vertical surfaces of the slot and curve in semicircular arcs around the bottom of the slot.
Experimental tests of the theory were performed, using 0.5–6 MHz ultrasonic pulses on duralumin cylinders with surface-breaking slots ranging from 0.5–2 mm in width and from 2–6 mm in depth. The numerical results were confirmed by these experiments.     

Transient and impulse responses of a one-dimensional linearly attenuating medium — I. Analytical results

Summary. The transient and impulse responses (Green's function) for onedimensional wave propagation in a standard linear solid are calculated using a Laplace Transform method. The spectrum of relaxation times is chosen so as to model a constant Q medium within an absorption band covering a broad frequency range which may be chosen so as to include the seismic frequencies. The inverse transform may be evaluated asymptotically in the limit of very long propagation times using the saddle point method. For shorter propagation times the method of steepest descent may be modified so as to yield an accurate first motion approximation. The character of the small amplitude precursor to the large amplitude Visible' signal is investigated analytically. It is shown that the signal velocity is intermediate between the high-frequency ('unrelaxed') and the low-frequency ('relaxed') limits of the phase velocity.     

Propagation of harmonic plane waves in a general anisotropic porous solid

Wave propagation is studied in a general anisotropic poroelastic solid. The presence of dissipation due to fluid-viscosity as well as hydraulic anisotropy of pore permeability are also considered. Biot's theory is used to derive a system of modified Christoffel equations for the propagation of plane harmonic waves in porous media. A non-trivial solution of this system is ensured by a determinantal equation. This equation is separated into two different polynomial equations. One is the quartic equation whose roots represent the complex velocities of four attenuating waves in the medium. The other is a eighth-degree polynomial whose roots represent the vertical slowness values for the four waves propagating upward and downward in a finite porous medium. Procedure is explained to associate the numerically obtained roots with the waves propagating in the medium. The slowness surfaces of waves reflected at the boundary of the medium are computed for a realistic numerical model. The behaviours of phase velocity surfaces are analysed with the help of numerical examples.     

Wavefield smoothing and the effect of rough velocity perturbations on arrival times and amplitudes

Geometric ray theory is an extremely efficient tool for modelling wave propagation through heterogeneous media. Its use is, however, only justified when the inhomogeneity satisfies certain smoothness criteria. These criteria are often not satisfied, for example in wave propagation through turbulent media. In this paper, the effect of velocity perturbations on the phase and amplitude of transient wavefields is investigated for the situation that the velocity perturbation is not necessarily smooth enough to justify the use of ray theory. It is shown that the phase and amplitude perturbations of transient arrivals can to first order be written as weighted averages of the velocity perturbation over the first Fresnel zone. The resulting averaging integrals are derived for a homogeneous reference medium as well as for inhomogeneous reference media where the equations of dynamic ray tracing need to be invoked. The use of the averaging integrals is illustrated with a numerical example. This example also shows that the derived averaging integrals form a useful starting point for further approximations. The fact that the delay time due to the velocity perturbation can be expressed as a weighted average over the first Fresnel zone explains the success of tomographic inversions schemes that are based on ray theory in situations where ray theory is strictly not justified; in that situation one merely collapses the true sensitivity function over the first Fresnel zone to a line integral along a geometric ray.     

An introduction to Maslov's asymptotic method

Summary. Familiar concepts such as asymptotic ray theory and geometrical spreading are now recognized as an asymptotic form of a more general asymptotic solution to the non-separable wave equation. In seismology, the name Maslov asymptotic theory has been attached to this solution. In its simplest form, it may be thought of as a justification of disc-ray theory and it can be reduced to the WKBJ seismogram. It is a uniformly valid asymptotic solution, though. The method involves properties of the wavefronts and ray paths of the wave equation which have been established for over a century. The integral operators which build on these properties have been investigated only comparatively recently. These operators are introduced very simply by appealing to the asymptotic Fourier transform of Ziolkowski & Deschamps. This leads quite naturally to the result that phase functions in different domains of the spatial Fourier transform are related by a Legendre transformation. The amplitude transformation can also be inferred by this method. Liouville's theorem (the incompressibility of a phase space of position and slowness) ensures that it is always possible to obtain a uniformly asymptotic solution. This theorem can be derived by methods familiar to seismologists and which do not rely on the traditional formalism of classical mechanics. It can also be derived from the sympletic property of the equations of geometrical spreading and canonical transformations in general. The symplectic property plays a central role in the theory of high-frequency beams in inhomogeneous media.     

Variationsl solution of long-period oscillations of the Earth

Summary The linearized equation of motion for the slightly elliptical rotating earth is obtained and using Phinney & Burridge's generalized spherical harmonics, the variational principle is derived for the normal mode oscillations of the Earth. The numerical solutions of two earth models 1066B and B1S6 are searched by minimizing the energy functional for the terrestrial spectral range longer than the lowest order free oscillation. The periods of core modes computed for the earth model B1S6, with stably stratified outer core, ranges from about 4 to 13hr and the periods for the 1066B are much more spread without clustering around the periods of 6 and 12 hr as in B1S6. The results for the earth model 1066B indicate that an outer core can support long-period oscillations even when it is not stably stratified. The Chandler wobble periods obtained are 402.3 day for B1S6 and 402.7 day for 1066B.     

Traveltimes for infrasonic waves propagating in a stratified atmosphere

The tau– p method of Buland & Chapman (1983) is reformulated for sound waves propagating in a stratified atmosphere under the influence of a height-dependent wind velocity profile. For a given launch angle along a specified azimuth, the ray parameter is redefined to include the influence of the horizontal wind component along the direction of wave propagation. Under the assumption of negligible horizontal wind shear, the horizontal wind component transverse to the ray propagation does not affect the direction of the wave normal, but displaces the reference frame of the moving wavefront, thus altering the observed incidence azimuth. Expressions are derived for the time, horizontal range, and transverse range of the arriving waves as a function of ray parameter. Algorithms for the location of infrasonic wave sources using the modified tau– p formulation in conjunction with regional atmospheric wind and temperature data are discussed.     

High-frequency radiation from crack (stress drop) models of earthquake faulting

We present a theory for the radiation of high-frequency waves by earthquake faults. We model the fault as a planar region in which the stress drops to the kinematic friction during slip. This model is entirely equivalent to a shear crack. For two-dimensional fault models we show that the high frequencies originate from the stress and slip velocity concentrations in the vicinity of the fault's edges. These stress concentrations radiate when the crack expands with accelerated motion. The most efficient generation of high-frequency waves occurs when the rupture velocity changes abruptly. In this case, the displacement spectrum has an ω^-2 behaviour at high frequencies. The excitation is proportional to the intensity of the stress concentration near the crack tips and to the change in the focusing factor due to rupture velocity. We extend these two-dimensional results to more general three-dimensional fault models in the case when the rupture velocity changes simultaneously on the rupture front. Results are similar to those described for two-dimensional faults. We apply the theory to the case of a circular fault that grows at constant velocity and stops suddenly. The present theory is in excellent agreement with a numerical solution of the same problem.
Our results provide upper bounds to the high-frequency radiation from more realistic models in which rupture velocity does not change suddenly. The ω^-2 is the minimum possible decay at high frequencies for any crack model of the source.     

The Acoustics of Pipe Arrays

Summary A theoretical analysis is given for the acoustical behaviour of the pipe-microbarograph systems used to detect acoustic gravity waves and other modes of infrasound. It is shown how to compute the response of the microbarograph to a fluctuating pressure at any one inlet port of the pipe and how the results of such computations may be used to calculate the response to a plane sound wave traversing the system.
The analysis is illustrated by numerical examples obtained by means of a computer program. These examples confirm that the tapered tube modelled after Daniels' line microphone has very good characteristics, but that good results may also be obtained using pipes of uniform bore. The work leans heavily on Benade's calculations of sound propagation in a circular conduit.     

Second-Order Correction to the Reed-Otterman Theory

Summary The theory of Reed and Otterman, which describes the propagation of a weak shock wave in an inhomogeneous medium, is extended to include the next higher-order term of relative overpressure. This extension allows a matching of theory to initial starting conditions of the shock at a point closer to the source than the Reed-Otterman theory itself. The modified theory gives a value of the shock pressure that agrees well with the numerical calculation performed by Greene and Whitaker for the upward moving shock wave from a low altitude detonation.