Constraining velocity and density contrasts across the crust–mantle boundary with receiver function amplitudes

This paper investigates the ability of P -wave receiver functions to constrain both the velocity and density contrasts across the Moho. Expressions for receiver function amplitudes corresponding to a layer over a half-space are extended to media with depth-dependent properties by explicitly accounting for free-surface reflection coefficients. Forward numerical computations show that receiver function amplitudes become frequency-dependent for depth-dependent structures and that, for a given frequency, wavelengths associated to Ps converted phases are larger than wavelengths associated to multiply reverberated phases. When consistent measurements of the amplitudes of Ps phase and multiples can be obtained, the extended expressions suffice to account for the observed amplitudes. Receiver function amplitudes are sensitive to both velocity and density contrasts across the Moho, and a two-step grid-search procedure is proposed to recover the contrasts from their measurement. The near-surface velocity is recovered in the first step from the amplitude of the direct P wave and then used in the second step to work out the relative density and S -wave velocity contrast from the amplitude of the secondary phases. Examples from central Spain (PAB) and the Indian Shield (HYB) are discussed and demonstrate that receiver function amplitudes can constrain contrasts across the Moho accurately enough to be utilized in geological interpretation.     

Polarization and amplitude attributes of reflected plane and spherical waves

The characteristics of a reflected spherical wave at a free surface are investigated by numerical methods; in particular, the polarization angles and amplitude coefficients of a reflected spherical wave are studied. The classical case of the reflection of a plane P wave from a free surface is revisited in order to establish our terminology, and the classical results are recast in a way which is more suited for the study undertaken. The polarization angle of a plane P wave, for a given angle of incidence, is shown to be 90° minus twice the angle of reflection of the reflected S wave. For a Poisson's ratio less than 1/3, there is a non-normal incident angle for which both amplification coefficients are 2 precisely; for this incident angle the direction of the particle motion at the free surface is also the direction of the incident wave. For a wave emanating from a spherical source, the polarization angle, for all angles of incidence, is always less than, or equal to, the polarization angle of a plane P wave. The vector amplification coefficient of a spherical wave, for all angles of incidence, is always greater than the vector amplification coefficient of a plane P wave. As expected, the results for a spherical wave approach the results for a plane P wave in the far field. Furthermore, there was a good agreement between the theoretical modelling and the numerical modelling using the dynamic finite element method (DFEM).     

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.     

An analysis of P-wave amplitudes recorded by seismological stations in the USSR

Summary. Short period P -wave amplitudes at Soviet seismic stations are analysed to determine an amplitude–distance curve in the range 0°–180° from the USSR. Station amplitude terms, corrected for crustal amplification effects, are determined from seismic sources located between 30° and 90°. A relationship between these station terms and regional P_n wave speed and regional heat flow is determined and is found to be in close agreement with similar observations in North America. Station amplitude terms were found to be reasonably uniform across the continental portion of the USSR but the Baikal Rift valley is characterized by very low P -wave amplitudes and is the most significant feature discemible from the P -wave amplitude terms.     

Three-dimensional geometrical ray theory and modelling of transmitted seismic energy of data from the Nevada Test Site

This paper presents a geometrically based algorithm for computing synthetic seismograms for energy transmitted through a 3-D velocity distribution. 3-D ray tracing is performed to compute the traveltimes and geometrical spreading (amplitude). The formulations of both kinematic and dynamic ray-tracing systems are presented. The two-point ray-tracing problem is solved by systematically updating the initial conditions and adjusting the ray direction until the ray intersects the specified endpoint. The amount of adjustment required depends on the derivatives of the position with respect to the given starting angles between consecutive rays. The algorithm uses derivatives to define the steepest-descent direction and to update the initial directions. The convergence rate depends on the complexity of the model.
Test seismograms compare favourably with those from a 2-D asymptotic ray theory algorithm and a 3-D Gaussian-beam algorithm. The algorithm is flexible in modelling arbitrary source and recorder geometries for various smoothly varying 3-D velocity distributions. The algorithm is further tested by simulating surface-to-tunnel vibroseis field data. Shear waves as well as compressional waves may be approximately included. Application of the algorithm to a data set from the Rainier Mesa of the Nevada Test Site produced a good fit to the transmitted (first arrival) traveltimes and amplitudes, with approximately 15 per cent variation in the local 3-D velocity.     

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.     

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.     

Modelling of scattering of seismic waves from a corrugated rough sea surface: a comparison of three methods

We compare three numerical methods to model the sea surface interaction in a marine seismic reflection experiment (the frequencies considered are in the band 10–100 Hz): the finite-difference method (FDM), the spectral element method (SEM) and the Kirchhoff method (KM). A plane wave is incident at angles of 0° and 30° with respect to the vertical on a rough Pierson–Moskowitz surface with 2 m significant wave height and the response is synthesized at 6, 10 and 50 m below the average height of the sea surface. All three methods display an excellent agreement for the main reflected arrival. The FDM and SEM also agree very well all through the scattered coda. The KM shows some discrepancies, particularly in terms of amplitudes.     

Phase-array decomposition of a seismic section

A seismic re fraction/wide-angle reflection profile is analysed for the presence of correlated events ('phases'). The correlation problem is formulated in terms of temporally, spatially and frequency-local complex covariances. For robustness, the method concentrates on phase rather than amplitude information. This allows a computationally efficient algorithm that can make allowance for signal correlation length and can model curved wavefronts. A statistical test based on residual phase misfit across the analysed subarray is used to assess the probability that a detected event represents a real correlated signal.
With our chosen analysis parameters and confidence level (over 99.9 per cent). 1222 events were detected in the data. Using simple techniques based on 1-D earth models, detected events are associated with a small number of particular wave types. In this way, we have succeeded in classifying almost 95 per cent of the detected events. Those that remain describe those components of the data that are inconsistent with our simple ray paths in the 1-D assumption and with our prescribed tolerance. These include reverberations, near-surface guided waves and reflected waves from strongly laterally inhomogeneous structures. According to our modelling, about 25 per cent of the detected events are consistent with simple P -wave reflected energy, and these are to a very large extent (over 85 per cent) distinct from all the other wave-type models we have used. A direct mapping of the detected events into the offset-depth domain reveals dear internal and external consistencies among the detections for the various wave types. Estimated earth structure is consistent with models from previous analyses based on much larger data sets.
We have thus succeeded in extracting correlated events from the data and decomposing these, approximately but meaningfully, into distinct classes (ray paths)     

Azimuthal variation of the P phase in Icelandic receiver functions

A curious observation has been made on radial receiver functions calculated from teleseisms recorded by 29 broad-band seismometers distributed over Iceland. The arrival time of the direct P phase of the radial receiver functions depends critically upon the azimuth of the teleseismic source. For a seismic station in West Iceland, the direct P  phase of the radial receiver function arrives consistently later for easterly source azimuths than for westerly source azimuths. The reverse applies for stations in East Iceland. In the original seismograms, the delayed P phase of the receiver function appears up to 450 ms later on the radial than on the vertical component. The seismometer locations in East and West Iceland are separated by the Neovolcanic Zone, a constructive plate boundary. The delayed P phases occur for seismic rays travelling across this zone. However, it is not obvious how wave propagation across the plate boundary zone could cause the observed delays. The tentative explanation proposed here involves the regional dip of the Icelandic lava sequences towards the Neovolcanic Zone. A dipping interface at shallow depth results in a P–S converted phase arriving shortly after the P phase. These phases cannot be separated in the radial receiver functions, given the bandwidth of the observed signals. However, a calculation of receiver functions from estimates of the P , SV and SH wavefields clearly reveals a P–S converted phase at about 500 ms for easterly source azimuths in West Iceland and for westerly source azimuths in East Iceland. The amplitudes of the direct P phase and the P–S phase converted at a dipping interface would be expected to vary strongly with azimuth in accordance with the observed behaviour.     

Anisotropic tomography of P-wave traveltimes

The ray path of a P -wave is specified in terms of the ray parameter and three Euler angles. the P -wave traveltime depends only on the ray parameter for a spherically symmetric earth. If we introduce an aspherical perturbation, including general ani-sotropy, the dependence on Euler angles can be expanded in terms of the rotation matrix for a fixed ray parameter. If the perturbation is isotropic, the expansion coefficients satisfy certain relations which may be used to obtain definite evidence for anisotropy rather than isotropic lateral heterogeneity.     

A global study of S-to-P and P-to-S conversions from the upper mantle transition zone

b
Long-period data of the Global Digital Seismograph Network (GDSN) recorded over the three-year period from 1984 to 1986 were studied for the occurrence of S-P and P-S conversions from the upper mantle transition zone that appear as precursors to teleseismic S arrivals. Conversions of this type were identified on a large number of single-station records. Simple stacking of many records enhanced the appearance of converted phases and demonstrated that no major lateral variations in the nature of the transition zone exist between various tectonic regions. S-P and P-S conversions from the 400 km discontinuity were best observed at distances between 70 and 85 while conversions from the 670 km discontinuity showed up best at distances beyond 87. The analysis of published source mechanisms and comparison with synthetic seismograms suggests that the appearance of converted phases is primarily governed by the earthquake radiation pattern. Phases that have undergone S-P conversions beneath the receiver are best observed from dip-slip events that radiate strong SV - and weak P -waves towards the station. P-S conversions beneath the source area, on the other hand, are frequently observed from events that radiate strong P and little SV energy towards the station, and also from some strike-slip events. Comparison of observed with synthetic seismograms suggests that the PREM model of Dziewonski & Anderson (1981) explains most of the observations. Observed S-P and P-S conversions from the 670 km discontinuity, however, often have larger amplitudes than in the synthetics. Constructive interference of converted waves with the P -wave coda, source radiation effects and a velocity contrast across the 670 km discontinuity which is higher than in PREM may all contribute to the discrepancy.     

Three-dimensional lithospheric modelling at NORSAR: linearity of the method and amplitude variations from the anomalies

Summary. A modification of the Aid et al . technique for three-dimensional lithospheric modelling is used to find smoothly varying models for the P -wave velocity structure beneath NORSAR. The method includes ray tracing and calculation of geometrical spreading in the anomalies. The results of linear inversion of the travel-time data compare well with those of previous investigators. The assumption of linearity, which removes the need to ray trace through the anomalies, is tested with iterative solutions for both synthetic and real data. A model with an rms velocity perturbation of 3 per cent, extending to 120 km depth, is found to be reasonably linear. In fact the procedure leads to two models which satisfy the same amount of the real data but which differ by far more than the standard errors. However, these differences are not significant once the imperfect resolution is accounted for by using the total estimation error of the stochastic inverse.
The depth of major anomalies appears to be greater than the array diameter and is therefore not well constrained. Comparing the geometrical spreading produced by these models with the amplitude variations observed at the array indicates that structure deeper than 120 km but shallower than 200 km makes an important contribution to the observations. None of the models used can produce variations as large as those in the amplitude data. For deep, essentially two-dimensional, anomalies the fit to these data is much better for sources to the NE of the array than for sources in other quadrants.     

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.     

Finite-frequency sensitivity of body waves to anisotropy based upon adjoint methods

We investigate the sensitivity of finite-frequency body-wave observables to mantle anisotropy based upon kernels calculated by combining adjoint methods and spectral-element modelling of seismic wave propagation. Anisotropy is described by 21 density-normalized elastic parameters naturally involved in asymptotic wave propagation in weakly anisotropic media. In a 1-D reference model, body-wave sensitivity to anisotropy is characterized by 'banana–doughnut' kernels which exhibit large, path-dependent variations and even sign changes. P -wave traveltimes appear much more sensitive to certain azimuthally anisotropic parameters than to the usual isotropic parameters, suggesting that isotropic P -wave tomography could be significantly biased by coherent anisotropic structures, such as slabs. Because of shear-wave splitting, the common cross-correlation traveltime anomaly is not an appropriate observable for S waves propagating in anisotropic media. We propose two new observables for shear waves. The first observable is a generalized cross-correlation traveltime anomaly, and the second a generalized 'splitting intensity'. Like P waves, S waves analysed based upon these observables are generally sensitive to a large number of the 21 anisotropic parameters and show significant path-dependent variations. The specific path-geometry of SKS waves results in favourable properties for imaging based upon the splitting intensity, because it is sensitive to a smaller number of anisotropic parameters, and the region which is sampled is mainly limited to the upper mantle beneath the receiver.     

2.5-dimensional crosshole acoustic response: a variable density approach

This paper presents the development of a 2.5-D simulation technique for acoustic wave propagation in media with variable density and velocity. A comparative study of the 2-D and 2.5-D responses of a model reveals the spatially and temporally damped nature of the 2.5-D acoustic wave equations. The simulated results for constant and variable density models show that the density variation affects only the reflectivity of the layer. The computational cost for variable density models is 2.17 and 2.26 times that for constant density models for the 2.5-D and 2-D cases, respectively. Furthermore, the 2.5-D computational cost in the time domain is only about 10–15 per cent more than that for two dimensions, so this modest increase in computational cost can avoid the exorbitant 3-D computational cost.
Snapshots for a crosshole geometry were computed at various times in order to study the effect of heterogeneity on the amplitude and shape of the wave front. Extensive analysis of an oil-bearing reservoir with and without the inclusion of a gas zone was performed using a point source as well as multiple sources. In addition, the effects of the thickness of a low-velocity layer (oil-bearing) and of the location of the source have been studied. It is concluded from the numerical response that the waveguide action of the low-velocity layer depends on its thickness in terms of the dominant wavelength. Trapping of waves was not observed when the source was outside the low-velocity layer. Furthermore, the presence of heterogeneity in the low-velocity layer contributes considerably to the leakage of energy in the adjacent layers due to scattering/diffraction. It was found that, in the 2.5-D numerical simulation, the stability condition and the requirement of the number of grid points per wavelength to avoid grid dispersion are the same as for the 2-D case.     

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.     

Three-dimensional seismic structure beneath the Australasian region from refracted wave observations

The earthquakes in the seismicity belt extending through Indonesia, New Guinea, Vanuatu and Fiji to the Tonga–Kermadec subduction zone recorded at the 65 portable broad-band stations deployed during the Skippy experiment from 1993–1996 provide good coverage of the lithosphere and mantle under the Australian continent, Coral Sea and Tasman Sea.
The variation in structure in the upper part of the mantle is characterized by deter-mining a suite of 1-D structures from stacked record sections utilizing clear P and S arrivals, prepared for all propagation paths lying within a 10° azimuth band. The azimuth of these bands is rotated by 20° steps with four parallel corridors for each azimuth. This gives 26 separate azimuthal corridors for which 15 independent 1-D seismic velocity structures have been derived, which show significant variation in P and S structure.
The set of 1-D structures is combined to produce a 3-D representation by projecting the velocity values along the ray path using a turning point approximation and stacking into 3-D cells (5° by 50 km in depth). Even though this procedure will tend to underestimate wave-speed perturbations, S -velocity deviations from the ak135 reference model exceed 6 per cent in the lithosphere.
In the uppermost mantle the results display complex features and very high S -wave speeds beneath the Precambrian shields with a significant low-velocity zone beneath. High velocities are also found towards the base of the transition zone, with high S -wave speeds beneath the continent and high P -wave speeds beneath the ocean. The wave-speed patterns agree well with independent surface wave studies and delay time tomography studies in the zones of common coverage.     

Simulating three-dimensional seismograms in 2.5-dimensional structures by combining two-dimensional finite difference modelling and ray tracing

Finite difference (FD) simulation of elastic wave propagation is an important tool in geophysical research. As large-scale 3-D simulations are only feasible on supercomputers or clusters, and even then the simulations are limited to long periods compared to the model size, 2-D FD simulations are widespread. Whereas in generally 3-D heterogeneous structures it is not possible to infer the correct amplitude and waveform from 2-D simulations, in 2.5-D heterogeneous structures some inferences are possible. In particular, Vidale & Helmberger developed an approach that simulates 3-D waveforms using 2-D FD experiments only. However, their method requires a special FD source implementation technique that is based on a source definition which is not any longer used in nowadays FD codes. In this paper, we derive a conversion between 2-D and 3-D Green tensors that allows us to simulate 3-D displacement seismograms using 2-D FD simulations and the actual ray path determined in the geometrical optic limit. We give the conversion for a source of a certain seismic moment that is implemented by incrementing the components of the stress tensor.
Therefore, we present a hybrid modelling procedure involving 2-D FD and kinematic ray-tracing techniques. The applicability is demonstrated by numerical experiments of elastic wave propagation for models of different complexity.     

Reflection—refraction of general P-and type-I S-waves in elastic and anelastic solids

Summary. The reflection and refraction of general (homogeneous or inhomo-geneous) plane P and type-I S ( SV ) body waves incident on plane boundaries are considered for general linear viscoelastic solids. Reflection—refraction laws, physical characteristics of the waves, and the nature of critical angles are examined in detail at welded boundaries and a free surface. General visco-elasticity with no low-loss approximations predicts that contrasts in intrinsic absorption at boundaries give rise to inhomogeneous reflected and refracted waves with elliptical particle motions, velocities and maximum attenuations that vary with frequency and angle of incidence, energy propagation at speeds and directions different from phase propagation, phase propagation that in general is parallel to the boundary for at most one angle of incidence, and reflection—transmission coefficients dependent on energy flow due to wave interaction. None of these physical characteristics are predicted for waves incident on boundaries that respond instantaneously.