Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space

Large scale seismic anisotropy in the Earth's mantle is likely dynamically supported by the mantle's deformation; therefore, tomographic imaging of 3-D anisotropic mantle seismic velocity structure is an important tool to understand the dynamics of the mantle. While many previous studies have focused on special cases of symmetry of the elastic properties, it would be desirable for evaluation of dynamic models to allow more general axis orientation. In this study, we derive 3-D finite-frequency surface wave sensitivity kernels based on the Born approximation using a general expression for a hexagonal medium with an arbitrarily oriented symmetry axis. This results in kernels for two isotropic elastic coefficients, three coefficients that define the strength of anisotropy, and two angles that define the symmetry axis. The particular parametrization is chosen to allow for a physically meaningful method for reducing the number of parameters considered in an inversion, while allowing for straightforward integration with existing approaches for modelling body wave splitting intensity measurements. Example kernels calculated with this method reveal physical interpretations of how surface waveforms are affected by 3-D velocity perturbations, while also demonstrating the non-linearity of the problem as a function of symmetry axis orientation. The expressions are numerically validated using the spectral element method. While challenges remain in determining the best inversion scheme to appropriately handle the non-linearity, the approach derived here has great promise in allowing large scale models with resolution of both the strength and orientation of anisotropy.     

Cross-dependence of finite-frequency compressional waveforms to shear seismic wave speeds

Seismic tomography has been one of the primary tools to image the interior of the earth and other elastic structures. To date the inversions of compressional ( P ) and shear ( S ) wave speeds have been carried out separately under the assumption that P traveltimes are affected only by the P wave speed of the elastic media and S traveltimes by the S wave speed. Using numerical and analytical solutions, we show that for finite-frequency seismic waves, S wave speed perturbations may have significant effects on P waveforms. This suggests that when waveform-derived traveltime and amplitude anomalies are used in tomographic inversions, the P -wave measurements should be related to not only P wave speed perturbations but also S wave speed perturbations.     

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.     

Crustal scattering at the KTB from a combined microearthquake and receiver analysis

In December of 1994 a fluid injection experiment which triggered several hundreds of microearthquakes was conducted at the KTB main borehole (Oberpfalz, Germany). These events were recorded with a temporal seismic network at the surface. Out of the complete data set, a cluster of five events recorded at four mini-arrays consisting of eight or nine stations was used to investigate the crustal scattering properties in the vicinity of the KTB. For this purpose, the 'Double Beam Method' (DBM; Krüger et al . 1993 , 1995 , 1996) and the 'Double Beam Imaging Method' (Scherbaum, Krüger & Weber 1997) were extended to curved wave fronts to drop the restriction of plane-wave propagation. This technique is used for imaging the crustal scattering strength using earthquake clusters recorded at close-by mini-arrays. The results of the array analysis show that the composition of the P coda is mainly affected by the site location of the arrays. Near-surface and deeper crustal scattering contribute in a very complicated pattern. Furthermore, with the present data set it was possible to identify reflections from the top of the Erbendorf Body. This is a very pronounced arrival in most of the recorded traces. In one of the arrays its amplitudes are even greater than the direct P phases. Five to eight coherent phases could be identified by the mini-arrays. Using only these phases, synthetic P -coda traces were constructed, which only contain the coherent part of the observed wavefield. By subtracting the synthetic coherent wavefield from the original traces we achieve a variance reduction in the P coda of up to 37 per cent. This leads to the conclusion that a large amount of the P coda at the KTB can be modelled by a simple deterministic single-scattering model using a small number of individual scatterers.     

Seismic waves in rocks with fluids and fractures

Seismic wave propagation through the earth is often strongly affected by the presence of fractures. When these fractures are filled with fluids (oil, gas, water, CO₂, etc.), the type and state of the fluid (liquid or gas) can make a large difference in the response of the seismic waves. This paper summarizes recent work on methods of deconstructing the effects of fractures, and any fluids within these fractures, on seismic wave propagation as observed in reflection seismic data. One method explored here is Thomsen's weak anisotropy approximation for wave moveout (since fractures often induce elastic anisotropy due to non-uniform crack-orientation statistics). Another method makes use of some very convenient crack/fracture parameters introduced previously that permit a relatively simple deconstruction of the elastic and wave propagation behaviour in terms of a small number of crack-influence parameters (whenever this is appropriate, as is certainly the case for small crack densities). Then, the quantitative effects of fluids on these crack-influence parameters are shown to be directly related to Skempton's coefficient B of undrained poroelasticity (where B typically ranges from 0 to 1). In particular, the rigorous result obtained for the low crack density limit is that the crack-influence parameters are multiplied by a factor  (1 − B )  for undrained systems. It is also shown how fracture anisotropy affects Rayleigh wave speed, and how measured Rayleigh wave speeds can be used to infer shear wave speed of the fractured medium in some cases. Higher crack density results are also presented by incorporating recent simulation data on such cracked systems.     

A Preliminary Investigation into the Relationship between Long-Period Seismic Noise and Local Fluctuations in the Atmospheric Pressure Field

Summary The displacement response of an elastic half space to a plane pressure wave is examined in order to establish the conditions under which sources of this type can contribute significantly to the long-period seismic noise field. The study is restricted to pressure waves which propagate at velocities well below the seismic wave velocities characteristic of the half space. The numerical studies indicate that pressure waves with amplitudes of 100 μbar or more can contribute significantly to the long-period vertical background noise observed at the surface, provided that the detectors are located on sections of alluvial fill or poorly to moderately indurated sandstones and shales whose thicknesses are greater than about a kilometre. These same waves can also create significant tilt noise on long-period horizontal seismographs located at or near the surface, regardless of the rock type. The seismic disturbances created by pressure waves decay rapidly away from the surface. Therefore, it appears that it may be possible to eliminate the effects of atmospherically generated noise by placing the detectors at moderate depths.     

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

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

Acoustoelasticity in cracked solids

A new model that accounts for the stress dependence of the phase velocity of elastodynamic waves propagating in a cracked solid under compression is presented. The phase velocities of longitudinal and shear waves are derived from the effective elastic properties of a cracked solid, which are evaluated within the framework of Kachanov's approach. Following Kachanov, the extra-compliance tensor of the cracked solid is related to the crack compliances, which display a marked non-linear behaviour when subjected to a compressive load. Such non-linear behaviour is shown to be derived from the elastic interaction between the contacting crack faces under compression. This work does not address the effect of mutual interaction among cracks and the generation of higher harmonics due to the medium non-linearity. Numerical examples are presented that illustrate the phase velocity changes occurring in a solid with a random distribution of parallel cracks as a function of an external compressive load. A distinctive feature of the acoustoelastic effect in solids with large parallel fractures and in solids with distributions of aligned microcracks is also illustrated.     

Flexure of the ocean lithosphere from island uplift, bathymetry and geoid height observations: the Society Islands

Summary. The response of a stratified elastic medium can be conveniently characterized by the Green's tensor for the medium. For coupled seismic wave propagation ( P—SV or fully anisotropic), the Green's tensor may be constructed directly from two matrices of linearly independent displacement solutions. Rather simple forms for the Green's tensor can be found if each displacement matrix satisfies one of the boundary conditions on the seismic field. This approach relates directly to 'reflection matrix' representations of the seismic field.
For a stratified elastic half space the Green's tensor is used to give a spectral representation for coupled seismic waves. By means of a contour integration a general completeness relation is obtained for the 'body wave' and 'surface wave' parts of the seismic field. This relation is appropriate for SH and P–SV waves in an isotropic medium and also for full anisotropy.     

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

Seismic evidence for a sharp lithospheric base persisting to the lowermost mantle beneath the Caribbean

Broad-band data from South American earthquakes recorded by Californian seismic networks are analysed using a newly developed seismic wave migration method—the slowness backazimuth weighted migration (SBWM). Using the SBWM, out-of-plane seismic P -wave reflections have been observed. The reflection locations extend throughout the Earth's lower mantle, down to the core–mantle boundary (CMB) and coincide with the edges of tomographically mapped high seismic velocities. Modelling using synthetic seismograms suggests that a narrow (10–15 km) low- or high-velocity lamella with about 2 per cent velocity contrast can reproduce the observed reflected waveforms, but other explanations may exist. Considering the reflection locations and synthetic modelling, the observed out-of-plane energy is well explained by underside reflections off a sharp reflector at the base of the subducted lithosphere. We also detect weaker reflections corresponding to the tomographically mapped top of the slab, which may arise from the boundary between the Nazca plate and the overlying former basaltic oceanic crust. The joint interpretation of the waveform modelling and geodynamic considerations indicate mass flux of the former oceanic lithosphere and basaltic crust across the 660 km discontinuity, linking processes and structure at the top and bottom of the Earth's mantle, supporting the idea of whole mantle convection.     

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.     

Using evidence of non-linear induced polarization for detecting extended ore mineralizations

The discrimination between electrolytic and electronic conductors is highly relevant to geological modelling as it allows conclusions to be drawn about the formation and mineral composition of rocks. The induced polarization (IP) method, which compares the electric current injected into the ground with the corresponding earth potential differences can be used for this purpose.
  This paper describes a new method based on the theory that non-linear electrochemical processes on the surface of electronic conductors are responsible for non-linear IP (NLIP) phenomena. This results in multiples of the fundamental frequency being observed in the telluric voltage spectra when a monochromatic current signal is fed into the ground. The non-linearity of the current–voltage characteristic is most effectively described by a spectral method.
  A laboratory experiment was carried out, using an electrolytic trough with a small graphite cylinder serving as an electronic conductor, which clearly demonstrated the validity of the method. A field experiment was undertaken at a borehole of approximately 450  m depth, located in the transition zone of the Tepla-Barrandium and Moldanubicum in East Bavaria. A sinusoidal current was injected into the ground using a logging tool at depths varying between 150 and 450  m. The corresponding potential differences were simultaneously observed along a profile on the surface. Field and laboratory results show a striking similarity. It can be concluded that an extensive electronic conductor—probably graphite—is steeply dipping southwards meeting the borehole at approximately 310  m depth.     

Scattering of waves in elastic media containing passive and active cracks

Scattering of wavefields in a 3-D medium that includes passive and/or active structures, is numerically solved by using the boundary integral equation method (BIEM). The passive structures are velocity anomalies that generate scattered waves upon incidence, and the active structures contain endogenous fracture sources, which are dynamically triggered by the dynamic load due to the incident waves. Simple models are adopted to represent these structures: passive cracks act as scatterers and active cracks as fracture sources. We form cracks using circular boundaries, which consist of many boundary elements. Scattering of elastic waves by the boundaries of passive cracks is treated as an exterior problem in BIEM. In the case of active cracks, both the exterior and interior problems need to be solved, because elastic waves are generated by fracturing with stress drop, and the growing crack boundaries scatter the incident waves from the outside of the cracks. The passive cracks and/or active cracks are randomly distributed in an infinite homogeneous elastic medium. Calculations of the complete waveform considering a single scatter show that the active crack has weak influence on the attenuation of first arrivals but strong influence on the amplitudes of coda waves, as compared with those due to the passive crack. In the active structures, multiple scattering between cracks and the waves triggered by fracturing strongly affect the amplitudes of first arrivals and coda waves. Compared to the case of the passive structures, the attenuation of initial phase is weak and the coda amplitudes decrease slowly.     

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.     

S to P scattering at the 650 km discontinuity

Summary. A search of seismograms recorded at the Warramunga seismic array (WRA) from events occurring below the Izu-Bonin Islands shows an arrival on some, but not all, of the records, with an onset at 25–30 s after P , which is not predicted by the standard travel-time tables. The slowness and azimuth of the phase show that it is generated almost in line with P , and the variation of arrival time with the hypocentral depth of the earthquake indicates that its origin lies on the receiver side of the source. It appears, in fact, to be an S to P conversion at a depth of 650–700 km, which is seen only when the receiver is close to a node of the P radiation pattern and an antinode for S so that its amplitude compared with that of P is at a maximum.
Finally, the duration of the phase indicates that it is not simply a refracted wave, but that it has a coda of scattered arrivals from lateral heterogeneity in the neighbourhood of 650 km below the Izu-Bonin Islands.     

Finite difference seismograms for laterally varying marine models

Summary. The method of finite differences is applied to the elastic wave equation to generate synthetic seismograms for laterally varying seafloor structures. The results are compared with borehole seismic data from the Gulf of California (Deep Sea Drilling Project Site 485) in which lines were shot over flat and rough topography. The significant new phenomenon observed in both the synthetic seismograms and the field data is the generation of a 'double head wave' due to the interaction of the incident wavefront with the side of a hill and the flat seafoor adjacent to the hill.
In these models the hills are on the order of a seismic wavelength in height and steep velocity gradients occur over distances comparable to wavelengths. Ray theoretical methods would not be suitable for studying such structures. True amplitude record sections are obtained by the finite difference method, which show for these models that the head wave generated at the flat seafloor adjacent to the hill is lower in amplitude than if the hill were not present and is lower in amplitude than the head wave generated at the hill.
A second feature which is important for borehole receivers is the existence of the 'direct wave root' in the upper basement. This energy occurs below the sharp interface when the direct wave impinges on the interface from above. There is no corresponding Snell's law ray path for this energy and the energy is evanescent with depth in the lower medium.
The properties of both the double head wave and the direct wave root are clearly demonstrated in the finite difference 'snapshot' displays.     

A data-mapping approach to 3-D wavefield redatuming

In case of a complex overburden, the seismic data can be greatly improved by applying a full wavefield redatuming procedure. In practice, the application of the redatuming process to 3-D data acquired by conventional acquisition designs is non-trivial. Because of the large amount of data involved in the 3-D redatuming process and because of the sparseness of these data, it is impossible to apply conventional wave equation datuming directly.
We present a data mapping approach to redatuming (DMR), which follows the concept of Kirchhoff data mapping. A simplified background medium where no ray bending occurs is assumed for the medium below the datum in order to map an input data set referenced to the acquisition surface to an output data set referenced to the new datum level. The DMR method can be interpreted as a simplified version of the Kirchhoff summation redatuming (KSR) method, where one of the 2-D integrals over the acquisition coordinates can be solved analytically. Consequently, in this approach fewer traces are involved in the computation of one time sample (a 2-D integral is computed instead of a 4-D integral), which makes it particularly attractive for the application to 3-D data sets.
In this paper the theory underlying data mapping redatuming is discussed and the proposed approach is tested on fully sampled 2-D and 3-D synthetic data from models with both simple and complex velocity distributions in the subsurface.
The tests clearly show that the objective of producing results that are comparable to the conventional KSR has been achieved. The redatumed traces are dynamically and kinematically correct. Furthermore, these results confirm that the dependency of the new approach on the assumed medium below the datum level is, indeed, weak because the assumption of a velocity medium where no ray bending occurs is already sufficient to produce correct results.     

Shear wave anisotropy in the crust around the San Andreas fault near Parkfield: spatial and temporal analysis

We systematically analysed shear wave splitting (SWS) for seismic data observed at a temporary array and two permanent networks around the San Andreas Fault (SAF) Observatory at Depth. The purpose was to investigate the spatial distribution of crustal shear wave anisotropy around the SAF in this segment and its temporal behaviour in relation to the occurrence of the 2004 Parkfield M 6.0 earthquake. The dense coverage of the networks, the accurate locations of earthquakes and the high-resolution velocity model provide a unique opportunity to investigate anisotropy in detail around the SAF zone. The results show that the primary fast polarization directions (PDs) in the region including the SAF zone and the northeast side of the fault are NW–SE, nearly parallel or subparallel to the SAF strike. Some measurements on the southwest side of the fault are oriented to the NNE–SSW direction, approximately parallel to the direction of local maximum horizontal compressive stress. There are also a few areas in which the observed fast PDs do not fit into this general pattern. The strong spatial variations in both the measured fast PDs and time delays reveal the extreme complexity of shear wave anisotropy in the area. The top 2–3 km of the crust appears to contribute the most to the observed time delays; however substantial anisotropy could extend to as deep as 7–8 km in the region. The average time delay in the region is about 0.06 s. We also analysed temporal patterns of SWS parameters in a nearly 4-yr period around the 2004 Parkfield main shock based on similar events. The results show that there are no appreciable precursory, coseismic, or post-seismic temporal changes of SWS in a region near the rupture of an M 6.0 earthquake, about 15 km away from its epicentre.