The S receiver functions: synthetics and data example

Recently, the S receiver function method has been successfully developed to identify upper mantle interfaces. S receiver functions have the advantage of being free of S -wave multiple reflections and can be more suitable than P receiver functions for studying mantle lithosphere. However, because of specific ray geometry and interference of diverse phases, the S receiver function method has some technical difficulties and limitations. We use synthetic seismograms to demonstrate the feasibility and limitations of S receiver functions for studying mantle structures. Full-wavefield seismograms were calculated using the reflectivity method and processed to generate synthetic S receiver functions for S , SKS and ScS waves. Results show that S receiver functions can be obtained from waveforms of S , SKS and ScS waves. The synthetic S receiver functions for these incident waves show S -to- P converted phases at all discontinuities in the crust and upper mantle. Useful ranges of epicentral distances for calculation of S receiver functions are: 55°–85° for S , >85° for SKS and 50°–75° for ScS waves. We apply both the S and P receiver function methods to data recorded at broadband station YKW3 in Northwest Canada. The study shows that there is significant agreement among different receiver function methods, and demonstrates the usefulness of S receiver functions for imaging the mantle lithosphere.     

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.     

A lower mantle S-wave triplication and the shear velocity structure of D"

Summary. A lower mantle S-wave triplication detected with short- and long-period WWSSN and CSN recordings indicates a substantial shear velocity discontinuity near 280 km above the core–mantle boundary. The triplication can be observed in rotated SH seismograms from intermediate and deep focus events throughout the distance range from 70° to 95°. Three distinct source region–receiver array combinations that have been investigated in detail demonstrate consistent travel time and relative amplitude behaviour of the triplication, with slight systematic shifts in the triplication indicating up to 40 km variations in the depth of the discontinuity. Modelling of the observations with synthetic seismograms produced with the Cagniard de Hoop and reflectivity methods constrains the shear velocity increase to be 235 ± 0.25 per cent, comparable to upper mantle discontinuities. Short-period observations indicate that the velocity increase may be a sharp first-order discontinuity, or may extend over a transition zone no more than 50 km thick. The shear velocity gradient below the discontinuity, within the D" layer, is not well-constrained by the SH data, but slightly positive or near zero velocity gradients are consistent with the long-period amplitude ratios of ScSH/SH .     

2-Dreflectivity method and synthetic seismograms for irregularly layered structures -I. SH-wave generation

Summary. The reflectivity method for complete SH seismograms has been extended to two-dimensionally layered structures. The Aki-Larner technique is generalized to solve the integral equations for 2-D boundary conditions, and propagator matrices are enlarged to express a total SH wavefield. Synthetic seismograms in a soft basin are calculated for an incident plane-wave. They compare favourably with the results of the finite-element and finite-difference methods even in the later portion where asymptotic ray and beam theories break down. Synthetic seismograms due to a line force and a point dislocation are also presented.     

On a generalization of Filon's method and the computation of the oscillatory integrals of seismology

Summary. We review Filon's method (FM) for the quadrature of oscillatory integrals and then introduce a generalization of Filon's method (the GFM) which enables us to treat a large class of oscillatory integrals to which FM cannot be directly applied. One member of this class is the integral ( p ) exp [ sg ( p )] dp which occurs in the spectral WKBJ and Cagniard-de Hoop methods of seismogram synthesis. Another large class of integrals can be treated directly with FM but is better treated with the GFM since, for a given error tolerance, the GFM is simpler and faster. This class consists of integrals of the form ( p ) J ( s, p ) dp in which J ( s, p ) is a special function with an asymptotic expansion valid for large s. Such integrals occur in the reflectivity method. In general, every non-Filon formula for the quadrature of integrals from either class has an associated GFM formula (called the GFM analogue) which reduces to the original formula as s approaches zero but is more efficient than the original formula wher, s is large. We show how the GFM can be applied to the computation of synthetic seismograms in the reflectivity method and the spectral WKBJ method.
Although reflectivity integrals can, in theory, be computed with FM the GFM is easier to code and more economical. For reflectivity computations where: (a) the source and receiver are many wavelengths apart, or (b) the depth to the reflectivity zone is much greater than its thickness, the GFM approach is much more efficient than any non-Filon quadrature technique. Some test calculations are presented for wavefields containing only body waves and for wavefields containing both body waves and locked modes.
In the spectral WKBJ method the GFM permits the use of a much greater step size in the quadrature than would otherwise be possible. Each quadrature step contains a stationary point so no advantages accrue from deforming the contour of integration over the saddle points of the integrand.     

A slowness approach to the reflectivity method of seismogram synthesis

b
Of the many schemes available for computing synthetic seismograms, the reflectivity method is probably the most widely used because of its ability to provide complete solutions. The method does, however, suffer the disadvantage that intermediate results are quite difficult to interpret. A new reflectivity technique, here called reflectivity-slowness, results if the original method is reformulated using a slowness rather than a spectral approach. The new procedure bears a strong similarity to the WKBJ method, but retains the ability to give complete solutions. The reflectivity-slowness and WKBJ methods share the property that intermediate results are readily interpreted; this feature may eventually be exploited in the solution of the inverse problem.     

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.     

Reflectivity of the crystalline crust: hypotheses and tests

Summary. The nature of reflectors within the crystalline basement remains the subject of inference except where reflections have been traced directly to outcrop. Geological models of basement reflectors need to be developed which incorporate geophysical constraints obtained from measurements on seismograms, but most geological information still comes from speculative interpretations of seismic experiments run in different regimes. Pronounced lower-crustal reflectivity, detected worldwide, is ascribed in various geological hypotheses to primary lithologic layering, to ductile strain banding, or to trapped fluids. A BIRPS deep crustal profile across the Atlantic continental margin suggests that the observed reflectivity is not related in any simple way to the amount of extensional strain undergone. Study of worldwide crustal profiles shows that exposed high-pressure terranes are not as reflective as in situ lower crust at high pressure, suggesting either that these granulite terranes are not representative of the lower crust or that physical properties, possibly the presence of fluids or thermally controlled ductile strain banding, are more likely responsible for observed reflectivity than are simple lithologic boundaries. The argument for the importance of physical properties in causing observed lower-crustal reflectivity is strengthened by an observed negative correlation between depth to the reflective lower crust and regional surface heat-flow.     

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.     

Least-squares fitting of marine seismic refraction data

Summary. An iterative procedure is presented for fitting waveform data from a marine seismic refraction experiment. The wavefunction from the explosive source is known and the crustal structure is refined using the damped least squares procedure. The damping parameter serves the dual purpose of stabilizing an under-constrained inversion and improving the linearity by suppressing high frequencies. The synthetic seismograms and their model differentials are calculated using the WKBJ seismogram method. Both the synthetic seismograms and the linear algebra are sufficiently straightforward that the computations can be performed on an array processor. The inversion procedure is then sufficiently rapid that interactive computations are possible. The technique is illustrated using the FF2 refraction data from the Fanfare cruise of the Scripps Institution of Oceanography. These data had been interpreted previously by trial-and-error using the reflectivity method. Starting from two different, simple models, the inversion procedure obtains essentially one unique model. Its features are very similar to the previous model.     

Efficient calculation of Gaussian-beam seismograms for two-dimensional inhomogeneous media

Summary. This paper discusses several aspects of the calculation of theoretical seismograms for two-dimensional inhomogeneous media with the method of Gaussian beams. The most important steps of this method, kinematic and dynamic ray tracing, can be performed very efficiently, if the model cross-section is subdivided into triangles with linear velocity laws. Each Gaussian beam is characterized by a complex beam constant ε which determines its width and phase-front curvature. Various possibilities to choose ε are discussed, including cases where beam properties at the beam endpoint (and not at the beginning) are prescribed; for instance, the beam width at the endpoint can be specified. In such cases the beam constant is a function of the radiation angle at the source, and the decomposition of a cylindrical wave into beams has to take this into account by weighting the beams differently, at least in principle. The exact weight function is derived and shown to be reasonably well approximated by the weight function, corresponding to angle-independent ε Theoretical seismograms are presented for a laterally heterogeneous model of the crust–mantle transition which is characterized by complications in the reflection from the transition and in the refraction from below. These complications are modelled by and large with success. The seismograms, however, depend to a certain extent on the choice of the beam constant. Moreover, according to the reciprocity principle calculations with source and receiver interchanged should have the same results as calculations for the original configuration. In practice this is not so, and the difference increases with the strength of lateral heterogeneities. Hence, for a successful application of Gaussian beams the model should not vary too strongly in lateral direction.     

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.     

Recurrence relations for computing complete P and SV seismograms

Summary. A set of recurrence relations which are computationally more efficient than those of the reflection matrix method of Kennett & Kerry is presented for P - and SV -wave generation in a ( n + 1) layered medium. The recurrence relations contain no growing terms and thus provide a stable algorithm for computing complete P and SV synthetic seismograms. Our algorithm requires a fewer algebraic operations for computing the reflectivity and transmissivity coefficients, ranging from 15 per cent less for a source in the half-space to 30 per cent less for a source in the top layer, than the reflection matrix method.     

Waveform modelling of teleseismic S, Sp, SsPmP, and shear-coupled PL waves for crust- and upper-mantle velocity structure beneath Africa

We describe a waveform modelling technique and demonstrate its application to determine the crust- and upper-mantle velocity structure beneath Africa. Our technique uses a parallelized reflectivity method to compute synthetic seismograms and fits the observed waveforms by a global optimization technique based on a Very Fast Simulated Annealing (VFSA). We match the S , Sp, SsPmP and shear-coupled PL phases in seismograms of deep (200–800 km), moderate-to-large magnitude (5.5–7.0) earthquakes recorded teleseismically at permanent broad-band seismic stations in Africa. Using our technique we produce P - and S -wave velocity models of crust and upper mantle beneath Africa. Additionally, our use of the shear-coupled PL phase, wherever observed, improves the constraints for lower crust- and upper-mantle velocity structure beneath the corresponding seismic stations. Our technique retains the advantages of receiver function methods, uses a different part of the seismogram, is sensitive to both P - and S -wave velocities directly, and obtains helpful constraints in model parameters in the vicinity of the Moho. The resulting range of crustal thicknesses beneath Africa (21–46 km) indicates that the crust is thicker in south Africa, thinner in east Africa and intermediate in north and west Africa. Crustal P - (4.7–8 km s⁻¹) and S -wave velocities (2.5–4.7  km s⁻¹) obtained in this study show that in some parts of the models, these are slower in east Africa and faster in north, west and south Africa. Anomalous crustal low-velocity zones are also observed in the models for seismic stations in the cratonic regions of north, west and south Africa. Overall, the results of our study are consistent with earlier models and regional tectonics of Africa.     

Numerical modelling of propagation of elastic waves in anisotropic inhomogeneous media for the half-space and the sphere

Summary. A method based on a combination of partial separation of variables and finite-difference method is used for the calculation of complete theoretical seismograms for inhomogeneous anisotropic media. Examples of theoretical seismograms for several anisotropic models are presented.     

Shear-coupled PL

Summary. Observed teleseismic shear-coupled PL -waves ( SPL ) display a variety of waveforms depending on factors such as source depth, source type and velocity structure. Using a WKBJ spectral method for SV -wave propagation, synthetic seismograms of SPL are produced to examine the factors important in SV and SPL excitation. Results show that SPL is preferentially excited by shallow sources compared to deep sources. This is due to large source area reverberations which consequently leak as SV -waves into the mantle. Interaction at the receiver area then sets up the classic prograde elliptical motion by which SPL can be identified. This is in accordance with the teleseismic observations and indicates that most previous models for the propagation of SPL were not appropriate for shallow source since emphasis was placed on wave interactions occurring only near the receiver.     

The effect of noise on seismograms

Summary. We describe a method which provides an estimate of the accuracy to which time-domain features of seismic signals can be measured in the presence of noise. Observed seismograms are simulated by adding random noise with the same frequency spectrum and signal-to-noise ratio to matching synthetic seismograms. The effect of noise on synthetic and observed P -wave first motions is used as an illustration. It is shown that the apparent reliability of such observations, as determined by visual estimation, is often illusory.     

Generation of the teleseismic P-wave coda from Aleutian earthquakes

We investigate large-amplitude phases arriving in the P -wave coda of broad-band seismograms from teleseisms recorded by the Gräfenberg array, the German Regional Seismic Network and the Global Seismic Network. The data set consists of all events m _b≤ 5.6 from the Aleutian arc between 1977 and 1992. Earthquakes with large-amplitude coda waves correlate with the presence of oceanic crust in the source region. The amplitudes sometimes approach those of the P wave, much larger than predicted by theory. Modelling indicates that phases in the P -wave coda cannot be P -wave multiples beneath the source and receiver, or underside reflections, which precede PP , from upper-mantle discontinuities. Among the events, seismograms are very similar, where the arrival times of the unusual phases agree approximately with the predicted times of S -to- P conversions from the upper-mantle discontinuities under the source. Because the large-amplitude phases in the P -wave coda have little, if any, dependence on event depth and have predominantly an SV -wave radiation pattern towards the receiver, we suggest that they originate as SV and/or Rayleigh waves and are enhanced by lateral heterogeneity and multipathing from the subducting Aleutian slab.     

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.     

P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation

P-SH conversion is commonly observed in teleseismic P waves, and is often attributed to dipping interfaces beneath the receiver. Our modelling suggests an alternative explanation in terms of flat-layered anisotropy. We use reflectivity techniques to compute three-component synthetic seismograms in a 1-D anisotropic layered medium. For each layer of the medium, we prescribe values of seismic velocities and hexagonally symmetric anisotropy about a common symmetry axis of arbitrary orientation. A compressional wave in an anisotropic velocity structure suffers conversion to both SV -and SH -polarized shear waves, unless the axis of symmetry is everywhere vertical or the wave travels parallel to all symmetry axes. The P-SV conversion forms the basis of the widely used 'receiver function' technique. The P-SH conversion occurs at interfaces where one or both layers are anisotropic. A tilted axis of symmetry and a dipping interface in isotropic media produce similar amplitudes of both direct ( P ) and converted ( Ps ) phases, leaving the backazimuth variation of the P-Ps delay as the main discriminant. Seismic anisotropy with a tilted symmetry axis leads to complex synthetic seismograms in velocity models composed of just a few flat homogeneous layers. It is possible therefore to model observations of P coda with prominent transverse components with relatively simple 1-D velocity structures. Successful retrieval of salient model characteristics appears possible using multiple realizations of a genetic-algorithm (GA) inversion of P coda from several backazimuths. Using GA inversion, we determine that six P coda recorded at station ARU in central Russia are consistent with models that possess strong (> 10 per cent) anisotropy in the top 5 km and between 30 and 43 km depth. The symmetry axes are tilted, and appear aligned with the seismic anisotropy orientation in the mantle under ARU suggested by SKS splitting.