Kirchhoff–Helmholtz reflection seismograms in a laterally inhomogeneous multi-layered elastic medium – II. Computations

Summary. High-frequency reflection and refraction seismograms for laterally variable multi-layered elastic media are computed by using the frequency domain elastic Kirchhoff–Helmholtz (KH) theory of Frazer and Sen. Both source and receiver wavefields are expanded in series of generalized rays and then elastic (KH) theory is applied to determine the coupling between each source ray and each receiver ray at each interface. The motion at the receiver is given as a series of integrals, one for each generalized ray. We use geometrical optics and plane wave reflection and transmission coefficients for rapid evaluation of the integrand. When the source or the receiver ray field has caustics on the surface of integration geometrical ray theory breaks down and this gives rise to singularities in the KH integrand. We repair this using methods suggested by Frazer and Sen.
Examples of reflection seismograms for 2-D structures computed by elastic KH theory are shown. Those for a vertical fault scarp structure are compared with the seismograms obtained by physical modelling. Then OBS data obtained from the mid-America trench offshore Guatemala area are analysed by computing KH synthetics for a velocity model that has been proposed for that area. Our analysis indicates the existence of a small low-velocity zone off the trench axis.
No head wave arrivals are obtained in our KH synthetics since we do not consider multiple interactions of a ray with an interface. The nearly discontinuous behaviour of elastic R/T coefficients near the critical angle causes small spurious phases which arrive later than the correct arrivals.     

Propagation of SH-waves in a layered medium

Summary. A set of recurrence relations similar to that of Kennett suitable for SH -wave generation in an ( n + l)-layered half-space is presented. The recurrence relations contain no growing terms, thus providing a stable and efficient algorithm for computing complete SH synthetic seismograms. The complete expansion of these recurrence relations gives the explicit form of the transfer function for SH -waves. The transfer function for a point source in layer s of the stratification is a series of 2 ⁿ terms in the denominator and a series of 2 ^{n−s +1} terms in the numerator. The result of Wang from ray summation is shown to be a special case of our general result. Numerical comparison of the algorithm of this paper with the propagator matrix method is also made.     

Computation of high-frequency seismic wavefields in 3-D laterally inhomogeneous anisotropic media

Summary. An algorithm for the computation of travel times, ray amplitudes and ray synthetic seismograms in 3-D laterally inhomogeneous media composed of isotropic and anisotropic layers is described. All 21 independent elastic parameters may vary within the anisotropic layers. Rays and travel times are evaluated by numerical solution of the ray tracing equations. Ray amplitudes are determined by evaluating reflection/ transmission coefficients and the geometrical spreading along individual rays. The geometrical spreading is computed approximately by numerical measurement of the cross-sectional area of the ray tube formed by three neighbouring rays. A similar approximate procedure is used for the determination of the coefficients of the paraxial ray approximation. The ray paraxial approximation makes computation of synthetic seismograms on the surface of the model very efficient. Examples of ray synthetic seismograms computed with a program package based on the described algorithm are presented.     

A method for the computation of finite frequency body wave synthetic seismograms in laterally varying media

Summary. Body wave synthetic siesmograms for laterally varying media are computed by means of a slowness implementation of the extended WKBJ (EWKBJ) theory of Frazer & Phinney. An EWKBJ seismogram is computed by first tracing rays through a particular model to obtain conventional ray information (travel time, ray end point, ray slowness) and then using these data in the finite frequency integral expression for the EWKBJ seismogram. The EWKBJ seismograms compare favourably to geometrical ray theory (GRT) seismograms but are significantly better because of the finite frequency nature of the EWKBJ calculation. More realistic behaviour is obtained with EWKBJ seismograms at normal seismic frequencies near caustics, where the GRT amplitude is infinite, and within geometrical shadow zones where GRT predicts zero amplitudes. In addition the EWKBJ calculation is more sensitive than GRT to focuses and defocuses in the ray field. The major disadvantage of the EWKBJ calculation is the additional computer time over that of GRT, necessary to calculate one seismogram although an EWKBJ seismogram costs much less to compute than a reflectivity seismogram. Another disadvantage of EWKBJ theory is the generation of spurious, non-geometrical phases that are associated with rapidly varying lateral inhomogeneities. Fortunately the amplitudes of these spurious phases are usually much lower than that of neighbouring geometrical phases so that the spurious phases can usually be ignored. When this observation is combined with the moderately increased computational time of the EWKBJ calculation then the gain in finite frequency character significantly outweighs any disadvantages.     

A modified first-motion approximation for the synthesis of body-wave seismograms

Summary. Modified first-motion approximations have been developed for the generation of synthetic body-wave seismograms using the Cagniard-de Hoop method. Comparisons are presented between classical first motion, modified first motion and full Cagniard treatments for problems involving a homogeneous sphere and a triplication in a realistic earth model. Results of these comparisons show that the modified first-motion approximations may be used for a wide variety of geophysically interesting problems with little loss of accuracy compared to the full Cagniard method.     

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.     

Hierarchical traveltime tomography

A new traveltime tomographic method was developed with hierarchical shape functions of the finite element method as slowness or velocity interpolation functions. The degree of the approximation of velocity modelling is adjusted by selecting a set of hierarchical shape functions in each element. The ray density parameter of each element controls the selection to make the approximation fine or coarse in the high- or low-ray-density area. The proposed method is applied to both synthetic traveltime data and actual data. The AIC is used to determine the number of model parameters. The result of the synthetic data shows that low-resolution model parameters can be eliminated by the ray density parameter. The result of the actual data shows that the velocity pattern is approximately the same in the fine approximation area and that the velocity fluctuation is suppressed in the coarse approximation area, compared with that obtained from a full set of hierarchical shape functions. The number of model parameters is drastically reduced. The resolution can be estimated by the checkerboard restoration test. The result of the real data set was compared with that of the linear velocity grid model.     

Born scattering of longperiod body waves

The Born approximation is applied to the modelling of the propagation of deeply turning longperiod body waves through heterogeneities in the lowermost mantle. We use an exact Green's function for a spherically symmetric earth model that also satisfies the appropriate boundary conditions at internal boundaries and the surface of the earth. The scattered displacement field is obtained by a numerical quadrature of the product of the Green's function, the exciting wavefield and structural perturbations. We study three examples: scattering of longperiod P waves from a plume rising from the coremantle boundary (CMB), generation of longperiod precursors to PKIKP by strong, localized scatterers at the CMB, and propagation of corediffracted P waves through largescale heterogeneities in D". The main results are as follows: (1) the signals scattered from a realistic plume are small with relative amplitudes of less than 2 per cent at a period of 20 s, rendering plume detection a fairly difficult task; (2) strong heterogeneities at the CMB of appropriate size may produce observable longperiod precursors to PKIKP in spite of the presence of a diffraction from the PKP B caustic; (3) corediffracted P  waves ( P _diff) are sensitive to structure in D" far off the geometrical ray path and also far beyond the entry and exit points of the ray into and out of D"; sensitivity kernels exhibit ringshaped patterns of alternating sign reminiscent of Fresnel zones; (4) P _diff also shows a nonnegligible sensitivity to shear wave velocity in D"; (5) down to periods of 40 s, the Born approximation is sufficiently accurate to allow waveform modelling of P _diff through largescale heterogeneities in D" of up to 5 per cent.     

A linear approximation to the solution of a one-dimensional Stefan problem and its geophysical implications

Summary. The motion of a phase boundary in the Earth caused by temperature and pressure excitations at the Earth's surface is determined under a linear approximation. The solution is found as a sum of convolutions of pressure and temperature Green's functions with the corresponding excitations. The Green's functions are given under the form of Laplace transforms that can be inverted either by numerical evaluation of a branch cut integral or by inversion of a series expansion. This solution is a generalization of a solution previously derived by Gjevik. This latter solution is the first term in the series expansion. The relaxation times associated with the phase boundary motion are of the order of 10⁵–10⁷yr for the olivine—spinel phase transition and of 10⁶–10⁷yr for the basalt—eclogite transition. The linear approximation remains valid for long times only if the phase boundary moves slowly.     

A procedure for estimating lateral variations from low-frequency eigenspectra data

Summary. Expressions are derived for the zeroth and first moments of an isolated normal mode multiplet in the presence of rotational anisotropy, lateral heterogeneity and attenuation. It is demonstrated that, to first order, the zeroth moment is independent of the heterogeneity and the first moment is a linear functional of the heterogeneity. The location of the multiplet, defined in terms of the ratio of the first moment to the zeroth moment, is used to set up a linear inverse problem for lateral variations. The approximations employed in deriving the differential kernels are uniformly valid for all angular orders. In addition, formula for the differential kernels asymptotically valid in the limit of large angular order are developed, and the correspondence between these formulae and the geometrical optics approximation is established. It is shown that, in the limit of weak heterogeneity and large angular order, the location of the multiplet is the average of the local perturbation to the eigenfrequency over the great circular path containing the source and receiver.     

Total locational surplus for facility users distributed continuously along a network

This paper proposes an appropriate index for evaluating the degree to which facilities are desirable to their users under realistic assumptions. We derive an analytical expression for the index and discuss its calculation. Throughout this paper, we make the following three assumptions: (1) the distance from a user to a facility is measured by the shortest path, (2) facility users are continuously distributed along a network, and (3) the choice behavior of the facility users follows the Huff model. Under the third assumption, the index, called the ‘total locational surplus,’ enables us to conduct evaluations based on the concept of consumer surplus in microeconomics. We argue that this index is considered the most reliable among possible indices. First, we express the total locational surplus under our three assumptions as a definite integral. We then express a part of the integrand in the definite integral as an infinite series of exponential functions, which allows us to derive the infinite series expression of the definite integral. The infinite series expression enables us to calculate an approximate value for the total locational surplus using the first several terms of the series. We describe the computational procedure for calculating this approximation and evaluate its time complexity. We also utilize the procedure to evaluate the spatial configuration of some shopping centers in the northern part of Japan.     

Generalized Born scattering of elastic waves in 3-D media

It is well known that when a seismic wave propagates through an elastic medium with gradients in the parameters which describe it (e.g. slowness and density), energy is scattered from the incident wave generating low-frequency partial reflections. Many approximate solutions to the wave equation, e.g. geometrical ray theory (GRT), Maslov theory and Gaussian beams, do not model these signals. The problem of describing partial reflections in 1-D media has been extensively studied in the seismic literature and considerable progress has been made using iterative techniques based on WKBJ, Airy or Langer type ansätze. In this paper we derive a first-order scattering formalism to describe partial reflections in 3-D media. The correction term describing the scattered energy is developed as a volume integral over terms dependent upon the first spatial derivatives (gradients) of the parameters describing the medium and the solution. The relationship we derive could, in principle, be used as the basis for an iterative scheme but the computational expense, particularly for elastic media, will usually prohibit this approach. The result we obtain is closely related to the usual Born approximation, but differs in that the scattering term is not derived from a perturbation to a background model, but rather from the error in an approximate Green's function. We examine analytically the relationship between the results produced by the new formalism and the usual Born approximation for a medium which has no long-wavelength heterogeneities. We show that in such a case the two methods agree approximately as expected, but that in a media with heterogeneities of all wavelengths the new gradient scattering formalism is superior. We establish analytically the connection between the formalism developed here and the iterative approach based on the WKBJ solution which has been used previously in 1-D media. Numerical examples are shown to illustrate the examples discussed.     

Properties of S waves near a kiss singularity: a comparison of exact and ray solutions

We have studied the properties of S waves generated by a point source in a homogeneous, transversely isotropic, elastic medium, propagating in directions close to that of a kiss singularity, which coincides with the symmetry axis of the medium. We have proved analytically as well as numerically that the ray solution can describe the S waves correctly far from the source in all directions, including that of the kiss singularity. We have found that, in contrast to the far-field P wave, which can be reproduced satisfactorily by the zeroth-order ray approximation in all directions from the source, the far-field S waves can be reproduced satisfactorily by the zeroth-order ray approximation only in directions far from the kiss singularity. In directions near the kiss singularity, the first-order ray approximation must also be considered, because the zeroth- order ray approximation yields distorted results. The first-order ray approximation can be of high frequency and can be detected in the far field.     

A method of comparison of exact and asymptotic wave field computations

Summary. A method of comparison of exact numerical computations with an asymptotic ray series expansion consisting of the two first terms is proposed. The method makes it unnecessary to derive complicated explicit expressions for the second leading term of the ray series.
As a practical example we consider the anomalous PS arrival generated in the case of a near-vertical incidence of a spherical P wave on a solid/solid boundary. The areas in which the PS wave may be described by two leading terms of the ray series expansion are marked and deviations from the ray theory are analysed.     

Ray-reflectivity method for SH-waves in stacks of thin and thick layers

Summary Reflectivity and ray theories are united to produce a hybrid technique of computing synthetic seismograms for a plane layered medium in subcritical regions. Numerical experiments have indicated that this technique is useful when the depth structure is one composed of thick layers separated by finely layered zones. As the theory for wave propagation in a plane layered medium is well known, the simple SH case is investigated so that the basic idea of the method may be conveyed without an excess of mathematics that would be necessitated if the P-SV problem were considered.
In computing the ray-reflectivity seismogram, the thick layers are treated using asymptotic ray theory while the thin-layered zones are treated as quasiinterfaces where analogues of reflection and transmission coefficients called reflectivities and transmittivities are calculated utilizing a Thomson-Haskell formulation. A stationary phase approximation is employed when evaluating the integral which gives the displacement due to an arbitrary ray propagating in the thick layers of the above-mentioned medium, and the validity of this approximation is discussed.
A comparison of ray, numerical integration (reflectivity) and ray-reflectivity synthetic sections indicates that this method yields quite acceptable results for subcritical reflection work and is suitable for application in seismic interpretation as individual arrivals associated with ray-paths in the thick layers may be identified. Furthermore, the method is quite cost efficient and may be extended to a medium where the thick layers are non-planar using asymptotic ray theory in these layers.     

Validity of the great circular average approximation for inversion of normal mode measurements

Summary. Synthetic seismograms based upon first-order perturbation theory are analysed to test the validity of assumptions which form the basis of current velocity inversion procedures. It is found that the lowest order geometrical optics approximation, namely that measured normal mode eigen-frequencies reflect the average structure underlying the source–receiver great circle path, becomes less valid near nodes in the source radiation pattern and near the surface wave foci at the source and its antipode. These failures are a consequence of singlet interference within an isolated normal mode multiplet. The technique of determing frequency by fitting a single resonance peak to a multiplet yields results which agree well with the first-order theory for slow and fast paths where excitation is dominated by one pair of singlets but on intermediate paths where singlet interference is more of a problem, agreement is not as good. Inversion of small data sets is particularly sensitive to frequency fluctuations near radiation nodes, while larger sets are influenced more by antipodal deviations from geometrical optics. The latter leads to inversions which fail to recover the short wavelength structure of the starting model. Basing inversions directly upon first-order theory shows promise of improving recovery of short wavelengths.     

Finite-frequency vectorial tomography: a new method for high-resolution imaging of upper mantle anisotropy

Amplitude measurements of the transverse component of SKS waves, the so-called splitting intensity, can be used to formulate a non-linear inverse problem to image the 3-D variations of upper mantle anisotropy. Assuming transverse isotropy (or hexagonal symmetry), one can parametrize anisotropy by two anisotropic parameters and two angles describing the orientation of the symmetry axis. These can also be written as two collinear pseudo-vectors. The tomographic process consists of retrieving the spatial distribution of these pseudo-vectors, and thus resembles surface wave vectorial tomography. Spatial resolution results from the sensitivity of low-frequency SKS waves to seismic anisotropy off the ray path. The expressions for the 3-D sensitivity kernels for splitting intensity are derived, including the near-field contributions, and validated by comparison with a full wave equation solution based upon the finite element method. These sensitivity kernels are valid for any orientation of the symmetry axis, and thus generalize previous results that were only valid for a horizontal symmetry axis. It is shown that both lateral and vertical subwavelength variations of anisotropy can be retrieved with a dense array of broad-band stations, even in the case of vertically propagating SKS waves.     

The theory of finite frequency body wave synthetic seismograms in inhomogeneous elastic media

b
A new method is presented by means of which one can compute finite frequency synthetic seismograms for media whose velocity and density are continuous functions of two or three spatial variables. Basically, the method is a generalization of the familiar phase integral method, to which it reduces in a stratified medium. For a given source location the travel-time and distance functions needed to compute synthetics are obtained by numerically tracing rays through the model. This information is then used to evaluate a double integral over frequency and take-off angle at the source. The solution obtained reduces to the geometrical optics solution wherever that is valid but it also works in shadows and at caustics without knowing explicitly where these may be located. The method can be used as a spectral method, in which the integral over take-off angle is evaluated first, or as a slowness method, in which the frequency integral is evaluated first.     

The construction of effective methods for electromagnetic modelling

Summary. This paper deals with the further development of finite-difference methods for electromagnetic field modelling in two-and three-dimensional cases. The main feature of the approach suggested here is the application of generalized asymptotic boundary conditions valid with the accuracy (1/ρ^N), where ρ is the distance from the heterogeneities. The finite-difference approximation of problems under solution is made using the balance method, which results in 5-point difference schemes in the 2-D case and 7-point difference schemes in the 3-D case. To solve the linear system of difference equations the successive over-relaxation (SOR) method is used, the relaxation factor being chosen during the iteration procedure. In view of the vectorial character of the problem for the 3-D case, a successive blocked over-relaxation method (SBOR) is applied.
The model's validity is based on the comparison of the fields accounted at the ground surface with those computed by the integral transformation of excessive currents, determined in the heterogeneity region using the finite-difference scheme.