Phase–velocity dispersion curves and small-scale geophysics using noise correlation slantstack technique

It has been demonstrated both theoretically and experimentally that the Green's function between two receivers can be retrieved from the cross-correlation of isotropic noise records. Since surface waves dominate noise records in geophysics, tomographic inversion using noise correlation techniques have been performed from Rayleigh waves so far. However, very few numerical studies implying surface waves have been conducted to confirm the extraction of the true dispersion curves from noise correlation in a complicated soil structure. In this paper, synthetic noise has been generated in a small-scale (<1 km) numerical realistic environment and classical processing techniques are applied to retrieve the phase velocity dispersion curves, first step toward an inversion. We compare results obtained from spatial autocorrelation method (SPAC), high-resolution frequency-wavenumber method (HRFK) and noise correlation slantstack techniques on a 10-sensor array. Two cases are presented in the (1–20 Hz) frequency band that corresponds to an isotropic or a directional noise wavefield. Results show that noise correlation slantstack provides very accurate phase velocity estimates of Rayleigh waves within a wider frequency band than classical techniques and is also suitable for accurately retrieving Love waves dispersion curves.     

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.     

Rayleigh-wave dispersion function for a transversely isotropic layered spherical earth using the Thomson-Haskell matrix method

We consider the two coupled differential equations of the two radial functions appearing in the displacement components of spheroidal oscillations for a transversely isotropic (TI) medium in spherical coordinates. Elements of the layer matrix have been explicitly written—perhaps for the first time—to extend the use of the Thomson-Haskell matrix method to the derivation of the dispersion function of Rayleigh waves in a transversely isotropic spherical layered earth. Furthermore, an earth-flattening transformation (EFT) is found and effectively used for spheroidal oscillations. The exponential function solutions obtained for each layer give the dispersion function for TI spherical media the same form as that on a flat earth. This has been achieved by assuming that the five elastic parameters involved vary as r ^p and that the density varies as r ^p-2, where p is an arbitrary constant and r is the radial distance. A numerical illustration with p = - 2 shows that, in spite of the inhomogeneity assumed within layers, the results for spherical harmonic degree n , versus time period T , obtained here for the Primary Reference Earth Model (PREM), agree well with those obtained earlier by other authors using numerical integration or variational methods. The results for isotropic media derived here are also in agreement with previous results. The effect of transverse isotropy on phase velocity for the first two modes of Rayleigh waves in the period range 20 to 240 s is calculated and discussed for continental and oceanic models.     

The explicit crustal transfer function for SH-waves by ray theory — the case of one and two layers

Summary. Previous theoretical studies on the transfer function of crustal plane surface layers have been primarily based on applying recursion formulae to the solutions of wave equations. In this way the detailed physical insight of the transfer function is often obscure and it is very difficult to express the transfer function in an explicit form. To show that this situation can be improved by an alternative approach, we demonstrate in this paper with models of one and two surface layers that the explicit transfer function can be derived by summing multiple-reflected rays between interfaces. Since the derived transfer function is related explicitly to some physical parameters, such as wave travel times between interfaces, reflection and transmission coefficients at the interfaces, and attenuation and dispersion of waves in each layer of the model, it is more flexible than the traditional recursion form when applied to different models. This suggests that the ray theory can play an important role in problems of layered media.     

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.     

Elastic wave propagation in model sediments – II

Summary. A fluid-saturated cubic packing of like elastic spheres is taken to be in equilibrium under the effect of gravity and the effects of a superimposed low-frequency elastic wave are considered. In the first place, expressions for the wave velocity, dispersion and attenuation are derived for the dry packing. This dynamic theory leads to the result that, for very low frequencies, the wave velocity is proportional to the third root of the depth and not the sixth root as is obtained by using the effective elastostatic modulus of the packing. For the fluid-saturated packing, two waves, termed respectively the 'solid wave' and the 'fluid wave', are found to propagate. The 'solid wave' has the characteristics of a wave propagating within a dry packing whose parameters differ in a specified way from those of the original packing, whereas the 'fluid wave' has those of a wave within a homogeneous fluid with similarly modified parameters.     

Gaussian beams for surface waves in laterally slowly-varying media

Summary. Asymptotic ray theory is applied to surface waves in a medium where the lateral variations of structure are very smooth. Using ray-centred coordinates, parabolic equations are obtained for lateral variations while vertical structural variations at a given point are specified by eigenfunctions of normal mode theory as for the laterally homogeneous case. Final results on wavefields close to a ray can be expressed by formulations similar to those for elastic body waves in 2-D laterally heterogeneous media, except that the vertical dependence is described by eigenfunctions of 'local' Love or Rayleigh waves. The transport equation is written in terms of geometrical-ray spreading, group velocity and an energy integral. For the horizontal components there are both principal and additional components to describe the curvature of rays along the surface, as in the case of elastic body waves. The vertical component is decoupled from the horizontal components. With complex parameters the solutions for the dynamic ray tracing system correspond to Gaussian beams: the amplitude distribution is bell-shaped along the direction perpendicular to the ray and the solution is regular everywhere, even at caustics. Most of the characteristics of Gaussian beams for 2-D elastic body waves are also applicable to the surface wave case. At each frequency the solution may be regarded as a set of eigenfunctions propagating over a 2-D surface according to the phase velocity mapping.     

Inferring the relative measure of principal stress components

The phase velocity and the attenuation coefficient of compressional seismic waves, propagating in poroelastic, fluid-saturated, laminated sediments, are computed analytically from first principles. The wavefield is found to be strongly affected by the medium heterogeneity. Impedance fluctuations lead to poroelastic scattering; variations of the layer compressibilities cause inter-layer flow (a 1-D macroscopic local flow). These effects result in significant attenuation and dispersion of the seismic wavefield, even in the surface seismic frequency range, 10–100 Hz. The various attenuation mechanisms are found to be approximately additive, dominated by inter-layer flow at very low frequencies. Elastic scattering is important over a broad frequency range from seismic to sonic frequencies. Biot's global flow (the relative displacement of solid frame and fluid) contributes mainly in the range of ultrasonic frequencies. From the seismic frequency range up to ultrasonic frequencies, attenuation due to heterogeneity is strongly enhanced compared to homogeneous Biot models. Simple analytical expressions for the P -wave phase velocity and attenuation coefficient are presented as functions of frequency and of statistical medium parameters (correlation lengths, variances). These results automatically include different asymptotic approximations, such as poroelastic Backus averaging in the quasi-static and the no-flow limits, geometrical optics, and intermediate frequency ranges.     

Calculation of the elasto-dynamic Green's function in layered media by means of a modified propagator matrix method

Summary. The calculation of the two-dimensional elasto-dynamic Green's function for a stratified medium is investigated. The solution is represented in the form of an inverse Fourier integral which is to be integrated along a properly chosen path in the complex wavenumber plane. The integrand is computed using a modified propagator matrix method.
This method is based on a mixed formulation using the propagator matrix and the matrix of minors of the propagator matrix (compound matrix). The major advantages of this approach are the elimination of the numerical loss of precision problems associated with the Thomson-Haskell formulation, without losing the attractive tractability and compactness of the propagator matrix method. This modified method is first mathematically derived, and theoretical seismograms are then presented for two examples.     

Multiplet-Multiplet Coupling Due to Lateral Heterogeneity: Asymptotic Effects On the Amplitude and Frequency of the Earth's Normal Modes

Summary. Starting with the first-order formulation of quasi-degenerate splitting theory for the normal modes of a laterally heterogeneous earth, we have obtained an asymptotic expression for the coupling terms corresponding to neighbouring multiplets along the same dispersion branch as the mode considered, valid to order 1/ℓ, where ℓ is the angular order of this mode (ℓ≥ 1).
We show that, to order zero, these coupling terms introduce a small shift in epicentral distance into the expression for the long period seismogram obtained by normal mode summation. This shift depends on the difference between the great circle and the minor arc averages of the local frequency. the coupling terms thus permit us to reconcile results obtained by normal-mode summation and by a propagating wave approach, as far as the dependence on structure of the phase of surface waves is concerned.
To order 1/ℓ, the coupling terms result in a perturbation in the amplitude of the mode considered, which depends on spatial derivatives of the local frequency and thus on the structure in the vicinity of the source station great circle path. We show that this term is equivalent to that which is found using ray perturbation methods for propagating surface waves. We compare and discuss the assumptions underlying both approaches and illustrate, by an example, the potential of the asymptotic normal-mode formulation for improved modelling of lateral heterogeneity in the earth.     

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.     

Global observation of off-great-circle propagation of Long-Period surface waves

In the current generation of global dispersion maps of surface waves, the long-wavelength structure seems to be very well determined. There is general agreement in the patterns of global phase velocity anomalies up to harmonic degree 16. However, the shorter-wavelength structure varies significantly between published maps, and it appears that this part of the models depends strongly on the inversion technique and on the data set of surface-wave dispersion (usually phase measurements). Polarization data depend on the lateral gradient of phase velocity and hence are more sensitive to shorter-wavelength structure than phase data; thus, including these data should enhance resolution. In this paper, I demonstrate that polarization data of long-period surface waves (80 s), as a function of frequency, can be reliably measured using a multitaper technique. the resulting off-great-circle arrival angles of the surface-wave packets are relatively easy to interpret within a ray-theoretical framework. Our data base of three-component recordings is now large enough to provide useful constraints on global dispersion maps, particularly on the shorter-wavelength parts. Apart from the phase velocity model itself, a possible misorientation of the horizontal components at each station is included in a non-linear inversion as an additional independent model parameter. This gives a significant improvement in the fit to the data. Misorientations of more than 3° are probable for at least four of the 37 stations investigated.     

A search for the oceanic Lg phase

Summary. The Lg phase has been shown previously to be a collection of higher-mode surface waves guided by the continental crust (Knopoff, Schwab & Kausel). A simple scaling between continental and oceanic crustal thicknesses suggests that a search for an oceanic Lg phase should be made in the period range from 1 to 2s. In a search for SH polarized Lg arrivals over oceanic paths, we found that in addition to the fundamental mode, seismo-grams at relatively short ranges in the Pacific showed the presence of only the first higher mode with group velocities on the steep portion of the dispersion curve rather than at the group velocity minimum as expected. Numerical model analysis indicates that, contrary to the continental case, there is no strong confluence of stationary phases of higher-mode crustal waves in the appropriate period range to produce Lg wave packets; this is due to small but significant differences in scaled crustal structures. Further, lateral variations in the thickness of oceanic sediments are sufficient to scatter most of the crustal surface-wave energy within a relatively short distance. Even were this thickness uniform, attenuation in the sediments would be strong enough to absorb the Lg stationary phases in a short distance.     

Seismic vertical array analysis for phase decomposition

We propose a vertical array analysis method that decomposes complex seismograms into body and surface wave time histories by using a velocity structure at the vertical array site. We assume that the vertical array records are the sum of vertically incident plane P and S waves, and laterally incident Love and Rayleigh waves. Each phase at the surface is related to that at a certain depth by the transfer function in the frequency domain; the transfer function is obtained by Haskell's matrix method, assuming a 1-D velocity structure. Decomposed P , S and surface waves at the surface are estimated from the vertical array records and the transfer functions by using a least-squares method in the frequency domain; their time histories are obtained by the inverse Fourier transform. We carried out numerical tests of this method based on synthetic vertical array records consisting of vertically incident plane P and S waves and laterally incident plane Love and Rayleigh waves. Perfect results of the decomposed P , S , Love and Rayleigh waves were obtained for synthetic records without noise. A test of the synthetic records in which a small amount of white noise was added yielded a reasonable result for the decomposed P , S and surface waves. We applied this method to real vertical array records from the Ashigara valley, a moderate-sized sedimentary valley. The array records from two earthquakes occurring at depths of 123 and 148 km near the array (epicentral distance of about 31 km) exhibited long-duration later phases. The analysis showed that duration of the decomposed S waves was a few seconds and that the decomposed surface waves appeared a few seconds after the direct S -wave arrival and had very long duration. This result indicated that the long-duration later phases were generated not by multireflected S waves, but by basin-induced surface waves.     

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.     

On the inversion of Love- and Rayleigh-wave dispersion and implications for Earth structure and anisotropy

Summary. Various factors can make it difficult to explain observations of Love- and Rayleigh-wave dispersion with the same relatively simple isotropic model. These factors include systematic errors which might occur in determinations of observed group and phase velocities, lateral variations in structure along the path of travel, and the attempt to explain observations with a model comprised of only a small number of thick layers. The last of these factors is illustrated by an inversion of dispersion data in the central United States where shear-wave anisotropy had previously been invoked as one way to explain incompatible Love- and Rayleigh-wave velocities. It is shown that the data can be satisfied equally well by an isotropic model consisting of several thin layers.
In cases where the incompatibility of Love- and Rayleigh-wave data might be produced by intrinsic anisotropy, it is necessary to invert those data using an anisotropic theory rather than by separate isotropic inversions of Love and Rayleigh waves. Inversions of fundamental-mode data for a region of the Pacific, assuming anisotropic media in which the layers are transversely isotropic with a vertical axis of symmetry, lead to models which are highly non-unique. Even if the inversions solve only for shear velocities in the litho-sphere and asthenosphere it is not possible, without supplementary information, to ascertain the depth interval over which anisotropy occurs or to determine the thickness of the lithosphere or asthenosphere with much precision.     

The Structure of Atmospheric Turbulence and its Application to the Design of Pipe Arrays

Summary Turbulent boundary layers at the surface of the Earth limit the detection of infrasonic waves with periods greater than 1 s. Pipe arrays designed to improve the signal-to-noise ratios of infrasonic waves usually assume that the background noise due to this turbulent boundary layer is incoherent between the array inlets. The power at various points on a surface was measured; coherences between these points were determined and they were found to be significant in the period range 1–100 s. Such coherent noise must be considered when pipe arrays are designed.     

RAYLEIGH WAVES IN A MULTI-LAYERED MEDIUM WITH APPLICATIONS TO MICROSEISMS

Summary. An analysis is presented for the propagation of Rayleigh waves in a multi-layered medium with water as the uppermost layer from which it is possible to obtain the group velocity as well as the phase velocity for the Rayleigh waves in an oceanic path with the ocean bed consisting of two upper layers each of uniform composition and finite thickness, underlain by ultrabasic rock extended to large depths. The results have been applied in the analysis of microseisms associated with cyclones and "norwesters". It has been shown that the period of the microseisms recorded at different stations depends critically on the structure of the continental shelf. The results obtained in the present paper, when compared with the nature of the recorded microseisms associated with cyclones, and their variations, give useful information about the depth and structure of the ocean beds.     

Implementation of perfectly matched layers in an arbitrary geometrical boundary for elastic wave modelling

The perfectly matched layer (PML) absorbing boundary condition is incorporated into an irregular-grid elastic-wave modelling scheme, thus resulting in an irregular-grid PML method. We develop the irregular-grid PML method using the local coordinate system based PML splitting equations and integral formulation of the PML equations. The irregular-grid PML method is implemented under a discretization of triangular grid cells, which has the ability to absorb incident waves in arbitrary directions. This allows the PML absorbing layer to be imposed along arbitrary geometrical boundaries. As a result, the computational domain can be constructed with smaller nodes, for instance, to represent the 2-D half-space by a semi-circle rather than a rectangle. By using a smooth artificial boundary, the irregular-grid PML method can also avoid the special treatments to the corners, which lead to complex computer implementations in the conventional PML method. We implement the irregular-grid PML method in both 2-D elastic isotropic and anisotropic media. The numerical simulations of a VTI lamb's problem, wave propagation in an isotropic elastic medium with curved surface and in a TTI medium demonstrate the good behaviour of the irregular-grid PML method.