Crosshole seismic tomography including the anisotropy effect

Waveform tomography is commonly tested using numerically generated synthetic seismic data, before the method is applied to field seismic data. However, there are often noticeable differences between idealized synthetic data and real field data, and many factors in the field data, such as noise, irregular source/receiver geometry, affect the inversion solutions. For exploring the potential of reflection seismic waveform tomography, we presented a more realistic test than the synthetic data test, by applying it to physical modelling data, to reconstruct a laboratorial model with complex velocity variation. First, we provided a formulation of the perfectly matched layer absorbing boundary condition, associated with the second-order acoustic wave equation, in order to suppress artificial reflections from subsurface model boundaries in seismic waveform simulation and tomography. Then, we demonstrated the successful implementation of a layer-striping inversion scheme applicable to reflection seismic waveform tomography. Finally, we confirmed the effectiveness of frequency grouping, rather than a single frequency at each iteration, a strategy specifically for the frequency-domain waveform tomography.

Seismic bulletins, with trustworthy phase picks, origin times, and source locations are key for regional seismic studies, such as travel-time (TT) tomography, attenuation tomography, and anisotropy studies. To lay the groundwork for such studies in Israel, we revised the seismic bulletin of Israel and the surrounding area and obtained a trustworthy TT data set. From the earthquake and explosion bulletins of the Geophysical Institute of Israel, we compiled a starting data set of about 123,000 earthquakes and explosions that occurred during the past 40 yr. After screening out the poorly recorded events, we were left with a data set of ∼38,000 well-recorded events. We then revised the remaining data set in two consecutive steps. In the first, we reviewed and updated station metadata, including changes in station metadata parameters over time. In the second step, we jointly relocated a list of selected seismic events, using the Bayesian hierarchical location software package (BayesLoc) of Myers et al. (2007) that performs joint relocation of multiple events. We observed striking dissimilarities between the spatial distributions of the newly relocated catalog and the initial locations. Although the depth distribution of the starting catalog is trimodal with peaks at 0, 5, and 10 km, the distribution in this study is unimodal, with a broad peak between 7.5 and 12.5 km. By differencing the observed arrival times and the origin times obtained through relocation with BayesLoc, we obtained a revised TT database that consists of 261,336 Pg, 132,876 Pn, 114,816 Sg, and 60,394 Sn arrivals, from a set of 30,458 jointly relocated seismic sources. We compared prerevision and postrevision TTs as a function of epicentral distance and concluded that the revised data set contains far fewer outliers and inconsistencies than the original data set. The revised TT data set may be used for seismic studies, such as TT tomography, attenuation tomography, and anisotropy studies.

We argue in this paper for a hierarchical approach for Full Waveform Inversion that begins with the establishment of “starting models” with traveltime tomography, proceeds to Acoustic Waveform Tomography, and will eventually lead on to Anisotropic, Visco-Elastic Waveform Tomography inversions. With each step increasingly accurate physical descriptions are invoked, more of the information content of the data is used, and the final results are correspondingly enhanced. The ultimate goal is “Full” Waveform Inversion, in which the full (anisotropic, 3D, visco-elastic) physics of the earth is incorporated and all parts of the waveform are utilized. Within the framework of Acoustic Waveform Tomography itself, we formulate 2 alternative approaches to the Laplace-Fourier domain formulation of waveform inversion: i) A conventional L2 misfit function, and ii) phase-only misfits derived from the logarithmic wavefield residuals. We apply these approaches to two real 2D data examples: Marine OBS data from the Nankai subduction zone, and land surface seismic data from the Canadian Foothills. We find in both data examples that the phase-only misfit outperforms the conventional L2 approach, and we suggest that this is due to increased ability of the phase-only approach to recover the low-wavenumber components of the model. We also find that a logarithmic phase-amplitude approach that provides a marginal improvement over the phase-only approach. Finally, we use point scattering responses and chequerboard tests to demonstrate the artefacts introduced by the use of an acoustic approach.

Detecting and imaging dense non-aqueous phase liquids (DNAPLs) in the subsurface is a challenging problem that is of interest to the Department of Energy, Department of Defense, and several local and state agencies. Current DNAPL detection techniques, such as Direct Push Probe Technologies (DPT) and In-Situ Tracers (IST) have risks and limitations. Cross-well radar (CWR) is a radar-based geophysical technique with low invasiveness for real time monitoring of DNAPLs. This technique uses electromagnetic waves transmitted and received through antennas in the subsurface. The computerized tomography is applied to a pilot-scale experimental facility constructed by the authors (referred as SoilBED) in this research. Cross-tomography data at multiple depths and locations were collected to simulate a 1/100 scaled contaminated soil problem. The experimental results are converted to the time domain signals. The resulting signal travel times are compared with the simulated results by FDTD. The transformed signals are used to prepare signal intensity and travel time tomograms of the SoilBED cross-sectional slices, for background and scattered field to study the signature of the scatterers. The results are able to visualize the dielectric objects in the saturated soil. Introduction Cross-well radar is a potential technique for DNAPL detection. Most of the research and field works done in Cross-well radar and other radar based geophysical techniques deal with a single frequency. Selecting one single frequency may be effective for one non-dispersive medium, but perhaps not for other non-dispersive or most dispersive media. This research covers a broadband of frequency from 400 MHz to 2.2 GHz. Having broadband frequency data helps to either use the entire frequency range and work in frequency domain or select a single frequency or narrower frequency range that seems to have enough detection potential out of the entire observed frequency response. It is important to consider cost and practicality of the system. Establishing practical and cost effective systems require strong theoretical and laboratory background conditions. Therefore, both theoretical and experimental aspects of the research are studied in parallel, to evaluate the feasibility of DNAPL detection using CWR and to study wave propagation behavior and forward problem. Collecting repeatable and reproducible frequency response data is a very challenging process. This difficultly achieved data can be manipulated first or directly used towards image processing. As any research subject, conducting both experimental and theoretical approaches is the most desirable way to advance fundamental knowledge in the field of interest. Frequency response and arrival time can be converted to each other. Thus, frequency domain results can be transformed to time domain and used towards techniques such as travel time tomography and vice versa. The experimental frequency domain results can be converted to CTF (Channel Transfer Functions) and compared to the simulated frequency domain results by the first order Born approximation (Farid 2004, Farid et al. 2005, and Zhan 2005). The FDTD simulation code (Farid 2002 and 2003) is available as the theoretical tool to perform the comparison in time domain and study travel time. If the experimental results are converted to time domain and transformed to travel time, they can be compared to the FDTD simulated results. Travel Time Computed by FDTD To establish a better knowledge of wave propagation in soil, both dry and saturated backgrounds are studied. The sandy soil in the SoilBED facility with 3.88% moisture content is referred to as dry soil.

Ray tracing is avoidable in seismic traveltime tomography. In this study, seismic traveltime tomography is reformulated as an eikonal equation‐constrained optimization problem solved by the ray‐free adjoint‐state method. The resultant approach is called adjoint‐state traveltime tomography. For completeness, an eikonal equation‐based earthquake location method is developed to locate the hypocenters of earthquakes when necessary. The multiple‐grid model parameterization is adopted to discretize the relative slowness perturbation . The step‐size‐controlled gradient descent method, which has an effect of damping regularization, is used to find optimal earthquake hypocenters and velocity models. The performances of the earthquake location and adjoint‐state traveltime tomography methods are tested in the source area of the 2020 Mw 4.9 Anza earthquake, a seismologically active place in southern California. The obtained high‐resolution P‐wave velocity model demonstrates that the 2020 Mw 4.9 Anza earthquake and other historic moderate‐sized Anza earthquakes occurred near or at the boundaries of low VP rocks, a typical seismogenic environment for moderate to large crustal earthquakes. Meanwhile, the source zone of the Anza earthquakes is characterized by high VP/VS, indicating the possible existence of crustal fluids as well as the critical role of fluids in the occurrence of those moderated‐sized Anza earthquakes. On the whole, the well‐performed eikonal equation‐based earthquake location method and adjoint‐state traveltime tomography method provide competent and attractive tools for hypocenter and tomographic inversions.

Seismic attenuation tomography and its application to rock mass evaluation

Time shift estimation is a key issue in many areas of seismic exploration such as time-lapse studies and traveltime tomography. Automated estimation is useful because it enables applications such as time strain analysis and automated velocity model building. A commonly used automated estimation algorithm is the windowed crosscorrelation. Unfortunately, this method can be inaccurate in areas where time shifts vary significantly over short intervals. Algorithms based on mismatch minimization have been proposed to remedy this problem. All the mentioned methods suffer from cycle skipping in the face of large time shifts. I, therefore, used the dynamic time warping (DTW) algorithm as an alternative method to estimate time shifts. I claimed the time shift estimates produced by DTW are of higher quality than those produced by the previous methods, especially with regard to the cycle-skipping problem. I found that DTW always is able to perfectly restore time shifts applied to traces when the changes are purely temporal. I investigated its accuracy when amplitude changes are present, and found that the algorithm is well suited for time-lapse analysis. Finally, we determined that DTW is resilient against cycle-skipping errors, in the presence of large time shifts and with significant changes in the data, and thus a candidate for use in automated velocity model-building methods, such as waveform tomography or full-waveform inversion.

A fast marching algorithm for the factored eikonal equation

The current rapid growth in the number of seismometers available to the research community, combined with increasing computer power, will allow improvement in the type and quality of seismic images of the crust and lithosphere. An example of improved imaging capability is the inversion of the full seismic waveform, rather than solely travel times, in controlled‐source surveys (seismic refraction or reflection using human‐induced ground shaking).At the 12th Deep Seismic Methods workshop in 2003, sponsored by the International Association of Seismology and Physics of the Earth's Interior (IASPEI), an analysis of computer‐generated data exemplified the potential of increased source and station density. The synthetic seismic data set was generated from a geologic model that includes large‐, medium‐, and small‐scale stochastic variation. The source and seismometer spacing mimic imminent community capabilities. The Earth model was kept secret, and the data were made available for analysis (http://crust. geol.vt.edu/hole/ccss/). Figure 1 illustrates the results of blind travel time and waveform tomography applied to the data.The images, described in more detail below, illustrate an excellent match to the true Earth model.

<p>The Ordos block is located on the west side of the North China Carton, adjacent to the northeastern part of the Tibetan Plateau. Affected by two tectonic movements, Ordos block internal structure remains relatively stable structure, but surrounded by active tectonic belts. With the development of the second and third part of the “China Seismological Science Array”, the distribution of seismic observation stations in Ordos region has been greatly improved. This study will use the new seismic observation data of "Array III", and combined with the phase observation data of "Array II" to form a more complete seismic phase travel time data set. The regional seismic body-wave travel time tomography will figure out a more reliable three-dimensional velocity structure of P waves in Ordos.</p><p>Our study area spans from 32°N to 42°N and 108°E to 114°E , which includes the Ordos block and its adjacent structures . The seismic data we used for inversion were recorded by 1244 stations including: 198 permanent stations and 1043 temporary stations (ChinArray II and III), from November 2013 to August 2017. After manual labeled the seismic phase, we select events with more than ten phase records of individual seismic events. The epicentral distance is less than 200km. Finally, we obtained about 22,500 phase records of 1882 local seismic events.</p><p>The preliminary results are consistent with previous studies and surface structures of a wide range of velocity distributions. However, in the middle-upper crust under the Liupan Mountain west, the low-speed anomaly extending downward is shown, which may be caused by the shallow crustal damage caused with the continuous eastward compression of asthenosphere in the northeastern margin of the Qinghai-Tibet Plateau during the Cenozoic. It is worth noting that there is an EW-trending low-velocity zone under the Dingbian-Suide fault beneath the Ordos Basin, with a depth form lower crust to 50 km in upper mantel. This low-velocity anomaly divides the high-speed disturbance in the Ordos block into two parts，indicate the depth of the fault can reach the upper mantel. In the Taihang Mountains in the west of the study area, low-velocity anomalies extending to the upper layer of the mantle are shown. We initially believe that this anomaly is related to the volcanic thermal motion that once existed on the area.</p>

The structure and evolution of the lithosphere and underlying mantle control tectonics, magmatism, and natural resource development and distribution. Seismic tomography offers essential information on mantle temperature, the thickness and strength of the lithosphere, and the convection below it. Here, we report on the completion of a series of continent- and plate-scale waveform tomography studies using the automated multimode inversion of surface and S-wave waveforms. Together, the models cover all of the Earth’s continents and their surroundings and reveal a fascinating diversity of structures, while also indicating common mechanisms of lithospheric dynamics and evolution. Each model used global data coverage that was maximised witin the hemisphere centred at the continent. Structural information was extracted from Rayleigh waves and S and multiple S waves on over a million vertical-component-seismogram waveform fits. The effects of errors were minimised by statistical and targeted outlier analyses and the removal of the least mutually consistent data. The models advance the resolution of the imaging compared to the state of the art at the scale of the continents.Cratonic lithosphere on all continents shows complexities and fragmentation. The lithosphere beneath most diamondiferous kimberlites—originally emplaced on thick cratonic lithosphere—is still thick at present. Relatively low seismic velocities at kimberlite locations are indicative of craton-lithosphere thinning and are detected at multiple locations in all continents with known kimberlites. Cenozoic basalts are found exclusively where the lithosphere is observed to be thin. Beneath some of the volcanic regions, low-velocity anomalies extend deep into the mantle, consistent with deep-mantle upwellings feeding the magmatism. Sediment-hosted metal deposits tend to be located near contrasts in the thickness of the lithosphere, including both craton boundaries and other substantial heterogeneities. Intraplate seismicity is controlled by plate-boundary stresses and by the lateral variations of the lithospheric thickness and strength. Areas with relatively thin lithosphere localise deformation and seismicity, provided that sufficient tectonic stress is transmitted into the plate interior.ReferencesChagas de Melo, B, Lebedev, S, Celli, NL, et al., 2025. The lithosphere of South America from seismic tomography: Structure, evolution, and control on tectonics and magmatism, Gondwana Research, GR Focus - Invited Review, 138, 139–167. doi:10.1016/j.gr.2024.10.012.Chua, EL, Lebedev, S. Waveform tomography of the Antarctic Plate, Geophysical Journal International, submitted.Celli, NL, Lebedev, S, Schaeffer, AJ, Gaina, C, 2021. The tilted Iceland Plume and its effect on the North Atlantic evolution and magmatism, Earth and Planetary Science Letters 569, doi:10.1016/j.epsl.2021.117048.Celli, NL, Lebedev, S, Schaeffer, AJ, et al., 2020. The upper mantle beneath the South Atlantic Ocean, South America and Africa from waveform tomography with massive data sets, Geophysical Journal International 221, 178–204. doi:10.1093/gji/ggz574.de Laat, JI, Lebedev, S, Celli, NL, et al., 2023. Structure and evolution of the Australian plate and underlying upper mantle from waveform tomography with massive data sets, Geophysical Journal International 234, 153–189. doi:10.1093/gji/ggad062Dou, H, Xu, Y, Lebedev, S, et al., 2024. The upper mantle beneath Asia from seismic tomography, with inferences for the mechanisms of tectonics, seismicity, and magmatism, Earth-Science Reviews, 255, doi:10.1016/j.earscirev.2024.10484

We investigate issues of convergence, resolution, and nonlinearity related to the feasibility of adjoint tomography in regional and global tomography and exploration geophysics. Most current methods of adjoint tomography, whether based on adjoint methods or other formulations, suffer from slow convergence in that only the gradient (not the Hessian) is readily available for computing model updates. As an alternative to working with the unpreconditioned gradients, we examine the speed-up offered by various preconditioners that can be computed in the framework of adjoint methods. We show that each preconditioner bears some similarity to the Hessian, thus motivating and justifying its use for accelerating convergence. Next, we examine the role of the Hessian in resolution analysis. Recalling that the action of the Hessian on an arbitrary model perturbation relates to the classical point spread function concept, we introduce a scalar quantity termed the average eigenvalue that provides a good overall representation of resolution. Whereas a pointspread function reveals the orientation of misfit contours, the average eigenvalue describes the sharpness of the misfit function along the direction of the chosen model perturbation. Finally, we provide an example in which we directly compare the results of travel time and waveform tomography, illustrating the resolution limits of the former and the nonlinearity pitfalls of the latter.

Traveltime tomography is applied to investigate seismic structures of the Earth's subsurface. An accurate tomographic velocity model is important for a high-resolution waveform velocity building and its availability is one of the main components to mitigate the nonlinear inverse problem. We present a new methodology of obtaining velocity models for traveltime tomography studies. We found a way to get a highly accurate first-arrival traveltime tomography in combination with global optimization. The role of global optimization is twofold: to find initial solutions that are close to ‘truth’, and to guide tomographic inversion towards a geologically consistent model that explains the data. The main advantage of our workflow is a data-driven approach avoiding the use of a conventional layer-based parameterization and incorporation of manual interpretations into the velocity model. To date, a few geophysical studies have been focused on developing data-driven and a labour non-intensive regional tomographic velocity model building workflow. In our study, we present the tomographic velocity model building workflow as a combination of first-arrival traveltime tomography and global optimization. Global optimization allows to search for velocity parameters and depth to interfaces in the larger search area with a higher chance of convergence. After defining the geometry of main layers and general velocity trends, traveltime tomography with a bi-cubic B-spline model parameterization can be fitted to further update the velocity model. Our approach allows obtaining a highly accurate velocity model which can be used for seismic depth migration and as a starting model for a FWI seismic imaging. The workflow is developed and applied to synthetic and field regional seismic datasets. The developed methodology is applied for a shallow seismic engineering data and regional Ocean Bottom Seismic data. We identify four key components that lead to building an accurate tomographic velocity model: (i) understanding prominent horizons and possible velocity distribution of a layer within the study area. (ii) Performing ray penetration test to define offset ranges which carry the velocity information for the defined layers. (iii) Determining inversion schema to a perform global search for the velocity trends and major boundaries, and a local search to update lateral velocity variation. (iv) Iteratively update a set of defined layers (i.e., sediment, igneous crust and basement) in a top-down manner. 

In this study the potential of an inversion approach based on hydraulic travel time and hydraulic attenuation tomography was assessed. Both hydraulic travel time and hydraulic attenuation tomography are based on the transformation of the transient groundwater flow equation into the eikonal equation using an asymptotic approach. The eikonal equation allows the calculation of pressure propagation and attenuation along trajectories, which is computationally efficient. The attenuation and travel time‐based inversion approaches are naturally complementary: hydraulic travel times are determined by the hydraulic diffusivity, a combination of hydraulic conductivity and specific storage, whereas the attenuation is determined solely by specific storage. The potential of our hydraulic tomographical approach was investigated at a well‐characterized sand and gravel aquifer located in the Leine River valley near Göttingen, Germany. The database for the hydraulic inversion consists of 392 cross‐well slug interference tests performed between five wells, in which the positions of the sources (injection ports) and the receivers (observation ports), isolated with double packer systems, were varied between tests. The results have shown that the combination of hydraulic travel time and hydraulic attenuation tomography allows the reconstruction of the diffusivity and storage distribution in two and three dimensions with a resolution and accuracy superior to that possible with type curve analysis.

We use a machine learning-based tomography method to obtain high-resolution subsurface geophysical structure in Long Beach, CA, from seismic noise recorded on a “large-N” array with 5204 geophones (~13.5 million travel times). This method, called locally sparse travel time tomography (LST) uses unsupervised machine learning to exploit the dense sampling obtained by ambient noise processing on large arrays. Dense sampling permits the LST method to learn directly from the data a dictionary of local, or small-scale, geophysical features. The features are the small scale patterns of Earth structure most relevant to the given tomographic imaging scenario. Using LST, we obtain a high-resolution 1 Hz Rayleigh wave phase speed map of Long Beach. Among the geophysical features shown in the map, the important Silverado aquifer is well isolated relative to previous surface wave tomography studies. Our results show promise for LST in obtaining detailed geophysical structure in travel time tomography studies.