Elastic scattered waves from a continuous and heterogeneous layer

Elastic scattering from a continuous and laterally unbounded heterogeneous layer has been formulated using the Born approximation. A general solution of the scattered wave equation for the above-stated medium has been given in terms of a Fourier integral over plane waves. Far-field asymptotic expressions for weak elastic scattering by a finite, continuous and inhomogeneous layer have been presented which agree with earlier results. For perturbations of the two elastic parameters and the density having the same form of spatial variation, the spectrum of plane waves scattered from a heterogeneous layer is expressed as a product of an 'elastic scattering factor'and a 'distribution factor'. As in earlier results for small-scale heterogeneity, the scattering pattern depends on various combinations of perturbations of elastic parameters and density. In order to show the general characteristics of the elastic wave scattering, some scattering patterns have been given.     

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.     

Exact synthetic seismograms for an inhomogeneity in a layered elastic half-space

Summary. The propagation of a pulsed elastic wave in the following geometry is considered. An elastic half-space has a surface layer of a different material and the layer furthermore contains a bounded 3-D inhomogeneity. The exciting source is an explosion, modelled as an isotropic pressure point source with Gaussian behaviour in time.
The time-harmonic problem is solved using the null field approach (the T matrix method), and a frequency integral then gives the time-domain response. The main tools of the null field approach are integral representations containing the free space Green's dyadic, expansions in plane and spherical vector wave functions, and transformations between plane and spherical vector wave functions. It should be noted that the null field approach gives the solution to the full elastodynamic equations with, in principle, an arbitrarily high accuracy. Thus no ray approximations or the like are used. The main numerical limitation is that only low and intermediate frequencies, in the sense that the diameter of the inhomogeneity can only be a few wavelengths, can be considered.
The numerical examples show synthetic seismograms consisting of data from 15 observation points at increasing distances from the source. The normal component of the velocity field is computed and the anomalous field due to the inhomogeneity is sometimes shown separately. The shape of the inhomogeneity, the location and depth of the source, and the material parameters are all varied to illustrate the relative importance of the various parameters. Several specific wave types can be identified in the seismograms: Rayleigh waves, direct and reflected P -waves, and head waves.     

Time-domain computation of scattering from a heterogeneous layer

Adopting Born and ray approximations, time-domain synthetic seismograms for P-P and P-S scattering from a plane wave incident on a thin, laterally heterogeneous layer are presented in this paper. The time-domain P coda is a convolution between a structure function and the second-order derivative of the time function of the incident P wave. Examples of synthetic seismograms are given using a time function from a computed short-period seismogram for a point explosive source in a half-space. These show that it is impossible, with realistic values of the parameters involved, to generate significant codas when only single scattering is involved.     

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.     

The Love wave scattering matrix for a continental margin (numerical)

Summary . This paper presents the numerical computation of the results previously obtained by the author through a scattering matrix formulation (together with plane wave and variational approximations) which describes the diffraction of plane, harmonic, monochromatic Love waves incident normally (from either side) upon the plane of discontinuity in a structure consisting of a half-space with a surface step — an idealized model of a continental margin. Magnitudes of reflection and transmission coefficients are computed numerically for different frequencies for a model which has been considered previously by Knopoff & Hudson and also by Alsop in their studies of the same problem. The results obtained under the plane wave approximation are compared with those obtained under the variational approximation in order to assess the effects of the body-wave contributions. Finally, the results of both approximations are compared with those obtained by previous authors.     

Time-lapse traveltime change of singly scattered acoustic waves

We present a technique based on the single-scattering approximation that relates time-lapse localized changes in the propagation velocity to changes in the traveltime of singly scattered waves. We describe wave propagation in a random medium with homogeneous statistical properties as a single-scattering process where the fluctuations of the velocity with respect to the background velocity are assumed to be weak. This corresponds to one of two end-member regimes of wave propagation in a random medium, the first being single scattering, and the second multiple scattering. We present a formulation that relates the change in the traveltime of the scattered waves to a localized change in the propagation velocity by means of the Born approximation for the scattered wavefield. We validate the methodology with synthetic seismograms calculated with finite differences for 2-D acoustic waves. Potential applications of this technique include non-destructive evaluation of heterogeneous materials and time-lapse monitoring of heterogeneous reservoirs.     

Full reconstruction of a layered elastic medium from P-SV slant-stack data

Summary. We develop a méthod of reconstructing the elastic paraméters as functions of depth, for a horizontally stratified, isotropic elastic half-space. Unlike previous schemes, which have been able to retrieve the shear wave speed and density from SH seismograms slant stacked at two angles, our méthod makes use of P - SV data at a single stacking paraméter to obtain all three elastic constants. The data required are the elements of the full reflection matrix at the surface, corresponding to measurements of two separate components of the response to two independent sources, one explosive, the other generating shear waves.
In developing this inverse scheme fundamental differences emerge between the acoustic or SH problem, and the coupled P - SV case, the most important being in the nature of the interfacial scattering matrix. We show that it is not possible to make use of the downward reflection data for an interface to determine directly the remaining reflection and transmission coefficients, but that the scattering data may be completed by applying a simple iterative procedure at each interface.
We show the result of applying our inverse scheme to seismograms generated for a six-layered model, including a low-velocity layer. We are able to reconstruct both wave speeds and the density as functions of depth, all quantities being in close agreement with the original model.     

Shear-wave polarizations on a curved wavefront at an isotropic free surface

Summary. We present polarization diagrams of the particle motions at the free surface of an isotropic half-space generated by incident shear waves from a local buried point source. The reflectivity technique is used to calculate synthetic seismograms from which the particle motions are plotted. The particle motions are examined over a range of epicentral distances in a uniform isotropic half-space for different source frequencies and polarization angles, and for different Poisson's ratios. The particle motions due to a curved wavefront possess different characteristics from those generated by plane wavefronts at corresponding incidence angles. A curved wavefront generates a local SP -phase: a P -headwave which propagates along the free surface, and arrives shortly before the direct S -wave. These two arrivals give rise to cruciform particle motions in the sagittal and horizontal planes, which could be misinterpreted as anisotropy-induced shear-wave splitting. An examination of the particle motion in the transverse plane, mutually orthogonal to the sagittal and horizontal planes, can be used to discriminate between isotropic and anisotropic interpretations. The amplitude of the SP -phase is enhanced when it propagates in a low-velocity surface layer overlying the source layer, and may then become the dominant phase on radial-component seismograms. The presence of even a single surface layer may introduce considerable complexity into the seismogram, and we examine the effects of layer thickness, velocity contrast, and source depth on the corresponding polarization diagrams. Reliable information on the source and propagation path characteristics of shear waves from a buried local point source can only be obtained from free-surface records if they are recorded within a very limited epicentral distance range.     

Fault plane solutions using relative amplitudes of P and pP

Summary. One way of finding the fault plane orientations of small shallow earthquakes is by the generation of theoretical P -wave seismograms to match those observed at several distant stations. Here, a technique for determining the uniqueness of fault plane solutions computed using the modelling method of Douglas et al . is described. Relative amplitudes of P and pP , and their polarities if unambiguous, are measured on the observed seismograms to be modelled, and appropriate confidence limits are assigned to each measurement. A systematic search is then made for all fault plane orientations which satisfy these observations.
Examples show that if P and pP are not severely contaminated by other arrivals, a well-defined and unique fault plane orientation can often be computed using as few as three stations well distributed in azimuth. Further, even if pP is not identifiable on a particular seismogram, then an upper bound on its amplitude – deduced from the observed coda – still places a significantly greater constraint on the fault plane orientation than would be provided by a P onset polarity alone. Modelling takes account of all such information, and is able to further eliminate incompatible solutions (e.g. by the correct simulation of sP ). It follows that if solutions can be found which satisfy many observed seismograms, this places high significance on the validity of the assumed double-couple source mechanism.
This relative amplitude technique is contrasted with the familiar first motion method of fault plane determination which requires many polarity readings, whose reliabilities are difficult to quantify. It is also shown that fault plane orientations can be determined for earthquakes below the magnitude at which first motion solutions become unreliable or impossible.     

Investigation of 3D wavefields of reflected shear-waves and converted waves: mathematical modelling and reflection data processing

Summary. Seismic investigations using shear-wave and converted wave techniques show that very often reflected PS - and SS -waves have anomalous polarizations ( accessory components ). This phenomenon cannot be explained in terms of isotropic models with dipping boundaries. Computations of synthetic seismograms of reflected PS - and SS -waves were made for different models of transversely isotropic media with dipping anisotropic symmetry axes not normal to the boundaries. Synthetic seismograms were computed by ray techniques using an optimization algorithm to construct all rays arriving at a given receiver. These computations indicate that accessory components arise when the medium above the boundary is anisotropic, where they are caused by the constructive interference of qSV - and qSH -waves. If a low-velocity layer is present, displacement vectors of both waves have horizontal projections which are approximately orthogonal. The algorithm for wave separation is presented and some results of its use are given.     

Reflection coefficients for weak anisotropic media

The interaction of plane elastic waves with a plane boundary between two anisotropic elastic half-spaces is investigated. The anisotropy dealt with in this study is of a general type. Explicit expressions for energy-related reflection and transmission coefficients are derived. They represent an approximation which is valid for a small deviation of the elastic parameters from isotropy.
Classical perturbation theory is applied on a 6times6 non-symmetric real eigenvalue problem to calculate first-order corrections for the polarization and stress of the plane waves. The explicit solution of the isotropic problem is used as a reference case. Degenerate perturbation theory is used to consider the splitting of the isotropic S -wave into two anisotropic qS-waves. The boundary conditions for two half-spaces in welded contact lead to a 6times6 system of linear equations. A correction to the isotropic solution is calculated by linearization. The resultant coefficients are functions of horizontal slowness, Lamé parameters and densities of the reference media, and of the perturbation of the elasticity tensors from isotropy.     

Observation and modelling of fault-zone fracture seismic anisotropy - I. P, SV and SH travel times

Summary. Three-component VSP borehole seismograms taken in the vicinity of an active normal fault in California show strong systematic shear-wave splitting that increases with proximity to the fault. Using Červený's method of characteristics for ray tracing in anisotropic heterogeneous media and Hudson's formulation of elastic constants for media-bearing aligned fractures, we have fitted a suite of P, SV and SH hanging-wall and foot-wall travel times with a simple model of aligned fractures flanking the fault zone. The dominant fracture set is best modelled as parallel to the fault plane and increasing in density with approach to the fault. The increase in fracture density is non-uniform (power law or Gaussian) with respect to distance to the fault. Although the hanging-wall and the foot-wall rock are petrologically the same unit, the fracture halo is more intense and extensive in the hanging wall than in the foot wall. Upon approach to the fault plane, the fracture density or fracture-density gradient becomes too great for the seismic response to be computed by Hudson–Červený procedures (the maximum fracture density that can be modelled is about 0.08). Within this 25 m fracture domain it appears more useful to model the fault and near field fractures as a low-velocity waveguide. We observe production of trapped waves within the confines of the intense fracture interval.     

A fault plane solution using theoretical P seismograms

We use the method of Hudson and Douglas, Hudson & Blarney to compute seismograms which simulate the codas of 10 short period P -wave seismograms from a shallow earthquake. The polarities and relative amplitudes of P and pP measured from seven of the observed seismograms are used to compute a fault plane solution with confidence limits, assuming that the source radiates as a double couple. This solution is in approximate agreement with that given for the same earthquake by Sykes & Sbar, who used only the onset polarities of short-period P waves. The small difference between the two solutions can be explained by interference between the true first motion of P and microseismic noise at two stations.
The results show that, for some shallow earthquakes, the relative amplitude method has the following advantages over the first motions method. First, a P/pP amplitude ratio (with appropriate confidence limits) can always be measured, even in seismograms which are so noisy that the first motion of P is uncertain. Second, the fault plane solutions obtained from relative amplitudes have known confidence limits. Finally, by using more information from each seismogram, the relative amplitude method requires considerably fewer seismograms than the first motions method.     

Numerical simulation of the propagation of P waves in fractured media

We study the propagation of P waves through media containing open fractures by performing numerical simulations. The important parameter in such problems is the ratio between crack length and incident wavelength. When the wavelength of the incident wavefield is close to or shorter than the crack length, the scattered waves are efficiently excited and the attenuation of the primary waves can be observed on synthetic seismograms. On the other hand, when the incident wavelength is greater than the crack length, we can simulate the anisotropic behaviour of fractured media resulting from the scattering of seismic waves by the cracks through the time delay of the arrival of the transmitted wave. The method of calculation used is a boundary element method in which the Green's functions are computed by the discrete wavenumber method. For simplicity, the 2-D elastodynamic diffraction problem is considered. The rock matrix is supposed to be elastic, isotropic and homogeneous, while the cracks are all empty and have the same length and strike direction. An iterative method of calculation of the diffracted wavefield is developed in the case where a large number of cracks are present in order to reduce the computation time. The attenuation factor Q ⁻¹ of the direct waves passing through a fractured zone is measured in several frequency bands. We observe that the attenuation factor Q ⁻¹ of the direct P wave peaks around kd = 2, where k is the incident wavenumber and d the crack length, and decreases proportionally to ( kd ) ⁻¹ in the high-wavenumber range. In the long-wavelength domain, the velocity of the direct P wave measured for two different crack realizations is very close to the value predicted by Hudson's theory on the overall elastic properties of fractured materials.     

End-point contributions to synthetic seismograms

Summary . Synthetic seismograms represented by integrals generally display signals associated with the limits of integration. Sometimes these 'end-point' contributions are spurious (e.g. in the WKBJ seismogram) and sometimes they are the main physical interest (e.g. the Kirchhoff integral for an edge). The end-point contributions may be asymptotically approximated using integration by parts or Laplace's method and it may then be possible to reduce them if desired. We describe examples in the WKBJ seism ogram for reflected or transmitted waves in homogeneous layers and for turning waves. We also study signals due to discontinuities in reflection coefficients, by partitioning the real slowness integral so that the discontinuities lie at end points. Examples are the head wave, which is a physically correct signal, and spurious diffractions caused by using plane-wave coefficients for grazing rays in the WKBJ seismogram.     

Generation of complete theoretical seismograms for SH— II.

Summary. We report the results of our continuing efforts to compute theoretical seismograms for direct comparison with the experimental time series obtained with the long-period instruments of the WWSSN. The entire theoretical seismogram — body waves and surface waves — is generated for realistic sources buried in a radially heterogeneous, anelastic, spherical earth. The results described in Paper I (Nakanishi, Schwab & Knopoff) are extended to include the summation of 11 modes; for each, the dispersion, attenuation, and excitation are computed down to a minimum period of 1 s. Examples of the theoretical seismograms, and the comparison with experimental results are presented, The results of this comparison indicate that our first application of combined body- and surface-wave generation will concern the investigation of the intrinsic anelasticity in the upper mantle. The indicated technique for such an investigation is based on body waves simply crossing the region of high attenuation a few times in passing from focus to recording station, while a guided surface wave such as Sa , experiences this anelasticity over the entire propagation path.     

3-D linearized scattering of surface waves and a formalism for surface wave holography

Summary. Scattering of surface waves by lateral heterogeneities is analysed in the Born approximation. It is assumed that the background medium is either laterally homogeneous, or smoothly varying in the horizontal direction. A dyadic representation of the Green's function simplifies the theory tremendously. Several examples of the theory are presented. The scattering and mode conversion coefficients are shown for scattering of surface waves by the root of an Alpine-like crustal structure. Furthermore a 'great circle theorem'in a plane geometry is derived. A new proof of Snell's law is given for surface wave scattering by a quarter-space. It is shown how a stationary phase approximation can be used to simplify the Fourier synthesis of the scattered wave in the time domain. Finally a procedure is suggested to do 'surface wave holography'.     

The accuracy of models derived by WKBJ waveform inversion

Summary . In this paper the accuracy of velocity-depth profiles derived by matching WKBJ seismograms to observations is quantitatively evaluated. Seismograms computed with the WKBJ method are generally quite reliable but possess predictable, systematic inaccuracies in the presence of strong velocity gradients. The effects of these inaccuracies on models derived through WKBJ waveform inversion are studied, using reflectivity seismograms as 'data'. The velocity structure used is an oceanic lithosphere model that contains several transition regions separated by relatively homogeneous layers, producing partially-reflected reverberations in the reflectivity synthetics that are absent from the WKBJ seismograms. The inversion incorporates the 'jumping' strategy to solve for the smoothest models consistent with the data. We find these solutions to be independent of the starting model and to have a stable basic structure that agrees well with the correct model. The differences, everywhere less than a seismic wavelength, depend on the frequency content of the seismograms. Reverberations in the reflectivity seismograms that are well separated from WKBJ arrivals are treated as 'noise' in the inversion.     

Synthetic body wave seismograms for three-dimensional laterally varying media

Summary. Two methods of computing body wave synthetic seismograms in three-dimensional laterally varying media are discussed. Both these methods are based on the summation of Gaussian beams. In the first, the initial beam parameters are chosen at the source, in the second at the beam endpoints. Both these variants eliminate the ray method singularities. The expansion of the wavefield into plane waves may be considered as the limiting case of the first approach and the Chapman–Maslov method as the limiting case of the second approach. Computer algorithms are briefly described and numerical examples presented. In the first numerical example, the comparisons of the two approaches, based on summing Gaussian beams, with the reflectivity method indicate that the computed synthetic seismograms are satisfactorily accurate even in the caustic region. The next example suggests that the two methods discussed can be simply and effectively applied to 3-D laterally inhomogeneous structures.