Seismic Ray Research Articles

Inversion of velocity models using genetic algorithm method with sigmoidal parameterization

A seismic traveltime inversion method is proposed for building smooth velocity models using traveltime observed on irregular surface. Model parameterization in this study is described by a piecewise constant velocity field on a rectangular grid parameterized by sigmodal functions, which is beneficial for the description of irregular surface with high degree of approximation. The velocity field is defined in the rectangular grid which is used for the description of velocity distribution everywhere in the model sigmoidal interpolation. In addition, we use the simple Genetic Algorithm for the inversion procedure. Through this global scope inversion method, we provide high-resolution estimates of the model parameter and ensure that the results obtained are in accordance with the actual data. Our method is validated with synthetic examples of heterogeneous isotropic media and compared to Simulating Annealing. The inverted velocity models and approximate ray paths obtained coincide well with the trajectories simulated using the seismic ray tracing in synthetic heterogeneous isotropic media.

Evaluation of 3D seismic survey design parameters through ray-trace modeling and seismic illumination studies: a case study

This paper describes a case study of wavefront construction-based ray-trace modeling to access the 3D seismic exploration parameters that significantly impact achieving the exploration target in seismic data acquisition. The traditional methods assessment is based on the horizontal reflector concept, which does not consider the subsurface's inhomogeneities. This case study provides a methodology by considering the effect of subsurface variations on the estimation of seismic survey parameters. As the first step in this methodology, an elastic earth model is created to propagate theoretical seismic rays from a 3D seismic survey to generate seismic ray properties. Then illumination maps and synthetic seismic sections are generated from these seismic ray attributes to evaluate seismic survey performance in seismic imaging targets. The results proved that the method helps estimate seismic survey efficiencies in seismic target imaging. Therefore, the ray-trace modeling methodology can help obtain the target fold coverage in complex geological settings by designing and verifying the seismic survey parameters.

Slab control on the mega-sized North Pacific ultra-low velocity zone

Ultra-low velocity zones (ULVZs) are localized small-scale patches with extreme physical properties at the core-mantle boundary that often gather at the margins of Large Low Velocity Provinces (LLVPs). Recent studies have discovered several mega-sized ULVZs with a lateral dimension of ~900 km. However, the detailed structures and physical properties of these ULVZs and their relationship to LLVP edges are not well constrained and their formation mechanisms are poorly understood. Here, we break the degeneracy between the size and velocity perturbation of a ULVZ using two orthogonal seismic ray paths, and thereby discover a mega-sized ULVZ at the northern edge of the Pacific LLVP. The ULVZ is almost double the size of a previously imaged ULVZ in this region, but with half of the shear velocity reduction. This mega-sized ULVZ has accumulated due to stable mantle flow converging at the LLVP edge driven by slab-debris in the lower mantle. Such flow also develops the subvertical north-tilting edge of the Pacific LLVP.

A strange symbiosis: Seismic ray theory and Saturday touch football — Donald V. Helmberger as a graduate advisor

A strange symbiosis: Seismic ray theory and Saturday touch football — Donald V. Helmberger as a graduate advisor

Ray velocity derivatives in anisotropic elastic media – Part I: general anisotropy

SUMMARY We present an original, generic and efficient approach for computing the first and second partial derivatives of ray (group) velocities along seismic ray paths in general anisotropic (triclinic) elastic media. As the ray velocities deliver the ray element traveltimes, this set of partial derivatives constructs the so-called kinematic and dynamic sensitivity kernels which are used in different key seismic modelling and inversion methods, such as two-point ray bending methods and seismic tomography. The second derivatives are useful in the solution of the above-mentioned kinematic problems, and they are essential for evaluating the dynamic properties along the rays (amplitudes and phases). The traveltime is delivered through an integral over a given Lagrangian defined at each point along the ray. In our approach, we use an arclength-related Lagrangian representing a reciprocal of the ray velocity magnitude. Although this magnitude cannot be explicitly expressed in terms of the medium properties and the ray direction components, its derivatives can still be formulated analytically using the corresponding arclength-related Hamiltonian that can be explicitly expressed in terms of the medium properties and the slowness vector components; this requires first to obtain (invert for) the slowness vector components, given the ray direction components. Computation of the slowness vector and the ray velocity derivatives is considerably simplified by using an auxiliary scaled-time-related Hamiltonian obtained directly from the Christoffel equation and connected to the arclength-related Hamiltonian by a simple scale factor. This study consists of two parts. In Part I, we consider general anisotropic (triclinic) media, and provide the derivatives (gradients and Hessians) of the ray velocity, with respect to (1) the spatial location and direction vectors and (2) the elastic model parameters. The derivatives are obtained for both quasi-compressional and quasi-shear waves, where other types of media, characterized with higher symmetries, can be considered particular cases. In Part II, we apply the theory of Part I explicitly to polar anisotropic media (transverse isotropy with tilted axis of symmetry, TTI), and obtain the explicit ray velocity derivatives for the coupled qP and qSV waves and for SH waves. The derivatives for polar anisotropy are simplified (as compared to general anisotropy), obviously yielding more effective computations. The ray velocity derivatives are tested by checking consistency between the proposed analytical formulae and the corresponding numerical ones.

Approximations to a generalized real ray-tracing method for heterogeneous anisotropic viscoelastic media

The real ray-tracing approach leads to an effective solution in the real space domain using a homogeneous ray velocity vector. However, it fails to yield solutions for quasi-S-waves, which suffer from triplications of the wavefronts. To address this challenging problem, a new real ray-tracing method and its two approximations are developed to solve the complex ray equation. The numerical results indicate that the new real ray-tracing method is superior to the common real ray-tracing method in the presence of triplications of the quasi-S-waves in the computation of ray velocity, ray attenuation, and ray quality factors, as well as the reflection and transmission coefficients in viscoelastic anisotropic media. Based on the assumptions of the real slowness direction and real polarization vectors, two new approximations of the new real ray-tracing method are developed for directly computing the homogeneous complex ray velocity vectors of three wave modes (qP, qS1, and qS2). These approximations significantly improve computational efficiency by avoiding the iterative process required by the new real ray-tracing method that is inherited from the common real ray-tracing method. The computational accuracies are verified through transversely isotropic and orthorhombic models with different strengths of attenuation and anisotropy. Incorporation of the new approximation into the shortest-path method turned out to be efficient and accurate for seismic ray tracing in heterogeneous viscoelastic and transversely isotropic media with a vertical axis of symmetry, even in the presence of strong attenuation and anisotropy.

Application of the Stacked Refraction Convolution Section to 2D Ocean Bottom Seismometer Wide-angle Seismic Data Along the Tamayo Through Basin, Gulf of California

The stacked refraction convolution section can be used as an interpretation tool in wide-angle refraction seismic data generated by air gun shooting and recorded by Ocean Bottom Seismometers (OBS). The refraction convolution section is a full-wave extension of the Generalized Reciprocal Method (GRM), a method frequently used in shallow refraction seismic interpretation, but not applied to deep crustal-scale studies. The sum of the travel times of the waves refracted in the same interface and recorded in a pair of forward and reverse profiles, time-corrected by the reciprocal time, is an estimation close to the two-way travel times of the multichannel seismic reflection sections, but with seismic rays illuminating the interfaces upwards. The sum of seismic traces is obtained with the convolution section. Furthermore, several pairs of convolved forward-reverse refraction recordings of the same area can be stacked together to improve the signal to noise ratio. To show the applicability of the refraction convolution section in OBS deep data, we interpreted the basement structure of the Tamayo Through Basin in the southern Gulf of California, offshore Mexico. We compared the results with both, a multichannel seismic section recorded in the same profile, and the previous interpretations of the same wide-angle seismic data modeled with ray tracing and tomography methods. The basement imaged by the stacked refraction convolution section is similar in geometry to that obtained by seismic reflection processing. The stacked refraction convolution section identifies the full extent of the basement and confirms the location of a nearly constant thickness volcanic layer in the northwestern half of the basin. However, only a small area of volcanic deposits is found in the shallower parts of the southwestern margin. We also show that the convolution process can be used to estimate the occurrence of lateral variations of seismic velocities in the basement, as a further application of the GRM to deep refraction data.

A systematic investigation of piercing-point-dependent seismic azimuthal anisotropy

SUMMARYThe vast majority of teleseismic XKS (including SKS, SKKS and PKS) shear wave splitting studies interpret the observed splitting parameters (fast orientation and splitting time) based on the assumption of a spatially invariant anisotropy structure in the vicinity of a recording station. For such anisotropy structures the observed splitting parameters are either independent of the arriving azimuth of the seismic ray paths if the medium traversed by the ray paths can be represented by a single layer of anisotropy with a horizontal axis of symmetry (i.e. simple anisotropy), or demonstrate a periodic variation with respect to the arriving azimuth for a more complicated structure of anisotropy (e.g. multiple layers with a horizontal axis of symmetry, or a single layer with a dipping axis). When a recording station is located near the boundary of two or more regions with different anisotropy characteristics, the observed splitting parameters are dependent on the location of the ray piercing points. Such a piercing-point dependence is clearly observed using a total of 360 pairs of XKS splitting parameters at three stations situated near the northeastern edge of the Sichuan Basin in central China. For a given station, the fast orientations differ as much as 90°, and the azimuthal variation of the fast orientations lacks a 90° or 180° periodicity which is expected for double-layered or dipping axis anisotropy. The observed splitting parameters from the three stations are spatially most consistent when they are projected at a depth of ∼250 km, and can be explained by shear strain associated with the absolute plate motion and mantle flow deflected by the cone-shaped lithospheric root of the Sichuan Basin.

An improved LTI algorithm for multiphase raytracing in undulating surface model

Abstract Based on the linear traveltime interpolation (LTI) algorithm, we propose an improved LTI (referred as ILTI) algorithm for multiphase seismic ray tracing, which uses a velocity node in model parameterization and introduces secondary nodes between adjacent velocity nodes. To better fit the undulating surface model, an irregular velocity node is used in near the irregular interface, and regular velocity nodes are still used in the region far away from the irregular interface. We derive an iterative fixed-point formula for calculating traveltime. By combining multistage computational technology and wavefront narrowband expansion technology, the proposed ILTI algorithm can efficiently trace the multiphase seismic raypath and compute the corresponding traveltime field. Through comparison and analysis with the traditional LTI algorithm, its computational accuracy can be highlighted by at least one order of magnitude. Compared with the popular fast marching method (FMM) and irregular shortest-path (ISPM) algorithms, it also has the advantages in terms of computational accuracy and efficiency. Numerical simulations in the Marmousi model show that the algorithm is also suitable for tracking multiphase seismic rays in the complex velocity model.

Dense seismic arrays deny a massive magma chamber beneath the Taipei metropolis, Taiwan

Several recent studies suggest that the Tatun Volcano Group (TVG) in the Taipei metropolis of Taiwan is still active with a mappable magma chamber beneath it. Here we report new seismic evidence from dense seismic arrays in northern Taiwan to refute the presence of a massive magma chamber. We investigated two near Taipei earthquakes with focal depths of ~ 140 km. We found that all the waveforms exhibited distinct S waves even when they traversed across the previously postulated magma chamber. Instead, the S-wave shadows found in the previous study may result from seismic waves traveling through a magma diapir above the subducting Philippine Sea Plate offshore northern Taiwan. Moreover, we found the P-wave delay increased with hypocentral distance when the seismic waves propagated through the footwall (west side) of the Shanchiao fault, regardless of whether they traversed across the postulated magma chamber. Our study results also indicate no abnormal attenuation when seismic rays traversed across the postulated magma chamber. Furthermore, the average {Q}_{P}/{Q}_{S} ratio around the TVG is less than 1, which implies that scattering attenuation is dominant. We conclude that a highly fractured rock body is beneath the TVG with a tiny fraction of magma instead of a massive magma chamber. Without sufficient magma supply, the TVG may stay dormant (except for small phreatic eruptions).

Efficient seismic ray tracing based on the shortest path method

SUMMARYLimiting the number of valid nodes around a reference ray can greatly reduce the calculation time of the shortest-path method (SPM). The calculation is executed by increasing the number of cells and/or nodes in the target area, step by step, until the rays converge. The ray obtained in the previous step is used as a reference ray, and the initial reference rays are given by the ordinary SPM. The Dijkstra algorithm and binary heap sorting method are used, as in the ordinary SPM. As the cell and node numbers increase, the calculation time for the modified SPM (mSPM) is reduced compared to that for the ordinary SPM. In the 3-D (100 × 100 × 100 km3) checkerboard velocity pattern model, the relative calculation time becomes two to four orders of magnitude smaller. The calculation time for the mSPM itself is approximately proportional to E log2(V), where E is the edge (ray path segment) number and V is the vertex (node) number, as seen in the heap sorting algorithm. The mean traveltime and ray path differences between the mSPM and pseudo-bending method (PBM) are small, less than 0.005 s and around 0.6 km, respectively, and slightly larger than those with respect to the ordinary SPM. The total differences from the exact solution are estimated to be less than 0.01 s and 1.0 km, which are sufficiently small for traveltime tomography. The traveltime and ray path can be improved by utilizing iterative calculations, shifting of the starting point and more neighbour nodes. The ray path obtained by the mSPM can be a local minimum, according to velocity models. The relative mean traveltime and ray path differences between adjacent cells and nodes generally show trends similar to those of the traveltime and ray path differences from the PBM. Hence, these relative differences can help reveal the behaviour of the differences from the PBM.

P-wave velocities of the upper mantle of the Tethysalpine geosynclines

The aim of the work is to calculate the velocity distribution of the longitudinal seismic waves (VP) in the upper mantle at depths from the discontinuity M to 400 km. The object is the territory of Tethys — the belt of alpine geosynclines crossing all of Eurasia from Gibraltar to the Indonesian archipelago.The model of the first approximation is constructed according to previous studies and our results on island arcs. It was possible to choose the velocity distribution in the upper mantle of Tethys, according to which the travel time close to the observed one was calculated. The degree of coordination is quite sufficient for the recognition of the used distribution of velocity as real. The data on approximately 18,000 earthquakes, the waves from which reached 27 seismic stations, are used. The location of the epicenters of earthquakes and stations ensured the passage of seismic rays precisely along the upper mantle of different Tethys regions. On the island of Sumatra a comparison was made of the constructed travel times with those obtained earlier in the study of the upper mantle of island arcs and coastal ridges of the Pacific Ring. Harmonization of the travel times is considered satisfactory. The VP distribution obtained in this way is compared with the calculated from the ideas of the advective-polymorphic hypothesis. The velocity model of the upper mantle of the Precambrian platform outside the zones of modern activation is complemented by the influence of the anomalous temperatures of both signs in the subsoil of the geosyncline. The result is somewhat different from the calculated earlier for the upper mantle island arc. There are two main reasons for this difference: probability of a somewhat higher radiogenic heat generation in the rocks of the upper mantle of the arc and somewhat smaller typical crustal thickness under the arc. The discrepancies between the models are small, averaging about 0.03 km/s. They do not exceed associated only with the errors of calculations. A noticeable excess of the magnitude of the discrepancies is detected only at a depth of 400 km, where the uncertainty of the calculation results is maximum.The discrepancy between the experimental and calculated velocity distributions in the upper mantle of Tethys is 0.05 km/s. They can be considered coincident.

Effective and efficient approaches for calculating seismic ray velocity and attenuation in viscoelastic anisotropic media

In viscoelastic anisotropic media, the elastic moduli, slowness vector, phase, and ray velocity are all complex-valued quantities in the frequency domain. Solving the complex eikonal equation becomes computationally complex and time-consuming. We have developed two approximate methods to effectively calculate the ray velocity vector, attenuation, and quality factor in viscoelastic transversely isotropic media with a vertical symmetry axis (VTI) and in orthorhombic (ORT) anisotropy. The first method is based on the perturbation theory (PER) under the assumption of a homogeneous complex ray vector, which is obtained by applying the elastic background and viscoelastic perturbations to the real and imaginary components of the modulus tensor, respectively. The perturbations of the slowness vectors of the three wave modes (qP, qSV, and qSH) are determined through the vanishing Hamiltonian function. The second method is derived by applying a real slowness direction (RSD) to the inhomogeneous complex slowness vector and then approximately calculating the complex ray velocity vector with the condition of the homogeneous complex vector. The numerical results verify that the two approaches can produce accurate ray velocity vector, attenuation, and quality factors of the qP-wave in viscoelastic VTI and ORT media. The RSD method can yield high accuracies of ray velocity for the qSV- and qSH-wave in viscoelastic VTI models even at triplication of the qSV wavefronts, as well as qS1 and qS2 in a weak ORT medium ([Formula: see text] > 20), except for near the cusp of the qS1 wavefronts (errors approximately 6%) where the PER has more than 10% error.

Seismic complex ray tracing in 2D/3D viscoelastic anisotropic media by a modified shortest-path method

In a viscoelastic anisotropic medium, velocity anisotropy and wave energy attenuation occur and are often observed in seismic data applications. Numerical investigation of seismic wave propagation in complex viscoelastic anisotropic media is very helpful in understanding seismic data and reconstructing subsurface structures. Seismic ray tracing is an effective means to study the propagation characteristics of high-frequency seismic waves. Unfortunately, most seismic ray-tracing methods and traveltime tomographic inversion algorithms only deal with elastic media and ignore the effect of viscoelasticity on the seismic raypath. We have developed a method to find the complex ray velocity that gives the seismic ray speed and attenuation in an arbitrary viscoelastic anisotropic medium, and we incorporate them with the modified shortest-path method to determine the raypath and calculate the real and imaginary traveltime (wave energy attenuation) simultaneously. We determine that the complex ray-tracing method is applicable to arbitrary 2D/3D viscoelastic anisotropic media in a complex geologic model and the computational errors of the real and imaginary traveltime are less than 0.36% and 0.59%, respectively. The numerical examples verify that the new method is an effective and powerful tool for accomplishing seismic complex ray tracing in heterogeneous viscoelastic anisotropic media.

Seismic velocity precursors to the 2016 Mw 6.5 Norcia (Italy) earthquake

AbstractEarthquakes occur as the result of long-term strain accumulation on active faults and complex transient triggering mechanisms. Although laboratory experiments show accelerating deformation patterns before failure conditions are met, imaging similar preparatory phases in nature remains difficult because it requires dense monitoring in advance. The 2016 Amatrice-Visso-Norcia (central Italy) earthquake cascade, captured by an unprecedented seismic network, provided a unique testing ground to image the preparatory phase of a large event. The crustal volume of the Norcia incipient fault was densely illuminated by seismic rays from more than 13,000 earthquakes that occurred within the 3 mo before the main shock nucleation. We performed seismic tomography in distinct time windows that revealed the precursory changes of elastic wave speed, signaling (1) the final locked state of the fault, and (2) the rapid fault-stiffness alterations near the hypocenter just a few weeks before the event. The results are the first instance where short-lived, hard-to-catch crustal properties shed light on evolving earthquake cascades.

Joint interpretation of seismic refraction tomography and electrical resistivity tomography by cluster analysis to detect buried cavities

In the last few years, the geophysical methods of seismic refraction tomography (SRT) and electrical resistivity tomography (ERT) are among the most used geophysical techniques for the reconstruction of subsoil geometries, for the investigation of underground cavities and also for the archaeological prospecting. However, the main disadvantage of each geophysical method is the difficulty of final interpretation of the data. In order to eliminate artifacts and generally improve the reliability and accuracy of geophysical interpretation, it is useful to perform a joint approach of different geophysical methods, also introducing the a priori information. In this work, it is shown the integrated study of seismic refraction tomography and electrical resistivity tomography techniques, the two geophysical methods are tested on both synthetic and real data and the integration of data is useful in detecting buried cavities and also evaluate their geometric characteristics. Likelihood parameters has been defined and tested, in order to help recognizing voids from other lithological structures. Finally, a statistical approach based on cluster analysis of the P-wave velocity, the density of the seismic rays and the electrical resistivity of the synthetic and experimental models was used. Multi-space cluster distribution maps were built, allowing to better define and interpret the anomalies of the subsoil.

The impact of the complex evaporites features’ on the seismic illumination of the underlying rocks: A case study in the Brazilian presalt

The oil industry utilizes seismic techniques to image the Earth's subsurface. The number of reflected rays generated at a single point that return to the receivers impacts directly the quality of seismic data and, in turn, seismic interpretation. Seismic illumination is used as a quality control measure to find poor coverage zones, which are commonly found in pre-salt reservoirs of eastern Brazilian marginal basins. The Santos Basin has a wide evaporite section featuring some complex structures such as overhangs, mini-basins, salt domes, salt tongues and walls more than 3 km thick. Those structures change the paths of the seismic rays and reduce the coverage of seismic imaging below the evaporite layer in some regions. In this study, a quantitative and qualitative illumination analysis of a pre-salt horizon was performed comparing four different velocity model scenarios and the behavior of the seismic ray path was analyzed. The results show that both stratification and the geometric structures existing in the salt layer have a direct impact on the imaging of the pre-salt reservoirs. Each factor influences imaging in a different way and, when associated, can generate regions of low illumination quality. Consequently, it is of extreme importance to obtain a suitable understanding of the upper salt layer as it will govern the ray paths towards the layers below the evaporites. Furthermore, the mapping of the anomalies generated by the stratification should also be considered in the creation of the velocity model.

A hybrid approach to compute seismic travel times in three‐dimensional tetrahedral meshes

ABSTRACTWe propose an optimized method to compute travel times for seismic inversion problems. It is a hybrid method combining several approaches to deal with travel time computation accuracy in unstructured meshes based on tetrahedral elementary cells. As in the linear travel time interpolation method, the proposed approach computes travel times using seismic ray paths. The method operates in two sequential steps: At a first stage, travel times are computed for all nodes of the mesh using a modified version of the shortest path method. The difference with the standard version is that additional secondary nodes (called tertiary nodes) are added temporarily around seismic sources in order to improve accuracy with a reasonable increase in computational cost. During the second step, the steepest travel time gradient method is used to trace back ray paths for each source–receiver pair. Travel times at each receiver are then recomputed using slowness values at the intersection points between the ray path and the traversed cells. A number of numerical tests with an array of different velocity models, mesh resolutions and mesh topologies have been carried out. These tests showed that an average relative error in the order of 0.1% can be achieved at a computational cost that is suitable for travel time inversion.

Magma accumulation beneath Santorini volcano, Greece, from P-wave tomography

AbstractDespite multidisciplinary evidence for crustal magma accumulation below Santorini volcano, Greece, the structure and melt content of the shallow magmatic system remain poorly constrained. We use three-dimensional (3-D) velocity models from tomographic inversions of active-source seismic P-wave travel times to identify a pronounced low-velocity anomaly (–21%) from 2.8 km to 5 km depth localized below the northern caldera basin. This anomaly is consistent with depth estimates of pre-eruptive storage and a recent inflation episode, supporting the interpretation of a shallow magma body that causes seismic attenuation and ray bending. A suite of synthetic tests shows that the geometry is well recovered while a range of melt contents (4%–13% to fully molten) are allowable. A thin mush region (2%–7% to 3%–10% melt) extends from the main magma body toward the northeast, observed as low velocities confined by tectono-magmatic lineaments. This anomaly terminates northwest of Kolumbo; little to no melt underlies the seamount from 3 to 5 km depth. These structural constraints suggest that crustal extension and edifice loads control the geometry of magma accumulation and emphasize that the shallow crust remains conducive to melt storage shortly after a caldera-forming eruption.

What can P-wave polarization direction data tell us about subsurface velocity structures?

SUMMARY As seismic waves propagate in the Earth, the directions of particle motions are affected by the media that they encounter, and thus seismic wave polarization direction carries the information on the media. So far there remains unclear about what can be inferred from the P-wave polarization direction data. For clarifying it, we discuss the mapping relation between polarization direction and velocity distribution. It is found that the velocity model cannot be derived uniquely from the polarization direction data. By analysing the relation between slowness vectors of the seismic ray at the source and the receiver, we find that relative velocity gradient is the physical quantity that describes the capability to deflect seismic rays in a continuous medium. The equation describing the relation between polarization direction and relative velocity gradient is given. For imaging relative velocity gradients, we derive the calculation formula for the partial derivative of polarization direction with respect to velocity gradient parameters. Synthetic experiments are conducted. The test results demonstrate that the absolute velocity model cannot be recovered from P-wave polarization direction data, but the relative velocity gradient model can. Polarization direction tomography gives a way to build gradient maps for the geometric characteristic of the subsurface velocity structures.