Windowed f-k spectra of a three-dimensional wavefield excited by a point source in a two-dimensional multilayered elastic medium

When full 3-D modelling is too costly or cumbersome, computations of 3-D elastic wave propagation in laterally heterogeneous, multilayered 2-D geological structures may enhance considerably our ability to predict strong ground motion for seismological and engineering purposes. Towards this goal, we extend the method based on the combination of the thin-layer finite-element and boundary-element methods (TLFE-BEM) and calculate windowend f - k spectra of the 3-D wavefield. The windowed f - k spectra are spatially localized spectra from which the local properties of the wavefield can be extracted. The TLFE-BEM is particularly suited for calculating the complete wavefield where surface waves are dominant in multilayered media. The computations are performed in the frequency domain, providing the f - k spectra directly. From the results for the 3-D wavefield excited by a point source in a 2-D multilayered, sloped structure, it can be said that the phase velocity of the fundamental-mode Rayleigh wave in a laterally heterogeneous multilayered medium, estimated from the windowed f - k spectra, varies with the location of the point source. For the model calculated in this article, the phase velocity varies between the value for the flat layered structure of the thick-layer side and that for the structure just under the centre of the window. The exact subsurface structure just under the centre of an array in a laterally heterogeneous medium cannot be obtained if we use the f - k spectral analysis assuming a flat layered structure.     

Two efficient algorithms for iterative linearized inversion of seismic waveform data

We present two equivalent algorithms for iterative linearized waveform inversion for 3-D Earth structure with respect to an arbitrary 3-D starting model; one is a matrix formulation, and the second is a wavefield formulation. Both algorithms require the computation of accurate synthetic seismograms, but neither requires that any particular method be used to compute the synthetics. The matrix formulation is equivalent to our previously published algorithm (Hara, Tsuboi & Geller 1991), but requires less than 10 per cent of the CPU time of the previous algorithm. The wavefield algorithm is equivalent to that of Tarantola (1986) and Mora (1987), but appears to be substantially more efficient.     

Inversion for laterally heterogeneous upper mantle S-wave velocity structure using iterative waveform inversion

In contrast to previous work, which treats the Earth's lateral heterogeneity as an infinitesimal perturbation to a spherically symmetrical starting model, we conduct iterative linearized waveform inversion for the Earth's laterally heterogeneous structure. We use the Direct Solution Method (DSM) (Geller et al. 1990a) to calculate synthetic seismograms and their partial derivatives for a laterally heterogeneous earth model. We invert surface-wave data from the IDA and GEOSCOPE networks. We expand the lateral heterogeneity of rigidity in spherical harmonics up to angular order number 8 and use three parameters to specify the depth dependence of each harmonic, giving us a total of 240 unknowns. The short-wavelength lateral heterogeneity (s = 4, 6 and 8) in the deeper part of the upper mantle obtained by this study differs significantly from M84A. The relative improvement in the variance reduction as compared with model M84A is about 20 per cent for the IDA data and more than 100 per cent for the GEOSCOPE data.     

Inversion of seismic refraction data in planar dipping structure

Summary. A new method is presented for the direct inversion of seismic refraction data in dipping planar structure. Three recording geometries, each consisting of two common-shot profiles, are considered: reversed, split, and roll-along profiles. Inversion is achieved via slant stacking the common-shot wavefield to obtain a delay time–slowness (tau– p ) wavefield. The tau– p curves from two shotpoints describing the critical raypath of refracted and post-critically reflected arrivals are automatically picked using coherency measurements and the two curves are jointly used to calculate velocity and dip of isovelocity lines iteratively, thereby obtaining the final two-dimensional velocity model.
This procedure has been successfully applied to synthetic seismograms calculated for a dipping structure and to field data from central California. The results indicate that direct inversion of closely-spaced refraction/wide-aperture reflection data can practically be achieved in laterally inhomogeneous structures.     

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.     

A coupled mode examination of irregular waveguides including the continuum spectrum

Summary. A coupled mode theory is used to examine surface wave propagation in a laterally inhomogeneous acoustic waveguide. The theory is developed from the equations of motion for the pressure and velocity fields. The presence of lateral inhomogeneities in the form of varying layer thickness causes coupling among the discrete modes of the waveguide and radiation to the continuum. Expressions for the coupling coefficients among all mode types including coupling to the continuum spectrum are derived. The coupling coefficients are proportional to the horizontal derivative of the function describing the interface between layers of constant material properties but varying thickness. The coupled mode equations are solved in approximation for the case of a sinusoidal boundary and a sloping boundary. The results for radiation losses due to interaction with the irregular boundary of the waveguide are presented in analytical form, which clearly show the primary physical effects on the wavefield of the interaction. The far field amplitude of the scattered modes, excited by the interaction of some incident signal with a weak boundary irregularity, is modulated by the spatial Fourier transform of the irregularity.     

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.     

Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis – II. Crustal and upper-mantle structure

We determine the 3-D shear wave speed variations in the crust and upper mantle in the southeastern borderland of the Tibetan Plateau, SW China, with data from 25 temporary broad-band stations and one permanent station. Interstation Rayleigh wave (phase velocity) dispersion curves were obtained at periods from 10 to 50 s from empirical Green's function (EGF) derived from (ambient noise) interferometry and from 20 to 150 s from traditional two-station (TS) analysis. Here, we use these measurements to construct phase velocity maps (from 10 to 150 s, using the average interstation dispersion from the EGF and TS methods between 20 and 50 s) and estimate from them (with the Neighbourhood Algorithm) the 3-D wave speed variations and their uncertainty. The crust structure, parametrized in three layers, can be well resolved with a horizontal resolution about of 100 km or less. Because of the possible effect of mechanically weak layers on regional deformation, of particular interest is the existence and geometry of low (shear) velocity layers (LVLs). In some regions prominent LVLs occur in the middle crust, in others they may appear in the lower crust. In some cases the lateral transition of shear wave speed coincides with major fault zones. The spatial variation in strength and depth of crustal LVLs suggests that the 3-D geometry of weak layers is complex and that unhindered crustal flow over large regions may not occur. Consideration of such complexity may be the key to a better understanding of relative block motion and patterns of seismicity.     

Scattering parameters of the lithosphere below the Massif Central, France, from teleseismic wavefield records

Teleseismic P -wave recordings are analysed in the frequency range 0.3–6  Hz to derive structural (statistical) parameters of the lithosphere underneath the French Massif Central. For this we analyse differences in frequency-dependent intensities of the mean wavefield and the fluctuation wavefield. It is possible to discriminate a weak fluctuation regime of the wavefield in the frequency range below 1  Hz and a strong fluctuation regime starting above 1  Hz and continuing to higher frequencies. The observed wavefield fluctuations in the frequency range 0.3–3  Hz can be explained by scattering of the teleseismic P wave front at elastic inhomogeneities in the lithosphere. A statistical distribution of the inhomogeneities is assumed and the concept of random media is applied. The lithospheric structure under the Massif Central can be described as a 70  km thick heterogeneous layer with velocity fluctuations of 3–7 per cent and correlation lengths of the heterogeneities of 1–16  km.     

Crustal structure of Iraq from receiver functions and surface wave dispersion: implications for understanding the deformation history of the Arabian–Eurasian collision

We report the crustal structure for two locations in Iraq estimated by joint inversion of P -wave receiver functions (RFs) and surface (Rayleigh) wave group velocity dispersion. RFs were computed from teleseismic recordings at two temporary broad-band seismic stations located in Mosul (MSL) in the Zagros Fold Belt and Baghdad (BHD) in the Mesopotamian Foredeep. Group velocity dispersion curves at the sites were derived from continental-scale tomography. The inversion results show that the crustal thicknesses are 39 km at MSL and 43 km at BHD. We observe a strong Ps _Moho at BHD consistent with a sharp Moho discontinuity. However, at MSL we observe a weak Ps _Moho suggesting a transitional Moho where crustal thickening is likely to be occurring in the deep crust. Both sites reveal low velocity surface layers consistent with sedimentary thickness of about 3 km at station MSL and 7 km at BHD and agreeing well with the previous reports. Ignoring the sediments, the crystalline crustal velocities and thicknesses are remarkably similar at both stations. The similarity of crustal structure suggests that the crust of the northeastern proto-Arabian Platform was uniform before subsidence and deposition of the sediments in the Cenozoic. If crystalline crustal structure is uniform across the northern Arabian Platform then crustal thickness variations in the Zagros Fold Belt and Thrust Zone should reveal the history of deformation and crustal shortening in the Arabian–Eurasian collision zone and not reflect pre-existing crustal thickness variations in the Arabian Plate.     

Guided waves in three-dimensional structures

A new formulation for the propagation of surface waves in three-dimensionally varying media is developed in terms of modal interactions. A variety of assumptions can be made about the nature of the modal field: a single set of reference modes, a set of local modes for the structure beneath a point, or a set of local modes for a laterally varying reference structure. Each modal contribution is represented locally as a spectrum of plane waves propagating in different directions in the horizontal plane. The influence of 3-D structure is included by allowing coupling between different modal branches and propagation directions. For anisotropic models, with allowance for attenuation, the treatment leads to a set of coupled 2-D partial differential equations for the weight functions for different modal orders.
The representation of the guided wavefield requires the inclusion of a full set of modes, so that, even for isotropic models, both Love and Rayleigh modes appear as different polarization states of the modal spectrum. The coupling equations describe the interaction between the different polarizations induced by the presence of the 3-D structure.
The level of lateral variation within the 3-D model is not required to be small. Horizontal refraction or reflection of the surface wavefield can be included by allowing for transfer between modes travelling in different directions. Approximate forms of the coupled equation system can be employed when the level of heterogeneity is small, for example the coupling between the fundamental mode and higher modes can often be neglected, or forward propagation can be emphasized by restricting the interaction to a limited band of plane waves covering the expected direction of propagation.     

Finite-frequency traveltime tomography of San Francisco Bay region crustal velocity structure

Seismic velocity structure of the San Francisco Bay region crust is derived using measurements of finite-frequency traveltimes. A total of 57 801 relative traveltimes are measured by cross-correlation over the frequency range 0.5–1.5 Hz. From these are derived 4862 'summary' traveltimes, which are used to derive 3-D P -wave velocity structure over a 341 × 140 km² area from the surface to 25 km depth. The seismic tomography is based on sensitivity kernels calculated on a spherically symmetric reference model. Robust elements of the derived P -wave velocity structure are: a pronounced velocity contrast across the San Andreas fault in the south Bay region (west side faster); a moderate velocity contrast across the Hayward fault (west side faster); moderately low velocity crust around the Quien Sabe volcanic field and the Sacramento River delta; very low velocity crust around Lake Berryessa. These features are generally explicable with surface rock types being extrapolated to depth ∼10 km in the upper crust. Generally high mid-lower crust velocity and high inferred Poisson's ratio suggest a mafic lower crust.     

Computation of Large Anisotropic Seismic Heterogeneities (CLASH)

A general tomographic technique is designed in order (i) to operate in anisotropic media; (ii) to account for the uneven seismic sampling and (iii) to handle massive data sets in a reasonable computing time. One modus operandi to compute a 3-D body wave velocity model relies on surface wave phase velocity measurements. An intermediate step, shared by other approaches, consists in translating, for each period of a given mode branch, the phase velocities integrated along ray paths into local velocity perturbations. To this end, we develop a method, which accounts for the azimuthal anisotropy in its comprehensive form. The weakly non-linear forward problem allows to use a conjugate gradient optimization. The Earth's surface is regularly discretized and the partial derivatives are assigned to the individual grid points. Possible lack of lateral resolution, due to the inescapable uneven ray path coverage, is taken into account through the a priori covariances on parameters with laterally variable correlation lengths. This method allows to efficiently separate the 2ψ and the 4ψ anisotropic effects from the isotropic perturbations. Fundamental mode and overtone phase velocity maps, derived with real Rayleigh wave data sets, are presented and compared with previous maps. The isotropic models concur well with the results of Trampert & Woodhouse. Large 4ψ heterogeneities are located in the tectonically active regions and over the continental lithospheres such as North America, Antarctica or Australia. At various periods, a significant 4ψ signature is correlated with the Hawaii hotspot track. Finally, concurring with the conclusions of Trampert & Woodhouse, our phase velocity maps show that Rayleigh wave data sets do need both 2ψ and 4ψ anisotropic terms.     

Two-dimensional PP/PS-stereotomography: P- and S-waves velocities estimation from OBC data

We present the extension of stereotomography to P - and S -wave velocity estimation from PP - and PS -reflected/diffracted waves. In this new context, we greatly benefit from the use of locally coherent events by stereotomography. In particular, when applied to S -wave velocity estimation from PS -data, no pairing of PP - and PS -events is a priori required. In our procedure the P -wave velocity model is obtained first using stereotomography on PP -arrivals. Then the S -wave velocity model is obtained using PS -stereotomography on PS -arrivals fixing the P -wave velocity model. We present an application to an 'ideal' synthetic data set demonstrating the relevance of the approach, which allows us to recover depth consistent P - and S -waves velocity models even if no pairing of PP - and PS -events is introduced. Finally, results to a real data set from the Gulf of Mexico are presented demonstrating the potential of the method in a noisy data context.     

From ACH tomographic models to absolute velocity models

The ACH method, a widely used tomographic inverse method, is characterized by the use of relative residuals in order to avoid possible biases coming from outside the target volume. The ACH method thus does not really retrieve the 3-D structure of the target volume, but instead leads to velocity contrasts relative to the layer average of the velocity, this average value remaining unknown ( Aki et al. 1977 ). Two artefacts derive from this particularity: (1) velocity contrasts are known only in the horizontal direction and it is not possible, in a strict mathematical sense, to estimate the contrasts in the vertical direction with ACH alone; (2) negative anomalies are often interpreted as low velocities, whereas negative anomalies may correspond to high velocities if the average value of the corresponding layer is sufficiently high. The converse is true of positive anomalies. We show with synthetic data how these artefacts can affect the interpretation of tomographic images. We propose to correct the artefacts by reintroducing the 1-D regional average model, and show in synthetic experiments how effective this correction can be.
  The application of this procedure to data recorded in the Kunlun region shows that the retrieval of the absolute values of the 3-D velocity model is helpful for interpreting the tomographic images and better defining which features are anomalous.     

Small-scale heterogeneities below the Gräfenberg array, Germany from seismic wavefield fluctuations of Hindu Kush events

Small-scale elastic heterogeneities (<5  km) are found in the upper lithosphere underneath the Gräfenberg array, southeast Germany. The results are based on the analysis of broadband recordings of 17 intermediate-depth (201–272  km) events from the Hindu Kush region. The wavefront of the first P arrival and the following 40  s coda are separated into coherent and incoherent (scattered) parts in the frequency range from 0.05 to 5  Hz. The frequency-dependent intensities of the mean and fluctuation wavefields are used to describe the scattering characteristics of the lithosphere underneath the receivers. It is possible to discriminate a weak-fluctuation regime of the wavefield in the frequency range below approximately 1.5–2.5  Hz and a strong-fluctuation regime starting at 2.0–2.5  Hz and continuing to higher frequencies. In order to explain the observed wavefield fluctuations, an approach with seismic scattering at random media-type structures is proposed. The preferred model contains heterogeneities with 3–7 per cent perturbations in seismic velocity and correlation lengths of 0.6–4.8  km in the crust. This is compatible with models from active seismic experiments. Scattering in the lithospheric mantle is not required, but cannot be excluded at weak velocity contrasts (<3 per cent).     

Scattering of surface waves modelled by the integral equation method

The integral equation method is used to model the propagation of surface waves in 3-D structures. The wavefield is represented by the Fredholm integral equation, and the scattered surface waves are calculated by solving the integral equation numerically. The integration of the Green's function elements is given analytically by treating the singularity of the Hankel function at   R = 0  , based on the proper expression of the Green's function and the addition theorem of the Hankel function. No far-field and Born approximation is made. We investigate the scattering of surface waves propagating in layered reference models imbedding a heterogeneity with different density, as well as Lamé constant contrasts, both in frequency and time domains, for incident plane waves and point sources.     

Time-lapse traveltime change of singly scattered acoustic waves

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

Computational fluid dynamics and the physical modelling of an upland urban river

This paper describes the application of a commercially available, three-dimensional computational fluid dynamic (CFD) model to simulate the flow structure in an upland river that is prone to flooding. Simulations use a rectangular channel geometry, smooth sidewalls and a bed topography obtained from the field site that contains a subdued pool–riffle sequence. The CFD model uses the RNG κ– turbulence closure scheme of Yakhot and Orszag (J. Sci. Comput. 1 (1986) 1), as implemented in FLUENT 4.4.4, with a free surface. Results are shown for numerical runs simulating a 1:100 year return interval flood. Output from the numerical model is compared to a physical model experiment that uses a 1:35 scale fibreglass mould of the field study reach and measures velocity using ultrasonic Doppler velocity profiling (UDVP). Results are presented from the numerical and flume models for the water surface and streamwise velocity pattern and for the secondary flows simulated in the numerical model. A good agreement is achieved between the CFD model output and the physical model results for the downstream velocities.Results suggest that the streamwise velocity is the main influence on the flow structure at the discharge and channel configuration studied. Secondary flows are, in general, very weak being below the resolution of measurement in the physical model and less than 10% of the streamwise velocity in the numerical model. Consequently, there is no evidence for a ‘velocity dip’. It is suggested that the subdued topography or inlet morphology may inhibit the development of secondary flows that have been recorded in previous flat-bed, rectangular open channel flows. A significant corollary of these results is that the morphological evolution of the pool–riffle sequence at high discharges may be controlled primarily by the downstream distribution of velocity and sediment transport with little role for lateral sorting and sediment routing by secondary flows. This paper also raises a number of issues that may be of use in future CFD modelling of three-dimensional flow in open channels within the geomorphological community.     

Wavepath traveltime tomography

The elastic-wave equation is used to construct sensitivity kernels relating perturbations in elastic parameters to traveltime deviations. Computation of the functions requires a correlation of the forward-propagating seismic wavefield with a backward propagation of the residual wavefield. The computation of the wavefields is accomplished using a finite difference algorithm and is efficiently executed on a CM-2 parallel processor. The source and receiver locations have maximum sensitivity to velocity structure. The sensitivity kernels or wavepaths are well suited for transmission traveltime inversion such as cross-borehole tomography and vertical seismic profiling. Conventional ray tomography and wavepath tomography are applied to a set of P -wave arrival times, from a cross-borehole experiment at Kesterson, California. Because the wavepaths have increased sensitivity near the source and receiver there are differences in resolution of the velocity structure. Both techniques recover the same relative variations in velocity where the coverage is adequate. The wavepath solution is more laterally continuous and the dominant variation is vertical, as is expected for the layered sediments in this region.