Long-period corrections to body waves: theory

Summary. The usual asymptotic methods used to correct the high-frequency solutions of the wave equation are unsatisfactory as they do not give the low-frequency, partial reflections expected from a region of high velocity gradient. A new iterative solution is obtained which uses the first term of the Langer asymptotic expansion as the zeroth iterate. This satisfactorily gives the partial reflections from a region of high velocity gradient, even when they are generated near the turning point of the ray. Although the results are somewhat complicated in the frequency domain, in the time domain all types of wave interaction are described by six universal time functions. For any problem, these functions are scaled in time according to the depth of the interaction, and in strength according to the magnitude of the coupling parameter. Numerical results and approximations are given for these functions. Coupling parameters are investigated for acoustic and elastic waves in a plane model, and acoustic and elastic-gravitational waves in a spherical model. The same universal time functions allow the excitation of elastic waves to be studied when the source is in a region of high velocity gradient or is near the wave's turning point. Results are given for a moment tensor, point source in plane and spherical models.     

A new class of propagation models based on a factorization of the Helmholtz equation

Summary. A new class of propagation models, which is based on a factorization of the Helmholtz equation, that is, in principle, exact for range-independent environments, is presented. A numerical algorithm which makes use of the forward propagation feature of the propagation models and 'marches' the radiation field over a sequence of incremented range planes is derived and implemented. Numerical results of transmission loss as a function of range are obtained for a number of sample wave speed profiles for two members of the class of propagation models, one being an often-used parabolic wave theory and one being a high-frequency model that is new to our study. These results are compared to each other, to results obtained from an 'exact' numerical integration of the complete Helmholtz theory, and to results obtained using a WKB approximation to the Helmholtz equation formulation. The results demonstrate both the computational viability of the class of propagation models using the derived marching algorithm and very promising numerical agreement between the 'exact' results and those obtained using the high-frequency model.     

A variational approach to elastic wave propagation in anomalous layers

Summary . A unified variational approach is derived for elastic wave propagation in infinite and semi-infinite stratified media involving anomalous layers. As a specific application we provide accurate approximations to the transition matrix which appears in the so-called amplification problem in earthquake engineering.     

Relation between beachface morphology and wave climate at Trafalgar beach (Cádiz, Spain)

Two years of offshore wave data and daily time exposure images from Trafalgar beach, a 2-km-long sandy beach located on the southwest coast of Spain that frequently exhibits rhythmic features, were used to (1) explore the variability of the beachface morphology and (2) determine environmental conditions associated with the different morphological states. The beachface morphology at three distinct alongshore sectors was analyzed and classified and five different morphological states were found that are related with the presence or absence of beach cusps and a berm: (1) large beach cusps, (2) small beach cusps, (3) low-tide terrace; (4) plane beach berm and (5) plane beach. The predominant beachface morphology is characterized by the presence of large beach cusps, and the main wave climate consisted of offshore significant wave heights ranging from 0.5–1 m and wave periods between 4 and 12 s. An alongshore variation of the morphology is found which might be related to the nearshore wave variability (SWAN wave model results). The morphologies are, in some cases, well-correlated with the daily offshore incident wave climate (described by the daily maximum significant wave height and the corresponding period), particularly for the moderate to high energy wave conditions. Small beach cusps appear under short period waves, whereas when the wave periods are longer the morphology tends to change to large beach cusps. This transition only occurs if the forcing is maintained as constant for a certain duration, which depends itself on the wave energy. It is concluded that correlations over 90% are only found for the highest wave energy conditions or under long wave periods. For the remainder, it is not possible to generally correlate the beachface morphology based only on the wave forcing because the previous morphological state cannot be ignored.     

Teleseismic waveform modelling with a one-way wave equation

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

Exact and Approximate Generalized Ray Theory in Vertically Inhomogeneous Media

The generalized ray method in a vertically inhomogeneous model is formulated without any approximation by homogeneous layers. The solution is obtained as an infinite series in multiply 'reflected' waves. Each term can be solved using the exact method or the plane-wave, first-motion or geometrical approximations. It is shown that the first-motion approximation of the series converges rapidly, the ratio of successive terms in the infinite series being-(2 l + 1)(2 l )(6/π)².
In addition it is shown that the first-motion approximation, which reduces to the geometrical approximation when the latter is valid, is a useful alternative to geometrical ray theory, being more generally valid and being almost as simple to compute.     

Gaussian beams, complex rays, and the analytic extension of the Green's function in smoothly inhomogeneous media

Summary. Some relations between Gaussian beams, complex rays and the analytic extension of the Green's function in smoothly inhomogeneous media are shown in this paper. It is found that: (1) a single Gaussian beam is a paraxial approximation of the analytical extension of the ray-approximated Green's function in smoothly inhomogeneous media by putting the source point into a complex space. The Gaussian beam approximation of the Green's function has an advantage in computational efficiency and stability and can avoid the singularity problems at caustics, but also introduces a parabolic approximation to the wavefront and an angle-dependent amplitude damping. Therefore the validity of the Gaussian beam approximation should be checked using other methods. (2) Complex-ray tracing, which does not involve the paraxial approximation, can also avoid the singularity problemsm though without the computational efficiency. Therefore, it should be used to verify the Gaussian beam approximation, whenever possible. (3) The decomposition of a plane wave into an ensemble of Gaussian beams is equivalent to approximating the Green's function (the kernel of the ray-Kirchhoff method) with a single Gaussian beam. This introduces a parabolic approximation to the wavefront and a Gaussian windowing for arrival angles which may cause some problems in modelling wave propagation and scattering and has no advantages over other methods. (4) The representation of a point source field by a superposition of Gaussian beams, on the other hand, is equivalent to approximating the Green's function with a bundle of overlapped Gaussian beams. This representation is similar to a Maslov uniform asymptotic representation. It has no caustics and has improved accuracies at the caustics for quasi-plane waves compared to the extended WKBJ method.     

Regional tomographic inversion of the amplitude and phase of Rayleigh waves with 2-D sensitivity kernels

In this study, we test the adequacy of 2-D sensitivity kernels for fundamental-mode Rayleigh waves based on the single-scattering (Born) approximation to account for the effects of heterogeneous structure on the wavefield in a regional surface wave study. The calculated phase and amplitude data using the 2-D sensitivity kernels are compared to phase and amplitude data obtained from seismic waveforms synthesized by the pseudo-spectral method for plane Rayleigh waves propagating through heterogeneous structure. We find that the kernels can accurately predict the perturbation of the wavefield even when the size of anomaly is larger than one wavelength. The only exception is a systematic bias in the amplitude within the anomaly itself due to a site response.
An inversion method of surface wave tomography based on the sensitivity kernels is developed and applied to synthesized data obtained from a numerical simulation modelling Rayleigh wave propagation over checkerboard structure. By comparing recovered images to input structure, we illustrate that the method can almost completely recover anomalies within an array of stations when the size of the anomalies is larger than or close to one wavelength of the surface waves. Surface wave amplitude contains important information about Earth structure and should be inverted together with phase data in surface wave tomography.     

Seismogram synthesis using normal mode superposition: the locked mode approximation

Summary. A normal mode superposition approach is used to synthesize complete seismic codas for flat layered earth models and the P-SV phases. Only modes which have real eigenwavenumbers are used so that the search for eigenvalues in the complex wavenumber plane is confined to the real axis. In order to synthesize early P -wave arrivals by summing a number of'trapped'modes, an anomalously high velocity cap layer is added to the bottom of the structure so that most of the seismic energy is contained in the upper layers as high-order surface waves. Causality arguments are used to define time windows for which the resulting synthetic seismograms are close approximations to the exact solutions without the cap layer. The traditional Thomson—Haskell matrix approach to computing the normal modes is reformulated so that numerical problems encountered at high frequencies are avoided and numerical results of the locked mode approximation are given.     

Waveform inversion of mantle Love waves:the Born seismogram approach

Summary. Normal mode theory, extended to the slightly laterally heterogeneous earth by the first-order Born approximation, is applied to the waveform inversion of mantle Love wave (200–500 s) for the Earth's lateral heterogeneity at l = 2 and a spherically symmétric anelasticity ( Q _μ) structure. The data are from the Global Digital Seismograph Network (GDSN). The l =2 pattern is very similar to the results of other studies that used either different méthods, such as phase velocity measurements and multiplet location measurements, or a different data set, such as mantle Rayleigh waves from different instruments. The results are carefully analysed for variance reduction and are most naturally explained by heterogeneity in the upper 420 km. Because of the poor resolution of the data set for the deep interior, however, a fairly large heterogeneity in the transition zones, of the order of up to 3.5 per cent in shear wave velocity, is allowed. It is noteworthy that Love waves of this period range cannot constrain the structure below 420 km and thus any model presented by similar studies below this depth are likely to be constrained by Rayleigh waves (spheroidal modes) only.
The calculated modal Q values for the obtained Q _μ model fall within the error bars of the observations. The result demonstrates the discrepancy of Rayleigh wave Q and Love wave Q and indicates that care must be taken when both Rayleigh and Love wave data, including amplitude information, are inverted simultaneously.
Anomalous amplitude inversions of G2 and G3, for example, are observed for some source-receiver pairs. This is due to multipathing effects. One example near the epicentral region, which is modelled by the obtained l = 2 heterogeneity, is shown.     

Receiver functions at the stations of the German Regional Seismic Network (GRSN)

Several years of broad-band teleseismic data from the GRSN stations have been analysed for crustal structure using P -to- S converted waves at the crustal discontinuities. An inversion technique was developed which applies the Thomson-Haskell formalism for plane waves without slowness integration. The main phases observed are Moho conversions, their multiples in the crust, and conversions at the base of the sediments. The crustal thickness derived from these data is in good agreement with results from other studies. For the Gräfenberg stations, we have made a more detailed comparison of our model with a previously published model obtained from refraction seismic experiments. The refraction seismic model contains boundaries with strong velocity contrasts and a significant low-velocity zone, resulting in teleseismic waveforms that are too complicated as compared to the observed simple waveforms. The comparison suggests that a significant low-velocity zone is not required and that internal crustal boundaries are rather smooth.     

Seismic reflection traveltimes in two-dimensional statistically anisotropic random media

Velocity estimation remains one of the main problems when imaging the subsurface with seismic reflection data. Traveltime inversion enables us to obtain large-scale structures of the velocity field and the position of seismic reflectors. However, as the media currently under study are becoming more and more complex, we need to know the finer-scale structures. The problem is that below a certain range of velocity heterogeneities, deterministic methods become difficult to use, so we turn to a probabilistic approach. With this in view, we characterize the velocity field as a random field defined by its first and second statistical moments. Usually, a seismic random medium is defined as a homogeneous velocity background perturbed by a small random field that is assumed to be stationary. Thus, we make a link between such a random velocity medium (together with a simple reflector) and seismic reflection traveltimes. Assuming that the traveltimes are ergodic, we use 2-D seismic reflection geometry to study the decrease in the statistical traveltime fluctuations as a function of the offset (the source–receiver distance). Our formulae are based on the Rytov approximation and the parabolic approximation for acoustic waves. The validity and the limits are established for both of these approximations in statistically anisotropic random media. Finally, theoretical inversion procedures are developed for the horizontal correlation structure of the velocity heterogeneities for the simplest case of a horizontal reflector. Synthetic seismograms are then computed (on particular realizations of random media) by simulating scalar wave propagation via finite difference algorithms. There is good agreement between the theoretical and experimental results.     

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.     

Partial derivatives of surface wave phase velocities for flat anisotropic models

Summary. Surface wave behaviour in flat anisotropic structures is first illustrated by performing an exact computation on a simple two-layer model. The variational procedure of Smith & Dahlen is then used to compute the partial derivatives of surface wave phase velocities with respect to the elastic parameters in more realistic earth models. Linear relationships between the partial derivatives for a general anisotropic structure and those for a transversely isotropic structure are derived. When considering waves propagating in a fixed direction, there are only four independent derivatives for Rayleigh waves, and two for Love waves. To avoid the lack of resolution in an inverse method, we propose to use physically constrained models. These results are illustrated by using a model with hexagonal symmetry and a symmetry axis oriented either vertically or horizontally. Quasi-Love- and quasi-Rayleigh-wave partial derivatives are computed for both axis orientations. Modes up to the second overtone and periods ranging between 45 and 130 s have been considered. Finally, anomalies of phase velocity are computed in an oceanic model made of 1/6 oriented olivine crystals with horizontal or vertical preferred orientations of the a -axis.     

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).     

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.     

The subseismic approximation in core dynamics—II. Love numbers and surface gravity

This paper extends our earlier examinations of the utility of various approximations for treating the dynamics of the Earth's liquid core on time-scales of the order of 10⁴ to 10⁸ s. We discuss the effects of representing the response of the mantle and inner core by static (versus dynamic) Love numbers, and of invoking the subseismic approximation for treating core flow, used either only in the interior of the liquid core (SSA-1) or also at the boundaries (SSA-2). The success of each approximation (or combinations thereof) is measured by comparing the resulting surface gravity effects (computed for a given earthquake excitation), and (for the Slichter mode) the distribution of translational momentum, with reference calculations in which none of these approximations is made. We conclude that for calculations of the Slichter triplet, none of the approximations is satisfactory, i.e. a full solution (using dynamic Love numbers at elastic boundaries and no core flow approximation) is required in order to avoid spurious eigenfrequencies and to yield correct eigenfunctions (e.g. conserving translational momentum) and surface gravity. For core undertones, the use of static Love numbers at rigid boundaries is acceptable, along with SSA-1 (i.e. provided the subseismic approximation is not invoked at the core boundaries). Although the calculations presented here are for a non-rotating earth model, we argue that the principal conclusions should be applicable to the rotating Earth. Shortcomings of the subseismic approximation appear to arise because both SSA-1 and SSA-2 lower the order of the governing system of differential equations (giving rise to a singular perturbation problem), and because SSA-2 overdetermines the boundary conditions (making it impossible for solutions to satisfy all continuity requirements at core boundaries).     

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.     

Approximation of Gaussian spatial autoregressive models for massive regular square tessellation data

Most of the literature to date proposes approximations to the determinant of a positive definite n × n spatial covariance matrix (the Jacobian term) for Gaussian spatial autoregressive models that fail to support the analysis of massive georeferenced data sets. This paper briefly surveys this literature, recalls and refines much simpler Jacobian approximations, presents selected eigenvalue estimation techniques, summarizes validation results (for estimated eigenvalues, Jacobian approximations, and estimation of a spatial autocorrelation parameter), and illustrates the estimation of the spatial autocorrelation parameter in a spatial autoregressive model specification for cases as large as n = 37,214,101. The principal contribution of this paper is to the implementation of spatial autoregressive model specifications for any size of georeferenced data set. Its specific additions to the literature include (1) new, more efficient estimation algorithms; (2) an approximation of the Jacobian term for remotely sensed data forming incomplete rectangular regions; (3) issues of inference; and (4) timing results.