Seismic Traveltime Tomography Research Articles

Smoothness: The key factor in well-log information-assisted PINNtomo

Seismic traveltime tomography is popularly used to provide a smooth background model for subsequent depth migration or full waveform inversion (FWI). The core of the traveltime tomography process requires an eikonal solver for seismic traveltime calculation. Most of the conventional eikonal solvers, suffer from first-order convergence errors and difficulties in dealing with irregular topography, which will lead to inaccurate seismic tomography results. To solve these issues, physics-informed neural networks (PINNs) have been developed and demonstrated remarkable success in simultaneous traveltime calculation and velocity construction, which is referred to as PINNtomo. It combines the physics constraint (eikonal equation) and data constraint (traveltime fitting) to form the loss function of the networks. Though PINNtomo can manage to recover a smooth background velocity, it will fail in the inversion for a complicated high-resolution velocity model. We take advantage of the well-log velocity as a prior information to provide an additional model constraint in the loss function to improve the performance of PINNtomo. We found that the well-log information will provide high-resolution components in the inverted velocity model. However, it will destroy the lateral structure of the true model. We propose to use a smoothed version of the well-log velocity as the prior model constraint in PINNtomo. As a result, the accuracy of the inverted background velocity will be improved, which will be beneficial to the subsequent FWI.

A matching pursuit approach to the geophysical inverse problem of seismic travel time tomography under the ray theory approximation

Summary Seismic travel time tomography is a geophysical imaging method to infer the 3-D interior structure of the solid Earth. Most commonly formulated as a linearized inverse problem, it maps differences between observed and expected wave travel times to interior regions where waves propagate faster or slower than the expected average. The Earth’s interior is typically parametrized by a single kind of localized basis function. Here we present an alternative approach that uses matching pursuits on large dictionaries of basis functions. Within the past decade the (Learning) Inverse Problem Matching Pursuits ((L)IPMPs) have been developed. They combine global and local trial functions. An approximation is built in a so-called best basis, chosen iteratively from an intentionally overcomplete set or dictionary. In each iteration, the choice for the next best basis element reduces the Tikhonov–Phillips functional. This is in contrast to classical methods that use either global or local basis functions. The LIPMPs have proven their applicability in inverse problems like the downward continuation of the gravitational potential as well as the MEG-/EEG-problem from medical imaging. Here, we remodel the Learning Regularized Functional Matching Pursuit (LRFMP), which is one of the LIPMPs, for travel time tomography in a ray theoretical setting. In particular, we introduce the operator, some possible trial functions and the regularization. We show a numerical proof of concept for artificial travel time delays obtained from a contrived model for velocity differences. The corresponding code is available online.

Bayesian Inversion, Uncertainty Analysis and Interrogation Using Boosting Variational Inference

AbstractGeoscientists use observed data to estimate properties of the Earth's interior. This often requires non‐linear inverse problems to be solved and uncertainties to be estimated. Bayesian inference solves inverse problems under a probabilistic framework, in which uncertainty is represented by a so‐called posterior probability distribution. Recently, variational inference has emerged as an efficient method to estimate Bayesian solutions. By seeking the closest approximation to the posterior distribution within any chosen family of distributions, variational inference yields a fully probabilistic solution. It is important to define expressive variational families so that the posterior distribution can be represented accurately. We introduce boosting variational inference (BVI) as a computationally efficient means to construct a flexible approximating family comprising all possible finite mixtures of simpler component distributions. We use Gaussian mixture components due to their fully parametric nature and the ease with which they can be optimized. We apply BVI to seismic travel time tomography and full waveform inversion, comparing its performance with other methods of solution. The results demonstrate that BVI achieves reasonable efficiency and accuracy while enabling the construction of a fully analytic expression for the posterior distribution. Samples that represent major components of uncertainty in the solution can be obtained analytically from each mixture component. We demonstrate that these samples can be used to solve an interrogation problem: to assess the size of a subsurface target structure. To the best of our knowledge, this is the first method in geophysics that provides both analytic and reasonably accurate probabilistic solutions to fully non‐linear, high‐dimensional Bayesian full waveform inversion problems.

Travel-time seismic tomography for the seismic stratigraphic interpretation of the crust around the San Rafael knickpoint at Coca River, Ecuador

Travel-time seismic tomography for the seismic stratigraphic interpretation of the crust around the San Rafael knickpoint at Coca River, Ecuador

Three-dimensional velocity structure of a coal mine revealed by induced microseismic traveltime data

Early detection of rockburst risk areas is a prerequisite for ensuring safe production in coal mines. A reliable and efficient detection method is essential for identifying stress distribution and disturbances in coal and rock masses during mining. In this study, we employed seismic traveltime tomography to study the three-dimensional velocity structure of a continuous mining and excavation area of a coal mine by inverting a number of P-wave arrival times from the induced micro-earthquakes (M < 2.0). The response characteristics and stress evolution of the inversion areas at different mining stages were investigated by examining the obtained time-lapse velocity structure. The microseismic (MS) dataset was partitioned into two clusters in a temporal sequence, with each cluster containing 367 relocated events. A series of checkerboard and restoring resolution tests, as well as field investigations, confirmed the resolution, robustness, and reliability of the tomographic results. In addition, absolute velocity change patterns were proposed to assess and predict rockburst risk areas caused by mining activities. The results showed that high-velocity zones could be attributed to roof subsidence, stress transfer after coal and rock mass rupture, and areas with special geological conditions in isolated working face situations. The coupled effects of stress evolution and geological conditions may lead to more pronounced stress concentrations during the mining process. Furthermore, induced events occur not only in areas between high and low stress but also near the isoline edges of pattern changes. Our established velocity patterns have significant potential for predicting hazardous areas with a timeliness of up to 10 days.

High-resolution lithospheric structure of continental China from joint inversion of surface wave and gravity data

The lithospheric structure of continental China has been previously determined by seismic travel time tomography, surface wave tomography, and joint inversion of body wave and surface wave data. However, due to the inherent limitations of seismic data, the lithospheric structure of continental China is still not well resolved in the shallow part and in some regions where the station coverage is relatively sparse. In this study, we aim at improving the lithospheric structure by joint inversion of seismic surface wave data and satellite gravity data to take advantage of the uniform distribution and complementary strength of the gravity data. The empirical relationship between velocity and density is used as a bridge for joint inversion of surface wave and gravity data. The joint inversion shear-wave velocity Vs (density) model, named as USTClitho1.0g, can fit both surface wave and gravity data well. This high-resolution Vs model can better fit the active airgun source seismic arrival times and better delineate some features in continental China, such as the velocity contrast across the north–south gravity lineament (NSGL), the lithosphere thinning in eastern China, the crustal footprint of the Hainan mantle plume, the likely magma chamber beneath volcanos in northeast China, the middle-lower crust low velocity layer beneath the Tibetan plateau, and the tearing of subducted Indian mantle lithosphere. Our joint inversion Vs model can provide a reference model for geosciences in continental China and surrounding areas.

2D and 3D modelling of Molucca Sea double subduction zone based on primary waves travel time seismic tomography inversion

Molucca Sea is one of the areas with a high level of seismicity in Indonesia due to its location between three plates. The activities of them forming a quite unique subduction zone. The Eurasian Plate (Sangihe Microplate) move to the East and The Philippine Sea Plate (Halmahera Microplate) continues to Westward caused The Molucca Sea Microplate to be pressed, and subducts in two opposite directions modelat the same time forming an inverted “U”. This research used 3663 hypocenters and 9 BMKG stations in the last 5 years (2016 - 2021) with a magnitude of 3 - 8 Mw and a depth of 3 - 630 kilometers. Seismic inversion of travel time tomography is used to produce 2D models of the subduction zone around The Molucca Sea. The distribution of high ΔVp values is ranging from 3.5% to 4% is associated with a solid medium, like The Sangihe Plate, Halmahera Plate, and The Molucca Sea Plate subduction slab. While the low ΔVp is ranging from 1.8% to 2% is associated with the destruction zone and the presence of thermal fluids such as magma or partial melting (the presence of volcanoes). The 3D model is based on the Gaussian Process Regression principle using magnitude and depth data as initial model parameters. Plotting of the thrust faults at the front of the two arcs, namely The West Sangihe Thrust Fault and The East Halmahera Thrust Fault. It can be seen that the subduction under The Sangihe Microplate is deeper (627.2 kilometers) has the average subduction angle of 45.3° compared to the subduction under The Halmahera Microplate (280.0 kilometers) has the average subduction angle of 35.8°

Deformation of the NE Tibetan Plateau revealed by velocity and azimuthal anisotropy structures

Because of the complex interaction among the Tibetan Plateau, the Alxa block and the Ordos block, the deformation mechanisms of the northeastern (NE) Tibetan Plateau have been controversial for decades. To shed light on these mechanisms, we obtain high-resolution 3-D images of the P-wave velocity and azimuthal anisotropy structures beneath the NE Tibetan Plateau and adjacent area. We invert high-quality manually picked first arrivals using an eikonal equation-based seismic traveltime tomography method. Tomographic results mainly present two features. First, low-velocity anomalies are evident in the lower crust of the Qilian block. In addition, weak azimuthal anisotropy with NW–SE-oriented fast wave directions (FVDs) is detected in the lower crust. These characteristics indicate compressional deformation of the Qilian block. Second, the lower crust and uppermost mantle below the eastern Qilian orogenic belt display high-velocity anomalies. These high-velocity anomalies are collected with the high-velocity body beneath the Alxa block. This feature may imply that the Alxa block is underthrusting beneath the eastern Qilian orogenic belt. Our tomographic results show that the 1920 Ms 8.5 Haiyuan earthquake occurred near the boundary of the high-velocity anomaly, where the FVDs change dramatically. In addition, we observe that the 1927 Ms 7.9 Gulang earthquake occurred in the normal upper crust and above the strong low-velocity anomalies. From our tomographic results, we propose that two different seismogenic structures were associated with these two great earthquakes.

Phase-preserving theory and its linearization approximation for forward scattering field of scalar acoustic wave equation

The conventional wave-equation linearization methods, such as the first-order Born or Rytov approximation, always implicitly imply a weak-scattering assumption, making it valid only for weak perturbation models. To extend the wave-equation linearization theory to strong perturbation models, we consider a scenario that the reference model is smooth within the scale of the incident wave length, and propose a phase-preserving method which can predict the phase perturbation of forward scattering wave field. First, we introduce the WKBJ approximation to the scattered- and incident wave fields so that the integral of the unknown solution (i.e. the scattered field) in the nonlinear Ricatti integral equation can be replaced by the integral of scattering-angle and model perturbation, yielding an explicit expression of the scattered field. Theoretical derivation shows that the proposed phase-preserving method can accurately predict the phase-perturbation of forward scattered wave field regardless of the strength of velocity perturbations for one-dimensional wave propagation problem. To apply the phase-preserving approximation to the inverse problem, we further consider a scenario of small-angle forward propagation. In this case, the phase-preserving approximation can be linearized by neglecting the influence of scattering angles, leading to a linear relation between the scattered field and the model perturbation, which we refer to as the generalized Rytov approximation. Numerical experiments demonstrate that the generalized Rytov approximation can predict the phase perturbation of the scattered field with higher accuracy for small-angle forward propagation, and is suitable for strong model perturbations. The generalized Rytov approximation extends the validity and the scope of application of the traditional Rytov approximation. In specific application fields such as the seismic traveltime tomography or medical ultrasonic transmission imaging, a new traveltime/phase sensitivity kernel can be derived by replacing the conventional Rytov approximation with the proposed method, which can increase the inversion accuracy and speed up the convergence.

Hyperparameter estimation using a resolution matrix for Bayesian sensing

Bayesian image processing has progressively increased in importance in various fields of the natural sciences. It utilizes prior knowledge and forward models of the observational processes through Bayes’ theorem, enabling the accurate estimation of model parameters that represent the physical quantities of the target. Moreover, using hyperparameter estimation, we can determine the hidden physical parameters that govern the processes in and the structure of the target and sensing systems, such as the spatial continuity of the model parameters and the magnitude of the observational noise. Such a general framework, which uses Bayesian estimation to understand the essential physics of a target system, can be called ‘Bayesian sensing’. This paper discusses the physical meaning of and the mechanism underlying Bayesian sensing using the concept of resolution in the spatial-inversion problem. The spatial resolution of the model parameters can be mapped using a resolution matrix, more rigorously, a model resolution matrix defined as a linear mapping from the true model parameters to the recovered model parameters. We formulate the resolution matrix for Bayesian image processing and also show that in terms of resolution, the optimal hyperparameters are obtained from internally consistent equations that connect the estimated optimal hyperparameters with the actual ones calculated from the estimated model parameters. In addition, we show the equivalence of the internally consistent equations to the expectation-maximization (EM) algorithm and formulate the confidence intervals for the estimated hyperparameters, which indicate the general convergence of the hyperparameter estimates. We also show the effectiveness of the proposed method by performing synthetic numerical tests for two inversion-problem settings: linear travel-time seismic tomography and image deblurring. The resulting equations can contribute to understanding the hidden physical processes in and the structure of the target and observation systems for various problems.

Rapid magma ascent beneath La Palma revealed by seismic tomography

For the first time, we obtained high-resolution images of Earth's interior of the La Palma volcanic eruption that occurred in 2021 derived during the eruptive process. We present evidence of a rapid magmatic rise from the base of the oceanic crust under the island to produce an eruption that was active for 85 days. This eruption is interpreted as a very accelerated and energetic process. We used data from 11,349 earthquakes to perform travel-time seismic tomography. We present high-precision earthquake relocations and 3D distributions of P and S-wave velocities highlighting the geometry of magma sources. We identified three distinct structures: (1) a shallow localised region (< 3 km) of hydrothermal alteration; (2) spatially extensive, consolidated, oceanic crust extending to 10 km depth and; (3) a large sub-crustal magma-filled rock volume intrusion extending from 7 to 25 km depth. Our results suggest that this large magma reservoir feeds the La Palma eruption continuously. Prior to eruption onset, magma ascended from 10 km depth to the surface in less than 7 days. In the upper 3 km, melt migration is along the western contact between consolidated oceanic crust and altered hydrothermal material.

Waveform Energy Focusing Tomography With Passive Seismic Sources

By taking advantage of the information carried by the entire seismic wavefield, Full Waveform Inversion (FWI) is able to yield higher resolution subsurface velocity models than seismic traveltime tomography. However, FWI heavily relies on the knowledge of source information and good initial models, and could be easily trapped into local minima caused by cycle skipping issue because of its high nonlinearity. To mitigate these issues in FWI, we propose a novel method called Waveform Energy Focusing Tomography (WEFT) for passive seismic sources. Unlike conventional FWI, WEFT back-propagates the seismic records directly instead of the data residuals, and updates the velocity models by maximizing the stacking energy for all the moment tensor components from back-propagated wavefields around the sources. Therefore, except for source locations and origin times, WEFT does not require other source attributes in advance for the inversion. Since WEFT does not aim at fitting synthetic and observed waveforms, it has lower nonlinearity and is less prone to the cycle skipping issue compared to FWI. For the proof of concept, we have validated WEFT using several 2D synthetic tests to show it is less affected by inaccurate source locations and data noise. These advantages render WEFT more applicable for tomography using passive seismic sources when the source information is generally not accurately known. Although the inverted model from WEFT is inevitably influenced by the source distribution as well as its radiation patterns, and its resolution is likely lower than that of FWI, it can act as an intermediate step between traveltime tomography and FWI by providing a more reliable and accurate velocity model for the latter.

Geophysical Imaging of the Critical Zone along the Eastern Betic Shear Zone (EBSZ), SE Iberian Peninsula

The critical zone (CZ) represents the most-shallow subsurface, where the bio-, hydro-, and geospheres interact with anthropogenic activity. To characterize the thickness and lateral variations of the CZ, here we focus on the Eastern Betic Shear Zone (EBSZ), one of the most tectonically active regions in the Iberian Peninsula. Within the EBSZ, the Guadalentín Depression is a highly populated area with intensive agricultural activity, where the characterization of the CZ would provide valuable assets for land use management and seismic hazard assessments. To achieve this, we have conducted an interdisciplinary geophysical study along the eastern border of the Guadalentín Depression to characterize the CZ and the architecture of the shallow subsurface. The datasets used include Electrical Resistivity Tomography (ERT), first-arrival travel time seismic tomography, and multichannel analysis of surface waves (MASW). The geophysical datasets combined help to constrain the high-resolution structure of the subsurface and image active fault systems along four transects. The resulting geophysical models have allowed us to interpret the first ~150 m of the subsurface and has revealed: (i) the variable thickness of the CZ; (ii) the CZ relationship between the fault zone and topographic slope; and (iii) the differences in CZ thickness associated with the geological units. Our results provide a method for studying the shallow subsurface of active faults, complementing previous geological models based on paleo-seismological trenches, and can be used to improve the CZ assessment of tectonically active regions.

Crustal structure and intraplate seismicity in Nordland, Northern Norway: insight from seismic tomography

SUMMARY The Nordland region, Northern Norway, situated in an intraplate continental setting, has the highest seismicity rate in mainland Norway. However, the exact cause of seismicity in this region is still debated. Better understanding of factors that influence the seismicity in Nordland can help increase knowledge of intraplate seismicity in general. Here, we address this problem with the aid of a new high-resolution 3-D VP and VP/VS ratio images of the crust in Nordland using seismic traveltime tomography. These images show the existence of a localized, 10–15 km Moho step that runs parallel to the coast. The north–south extent of this step coincides with the region that exhibits the highest rates of seismicity. Focal mechanisms of selected earthquakes computed in this study are dominated by normal and oblique-normal, indicating a coast-perpendicular extension. The coast-perpendicular extensional stress regime deviates from the regional compression imposed by the ridge push from the North Atlantic. This deviation is thought to stem from the additional interference with local flexural stress caused by sediment redistribution and glacial isostatic adjustment, and possibly exacerbated by gravitational potential energy stress associated with the Moho step. The deformation due to the extensional regime is localized on pre-existing faults and fractures along the coastline. The tomography result shows that two distinct seismic swarms occurred in the coastal area with low VP and variable VP/VS ratio anomalies, pointing towards fractured crust and possibly the presence of fluids. The existence of fluids here can change the differential stress and promote seismic rupture.

MFAST: A MATLAB toolbox for ocean bottom seismometer refraction first-arrival traveltime tomography

First-arrival seismic traveltime tomography (FAST) is a well-established technique to estimate subsurface velocity structures. Although several existing open-source packages are available for first-arrival traveltime tomography, most were written in compiled languages and lack sufficient extendibility for new algorithms and functionalities. In this work, we develop an open-source, self-contained FAST package based on MATLAB, one of the most popular interpreted scientific programming languages, with a focus on ocean bottom seismometer refraction traveltime tomography. Our package contains a complete traveltime tomography workflow, including ray-tracing-based first-arrival traveltime computation, linearized inversion, quality control, and high-quality visualization. We design the package as a modular toolbox, making it convenient to integrate new algorithms and functionalities as needed. At the current stage, our package is most efficient for performing FAST for two-dimensional ocean bottom seismometer surveys. We demonstrate the efficacy and accuracy of our package by using a synthetic data example based on a modified Marmousi model.

On the Potential of 3D Transdimensional Surface Wave Tomography for Geothermal Prospecting of the Reykjanes Peninsula

Seismic travel time tomography using surface waves is an effective tool for three-dimensional crustal imaging. Historically, these surface waves are the result of active seismic sources or earthquakes. More recently, however, surface waves retrieved through the application of seismic interferometry have also been exploited. Conventionally, two-step inversion algorithms are employed to solve the tomographic inverse problem. That is, a first inversion results in frequency-dependent, two-dimensional maps of phase velocity, which then serve as input for a series of independent, one-dimensional frequency-to-depth inversions. As such, a set of localized depth-dependent velocity profiles are obtained at the surface points. Stitching these separate profiles together subsequently yields a three-dimensional velocity model. Relatively recently, a one-step three-dimensional non-linear tomographic algorithm has been proposed. The algorithm is rooted in a Bayesian framework using Markov chains with reversible jumps, and is referred to as transdimensional tomography. Specifically, the three-dimensional velocity field is parameterized by means of a polyhedral Voronoi tessellation. In this study, we investigate the potential of this algorithm for the purpose of recovering the three-dimensional surface-wave-velocity structure from ambient noise recorded on and around the Reykjanes Peninsula, southwest Iceland. To that end, we design a number of synthetic tests that take into account the station configuration of the Reykjanes seismic network. We find that the algorithm is able to recover the 3D velocity structure at various scales in areas where station density is high. In addition, we find that the standard deviation of the recovered velocities is low in those regions. At the same time, the velocity structure is less well recovered in parts of the peninsula sampled by fewer stations. This implies that the algorithm successfully adapts model resolution to the density of rays. It also adapts model resolution to the amount of noise in the travel times. Because the algorithm is computationally demanding, we modify the algorithm such that computational costs are reduced while sufficiently preserving non-linearity. We conclude that the algorithm can now be applied adequately to travel times extracted from station–station cross correlations by the Reykjanes seismic network.

A preliminary seismic travel time tomography beneath Ecuador from data of the national network

The subduction of the Nazca plate has specific features that impact the geophysical structure beneath Ecuador. In addition to the convexity of the trench between Peru and Ecuador, the Carnegie ridge and the Grijalva scarp are associated with strong heterogeneities in the plate that interfere with the subduction process. We have taken advantage of the large amount of manually picked P and S wave arrival times accumulated over time thanks to the national RENSIG network to perform an inversion for both event locations and seismic velocities over Ecuador within lat. 1.5∘N and 5∘S and long. 77∘W and 82∘W. After data filtering, the model presented in this article results from the inversion of 335,498 P and 111,457 S arrival times, corresponding to 25,410 events between 1988 and 2016. The intermediate depth seismicity outlines a continuous Wadati-Benioff zone in southern Ecuador that clearly defines the topography of the Farallon plate and its plunge towards the north-east at the Puyo nest, whereas in the northern part of the slab its pattern is more heterogeneous and is dominated by three nests at depths ranging between 75 km and 110 km, the Guayaquil, La Maná and Maldonado nests. The velocity model reveals a discontinuity of the Nazca plate along an axis oriented N110∘E, starting in the seismogenic zone at the southern limit of the occurrence of large earthquakes. This tear is associated with an overlap of the southern Farallon part of the slab by its younger northern part. We propose that it is the result of plate buckling due to lateral compression at depth caused by the sharp bend of the trench line between Peru and Ecuador.

Adjoint‐State Traveltime Tomography: Eikonal Equation‐Based Methods and Application to the Anza Area in Southern California

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

Subsurface space investigation jointly using seismic reflection and refraction tomography in urban area: a case study of Hangzhou, China

Seismic reflection method is one of the commonly used techniques in urban underground space exploration. However, in some areas, the effect of reflection seismic exploration is limited because of the difficulty of interpretation. In this paper, seismic reflection exploration and refraction travel time tomography are jointly used in Hangzhou to increase detection accuracy. The results show that different layers can be discerned for Quaternary sedimentary. The combined use of the two methods has obvious advantages and good practical value.

Characterizing a decametre-scale granitic reservoir using ground-penetrating radar and seismic methods

Abstract. Ground-penetrating radar (GPR) and seismic imaging have proven to be important tools for the characterization of rock volumes. Both methods provide information about the physical rock mass properties and geological structures away from boreholes or tunnel walls. Here, we present the results from a geophysical characterization campaign that was conducted as part of a decametre-scale hydraulic stimulation experiment in the crystalline rock volume of the Grimsel Test Site (central Switzerland). For this characterization experiment, we used tunnel-based GPR reflection imaging as well as seismic travel-time tomography to investigate the volumes between several tunnels and boreholes. The interpretation of the GPR data with respect to geological structures is based on the unmigrated and migrated images. For the tomographic analysis of the seismic first-arrival travel-time data, we inverted for an anisotropic velocity model described by the Thomsen parameters v0, ϵ and δ to account for the rock mass foliation. Subsequently, the GPR and seismic images were interpreted in combination with the geological model of the test volume and the known in situ stress states. We found that the ductile shear zones are clearly imaged by GPR and show an increase in seismic anisotropy due to a stronger foliation, while they are not visible in the p-wave (v0) velocity model. Regions of decreased seismic p-wave velocity, however, correlate with regions of high fracture density. For geophysical characterization of potential deep geothermal reservoirs, our results imply that wireline-compatible borehole GPR should be considered for shear zone characterization, and that seismic anisotropy and velocity information are desirable to acquire in order to gain information about ductile shear zones and fracture density, respectively.