Reflection Traveltime Tomography Research Articles

Eikonal-equation-based characteristic reflection traveltime tomography

Automatically picking reflection traveltimes from the prestack shot gathers and positioning and updating the reflector depths during inversion are challenging in conventional eikonal-equation-based adjoint-state reflection traveltime tomography methods. To solve these problems, we develop an eikonal-equation-based adjoint-state characteristic reflection traveltime tomography (ASCRT) method. This method automatically extracts the reflection traveltimes by sequentially applying characteristic reflector picking in the depth migration section, followed by eikonal-equation-based traveltime calculation and an event-tracking algorithm. We also develop an efficient eikonal-equation-based kinematic imaging method to rapidly position and update the reflector depths. With our ASCRT method, the characteristic reflector positions and the velocity above the reflector are alternately updated. Our ASCRT method is performed with a layer-stripping strategy that sequentially uses selected characteristic reflectors to update the velocity model from shallow to deep. Compared with the wave equation-based reflection traveltime inversion methods, the efficiency of the ASCRT method is improved by several orders of magnitude, and the memory consumption is remarkably reduced. Synthetic and field data examples are presented to demonstrate the effectiveness and efficiency of our method.

Second-order optimization for multiparameter reflection waveform inversion in acoustic VTI media

SUMMARY The seismic wave velocity in subsurface formations typically varies both spatially and directionally, resulting in a heterogeneous and anisotropic medium. To successfully image complex structures and potential reservoirs using multichannel seismic data, high-resolution interval velocity models are required. However, most controlled-source seismic recordings lack wide-angle signals, such as diving and refracted waves that propagate through deep formations. Consequently, geophysicists rely on the traveltimes and waveforms of reflections or back-scatterings in finite-frequency and finite-offset seismograms to estimate velocity and account for possible anisotropy. Recently, wave-equation based reflection waveform inversion (RWI) has emerged as an active research topic due to its ability to recover intermediate-wavelength model components that cannot be retrieved in reflection traveltime tomography. Nevertheless, the widely used gradient-type local optimization suffers from multiparameter trade-off and slow convergence of RWI in anisotropic media. Analysis of the Fréchet derivative reveals that the sensitivity of reflections to changes in long-to-intermediate wavelengths of the anisotropic velocities is jointly controlled by the specular reflection on the interface and the radiation pattern of scattered fields concerning the model parameters along the wave paths. Anatomy of the approximate Hessian reveals characteristics of parameter coupling and spatial resolution in the context of multiparameter RWI, inspiring the development of a matrix-free Gauss–Newton algorithm based on the second-order adjoint-state method for large-scale applications. Synthetic and real data examples demonstrate that the proposed approach can appropriately handle overburden heterogeneities and anisotropies, and thus improve imaging of underlying complex structures such as dipping and small-scale faults in the deep parts.

Wide‐spectrum reconstruction of a velocity model based on the wave equation reflection inversion method and its application

AbstractMost seismic data used for conventional exploration lack far‐offset signals and effective low‐frequency components. This makes it difficult for conventional full waveform inversion methods to recover the low‐wavenumber components of mid‐ and deep‐layer models. To resolve the bottleneck of ray‐theoretical reflection traveltime tomography, wave equations, such as reflection inversion methods, are receiving increasing attention. We reconstructed a wide‐spectrum velocity model based on the multiscale wave‐equation reflection inversion method for model decomposition. First, using the dynamic image warping method to obtain reflection traveltime residuals, the wave‐equation reflection traveltime inversion was used to recover the low‐wavenumber components of the background model, which also solved the cycle‐skipping problem. Second, the medium‐wavenumber component of the model is then supplemented by the reflection waveform inversion, while the reflectivity model is obtained relying on the least‐squares reverse time migration method. Based on this method, the background and perturbation models are updated alternately and iteratively. At the same time, the stratum structural tensor information was extracted using the perturbed model image to construct a preconditioned operator for the stratum structural constraint, suppress unreasonably high‐wavenumber components in the generalized gradient and improve the geological consistency of the inversion results. Testing of the model using the Sigsbee2b model and seismic field data from the East China Sea showed that, compared with the conventional reflection traveltime tomography method based on prestack depth migration, the wave‐equation reflection inversion strategy with well‐matched information from traveltimes to waveforms significantly improved the accuracy of the middle and deep velocity modelling and enhanced the imaging quality of the whole body.

Neogene stratigraphic architecture and three-dimensional velocity structure beneath the Gulf of İzmir (western Anatolia) from reflection traveltime tomography

Neogene stratigraphic architecture and three-dimensional velocity structure beneath the Gulf of İzmir (western Anatolia) from reflection traveltime tomography

Seismic reflection waveform inversion based on Gauss–Newton optimization

Abstract Owing to the lack of long-offset signals and effective low-frequency components, it is extremely difficult for classical full-waveform inversion (FWI) to obtain the long-wavelength components of medium and deep elastic parameter models. Considering ray theory, reflection traveltime tomography is restricted by resolution, and the reflection waveform inversion (RWI) method of the wave equation, and its practicability, has received greater attention. However, previous studies into RWI have primarily used first-order optimization methods, in which convergence and resolution still need substantial improvement. Based on the second-order optimization inversion theory, we deduced the reflection wave-sensitive kernel, objective function gradient and Hessian operator of background and perturbation model. Additionally, we revealed the de-fuzzification effect of the Hessian matrix on the gradient of inversion and the working mechanism of background velocity updating on direction optimization and inversion resolution improvement. Experiments on the SEG/EAGE overthrust model showed that an RWI with the Gauss–Newton method using an approximate Hessian matrix significantly improved convergence and wide-spectrum modeling capability compared with the conjugate gradient method. On the streamer seismic data of the East China Sea, the second-order optimization RWI method surpassed the commonly used reflection traveltime tomography based on pre-stack depth migration. In addition, the RTM of inversion model can depict the complex fault system in the Changjiang sag with a high resolution, which improves the imaging quality of the deep basement.

Crustal and uppermost mantle structure of the Porcupine Basin west of Ireland from seismic and gravity methods

The Porcupine basin is a failed rift that formed during the opening of the North Atlantic Ocean. Here, we provide additional insights into the crustal structure of the northern Porcupine Basin and assess the degree of crustal extension from long streamer seismic data and free-air gravity modelling. Recent re-processing of the 2013/2014 streamer data involved de-multiple and de-noise processing with pre-stack-depth migration using velocity models derived from reflection travel time tomography and full waveform inversion. Forward modelling of satellite gravity data was used to constrain the velocity and mass density of the basement and lower crust, where the seismic constraints are limited. Prior knowledge, from previous geophysical and geological studies, was included to guide the modelling and reconcile the gravity data with the seismic sections. It also served to minimize the non-uniqueness inherent in the forward problem. Crustal thinning occurs across the basin axis and increases southwards. We infer maximum stretching factor (β) at the basin axis, increasing from ∼4 to 6 in the North Porcupine at 52.5°N to >10 at ∼51°N in the Central Porcupine. The gravity modelling supports the presence of mantle serpentinisation, and the high β implies that crustal break up and mantle exhumation may have locally occurred in the Central Porcupine Basin, in agreement with previous seismic tomographic studies. Basaltic melt generation in the centre of the basin where the maximum bulk stretching factor exceeds10 in the hyper-extended crust may have led to localized sea-floor-spreading.

Robust full-waveform inversion with graph-space optimal transport: Application to 3D ocean-bottom cable Valhall data

Improving full-waveform inversion robustness to cycle-skipping has been the subject of a large number of studies. From the several families of approaches developed, one of the most documented consists in modifying the least-squares distance misfit function. From all the propositions made to improve and replace the least-squares distance, only a few of them have been applied to field data. One of the methods proposed recently, the graph space optimal transport distance, presents appealing properties for field data applications. We compare it with the least-squares distance in an analysis performed on the three-dimensional ocean bottom cable data from the Valhall field. We first perform this comparison starting the inversion from the reflection traveltime tomography model used in previous studies. We then perform a second comparison from a crude, linearly varying in-depth one-dimensional velocity model. Starting from this model, least-squares-based full-waveform inversion fails to provide a meaningful estimate of the pressure-wave velocity model due to cycle skipping. We illustrate how the graph-space optimal transport-based full-waveform inversion mitigates this issue. A meaningful estimate of the pressure-wave velocity model is obtained in the zone sampled by both diving and reflected waves, down to almost two kilometers depth.

Data-Domain Traveltime Inversion of Reflected Waves Using Segment Dynamic Image Warping

Reflection traveltime information is generally used for producing a kinematically accurate velocity model for seismic imaging. Accurately estimating the traveltime difference between the demigrated and observed data is crucial for data-domain reflection traveltime tomography. To enhance the robustness of signal's alignment, we propose a new method of segment dynamic image warping (SDIW). It first exploits windowed polynomial fitting for signal processing and then aligns the signals based on point-wise segment-to-segment matching. Compared with the conventional point-to-point matching strategy, our method is more robust and insensitive toward strong random noise. By minimizing the traveltime difference computed by SDIW, data-domain reflection traveltime inversion (DRTI) is then used to find an accurate background velocity model. Synthetic tests show that SDIW provides reliable time shifts between the demigrated and observed data even if the signal-to-noise ratio of input data is rather low. This enables DRTI to automatically reconstruct the deep background velocity model. The subsequent full waveform inversion (FWI) gradually added short-wavelength velocity structures and finally recovers a high-resolution and high-fidelity model. The results demonstrate that the DRTI provides a kinematically accurate background velocity model sufficient for following FWI to fully retrieve the detail of the subsurface structures even when the low frequencies are missing.

Application of Snell's law in reflection raytracing using the multistage fast marching method

Accurate calculations of travel times and raypaths of reflection waves are important for reflection travel time tomography. The multistage shortest path method (MSPM) and multistage fast marching method (MFMM) have been widely used in reflection wave raytracing, and both of them are characterized by high efficiency and accuracy. However, the MSPM does not strictly follow Snell's law at the interface because it treats the interface point as a sub-source, resulting in a decrease in accuracy. The MFMM achieves high accuracy by solving the Eikonal equation in local triangular mesh. However, the implementation process is complex. Here we propose a new method which uses linear interpolation to compute the incident travel time of interface points and then using Snell's law to compute the reflection travel time of grid points just above the interface. Our new method is much simpler than the MFMM; furthermore, numerical simulations show that the accuracy of the MFMM and our new method are basically the same, thus the reflection tomography algorithms which use our new method are easier to implement without decreasing accuracy. Besides, our new method can be extended easily to other grid-based raytracing methods.

Determination of a stress-dependent rock-physics model using anisotropic time-lapse tomographic inversion

In the petroleum industry, time-lapse (4D) studies are commonly used for reservoir monitoring, but they are also useful to perform risk assessment for potential overburden deformations (e.g., well shearing, cap-rock integrity). Although complex anisotropic velocity changes are predicted in the overburden by geomechanical studies, conventional time-lapse inversion workflows only deal with vertical velocity changes. To retrieve the geomechanically induced anisotropy, we have adopted a reflection traveltime tomography method coupled with a time-shift estimation algorithm of prestack data of the baseline and monitor simultaneously. For the 2D approach, we parameterize the anisotropy using five coefficients, enough to cover any type of anisotropy. Before applying the workflow to a real data set, we first study a synthetic data set based on the real data set and include velocity variations between baseline and monitor found in the literature (vertical P-wave velocity decrease in the cap rock and isotropic P-wave velocity change in the reservoir). On the synthetics, we measure the angular ray coverage necessary to retrieve the target anisotropy and observe that the retrieved anisotropies depend on the offset range. Based on a synthetic experiment, we believe that the acquisition of the real case study is suitable for performing tomographic inversion. The anisotropic velocity changes obtained on three inlines separated by 375 m are consistent and show a strong positive anomaly in the cap rock along the 45° direction (the [Formula: see text] parameter in Thomsen notation), whereas the vertical velocity change is surprisingly almost negligible. We adopt a rock-physics explanation compatible with these observations and geologic considerations: a reactivation of water-filled subvertical cracks.

Ray‐based reflection traveltime tomography using approximate stationary points

ABSTRACTReflection traveltime tomography has been used to describe subsurface velocity structures which, in practice, can be used as a background or initial model for pre‐stack depth migration or full waveform inversion. Conventional reflection traveltime tomography is performed by solving an optimization problem based on a ray‐tracing method. As a result, reflection traveltime tomography requires heavy computational efforts to carry out ray tracing and solve a large matrix equation. In addition, like most data‐domain tomography methods, reflection traveltime tomography depends on initial guesses and suffers from non‐uniqueness and uncertainty of solutions. In this research, we propose a deterministic ray‐based reflection traveltime tomography method by applying seismic interferometry. This method does not suffer from the non‐uniqueness problem and does not require a priori information on subsurface media. By adding a virtual layer (whose properties are known) on top of the real surface and applying convolution‐type interferometry, we approximately determine the stationary points (i.e., incident raypaths in the virtual layer). Then, we generate reflection points for a range of assumed velocities and estimate the velocity by considering the number of reflection points and the traveltime difference between the observed and calculated data. The reflection surface can then be recovered by using the estimated velocity. Once the first target layer is resolved, we can recover the whole media by recursively applying the same method to the lower layers. Numerical examples using surface seismic profile data for homogeneous‐layer (with a low‐velocity layer) and inhomogeneous‐layer models and real field data experiments on the Congo data set demonstrate that our method can successfully recover the velocities and depths of subsurface media without initial guesses. However, our method has some limitations for multi‐layer models because the method does not have sufficient reflection points for the deeper layers.

Parsimonious slope tomography based on eikonal solvers and the adjoint-state method

Velocity macromodel building is an essential step of the seismic imaging workflow. Indeed, obtaining acceptable results through migration or full waveform inversion is highly dependent on the kinematic accuracy of the background/initial velocity model. Two decades ago, stereotomography was proposed as an alternative to reflection traveltime tomography, the first relying on semi-automatic picking of locally coherent events associated with small reflection or diffraction segments tied to scatterers in depth by a pair of rays, while the latter on interpretive picking of laterally continuous reflections. The flexibility of stereotomography paved the way for many developments that have shown the efficiency of the method whilst emphasizing on the complementary information carried out by traveltimes and slopes of locally coherent events. A recent formulation recast stereotomography under a matrix-free formulation based on eikonal solvers and the adjoint-state method. In the latter, like in the previous works, the scatterer positions and the velocity field are updated jointly to tackle the ill-famed velocity–position coupling in reflection tomography. Following on from this adjoint-state formulation, we propose a new parsimonious formulation of slope tomography that offers the chance to restrain the problem to minimizing the residuals of a single data class being a slope, in search of a sole parameter class being the subsurface velocity field. This parsimonious formulation results from a variable projection, which is implemented by enforcing a consistency between the scatterer coordinates and the velocity macromodel through migration of kinematic attributes. We explain why the resulting reduced-parametrization inversion is more suitable for tomographic problems than the most common joint inversion strategy. We benchmark our method against the complex Marmousi model along with a validation through time domain full waveform inversion and then present the results of a field data case study.

Fast 3D frequency-domain full-waveform inversion with a parallel block low-rank multifrontal direct solver: Application to OBC data from the North Sea

Wide-azimuth long-offset ocean bottom cable (OBC)/ocean bottom node surveys provide a suitable framework to perform computationally efficient frequency-domain full-waveform inversion (FWI) with a few discrete frequencies. Frequency-domain seismic modeling is performed efficiently with moderate computational resources for a large number of sources with a sparse multifrontal direct solver (Gauss-elimination techniques for sparse matrices). Approximate solutions of the time-harmonic wave equation are computed using a block low-rank (BLR) approximation, leading to a significant reduction in the operation count and in the volume of communication during the lower upper (LU) factorization as well as offering great potential for reduction in the memory demand. Moreover, the sparsity of the seismic source vectors is exploited to speed up the forward elimination step during the computation of the monochromatic wavefields. The relevance and the computational efficiency of the frequency-domain FWI performed in the viscoacoustic vertical transverse isotropic (VTI) approximation was tested with a real 3D OBC case study from the North Sea. The FWI subsurface models indicate a dramatic resolution improvement relative to the initial model built by reflection traveltime tomography. The amplitude errors introduced in the modeled wavefields by the BLR approximation for different low-rank thresholds have a negligible footprint in the FWI results. With respect to a standard multifrontal sparse direct factorization, and without compromise of the accuracy of the imaging, the BLR approximation can bring a reduction of the LU factor size by a factor of up to three. This reduction is not yet exploited to reduce the effective memory usage (ongoing work). The flop reduction can be larger than a factor of 10 and can bring a factor of time reduction of around three. Moreover, this reduction factor tends to increase with frequency, namely with the matrix size. Frequency-domain viscoacoustic VTI FWI can be viewed as an efficient tool to build an initial model for elastic FWI of 4C OBC data.

Imaging through shallow gas: Integrating broadband acquisition, processing and high-end model building for improved imaging of deeper targets

SUMMARY We present a case study offshore Malaysia, shallow gas features in the overburden distort the seismic imaging at the target level. While a multifaceted approach involving a combination of seismic acquisition and processing strategies were used to improve the bandwidth of the seismic data, particularly for the low-frequency content of the seismic image, several distortions still existed at the target level. The prominent structural sag evident at the reservoir level is a typical indication that the overlying shallow gas velocity model needed to be resolved and incorporated into a depth migration algorithm. To resolve the transversally and laterally variant velocity features in the shallow gas areas, a solution that consisted of full waveform inversion (FWI) and high-resolution reflection traveltime tomography was utilized to produce an accurate compressional velocity model. To further resolve the amplitude and phase distortions at the reservoir level due to shallow gas effects, Q tomography was incorporated into the model building phase to derive a space-variant 1/Qmodel and Q compensation was integrated within depth migration. The integrated approach of broadband receiver acquisition, data processing strategies and high-end Earth model building has cumulatively improved the imaging of the reservoir below the shallow gas anomalies.

TOMO3D: 3-D joint refraction and reflection traveltime tomography parallel code for active-source seismic data—synthetic test

We present a new 3-D traveltime tomography code (TOMO3D) for the modelling of active-source seismic data that uses the arrival times of both refracted and reflected seismic phases to derive the velocity distribution and the geometry of reflecting boundaries in the subsurface. This code is based on its popular 2-D version TOMO2D from which it inherited the methods to solve the forward and inverse problems. The traveltime calculations are done using a hybrid ray-tracing technique combining the graph and bending methods. The LSQR algorithm is used to perform the iterative regularized inversion to improve the initial velocity and depth models. In order to cope with an increased computational demand due to the incorporation of the third dimension, the forward problem solver, which takes most of the run time (∼90 per cent in the test presented here), has been parallelized with a combination of multi-processing and message passing interface standards. This parallelization distributes the ray-tracing and traveltime calculations among available computational resources. The code's performance is illustrated with a realistic synthetic example, including a checkerboard anomaly and two reflectors, which simulates the geometry of a subduction zone. The code is designed to invert for a single reflector at a time. A data-driven layer-stripping strategy is proposed for cases involving multiple reflectors, and it is tested for the successive inversion of the two reflectors. Layers are bound by consecutive reflectors, and an initial velocity model for each inversion step incorporates the results from previous steps. This strategy poses simpler inversion problems at each step, allowing the recovery of strong velocity discontinuities that would otherwise be smoothened.

Data-domain wave equation reflection traveltime tomography

Estimation of an accurate macro velocity model plays an important role in seismic imaging and model parameter inversion. Full waveform inversion (FWI) is the classical data-domain inversion method. However, the misfit function of FWI is highly nonlinear, and the local optimization cannot prevent convergence of the misfit function toward local minima. To converge to the global minimum, FWI needs a good initial model or reliable low frequency component and long offset data. In this article, we present a wave-equation-based reflection traveltime tomography (WERTT) method, which can provide a good background model (initial model) for FWI and (least-square) pre-stack depth migration (LS-PSDM). First, the velocity model is decomposed into a low-wavenumber component (background velocity) and a high-wavenumber component (reflectivity). Second, the primary reflection wave is predicted by wave-equation demigration, and the reflection traveltime is calculated by an automatic picking method. Finally, the misfit function of the l2-norm of the reflection traveltime residuals is minimized by a gradient-based method. Numerical tests show that the proposed method can invert a good background model, which can be used as an initial model for LS-PSDM or FWI.

Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation

Computationally efficient 3-D frequency-domain full waveform inversion (FWI) is applied to ocean-bottom cable data from the Valhall oil field in the visco-acoustic vertical transverse isotropic (VTI) approximation. Frequency-domain seismic modelling is performed with a parallel sparse direct solver on a limited number of computer nodes. A multiscale imaging is performed by successive inversions of single frequencies in the 3.5–10 Hz frequency band. The vertical wave speed is updated during FWI while density, quality factor QP and anisotropic Thomsen's parameters δ and ϵ are kept fixed to their initial values. The final FWI model shows the resolution improvement that was achieved compared to the initial model that was built by reflection traveltime tomography. This FWI model shows a glacial channel system at 175 m depth, the footprint of drifting icebergs on the palaeo-seafloor at 500 m depth, a detailed view of a gas cloud at 1 km depth and the base cretaceous reflector at 3.5 km depth. The relevance of the FWI model is assessed by frequency-domain and time-domain seismic modelling and source wavelet estimation. The agreement between the modelled and recorded data in the frequency domain is excellent up to 10 Hz although amplitudes of modelled wavefields propagating across the gas cloud are overestimated. This might highlight the footprint of attenuation, whose absorption effects are underestimated by the homogeneous background QP model (QP = 200). The match between recorded and modelled time-domain seismograms suggests that the inversion was not significantly hampered by cycle skipping. However, late arrivals in the synthetic seismograms, computed without attenuation and with a source wavelet estimated from short-offset early arrivals, arrive 40 ms earlier than the recorded seismograms. This might result from dispersion effects related to attenuation. The repeatability of the source wavelets inferred from data that are weighted by a linear gain with offset is dramatically improved when they are estimated in the FWI model rather than in the smooth initial model. The two source wavelets, estimated in the FWI model from data with and without offset gain, show a 40 ms time-shift, which is consistent with the previous analysis of the time-domain seismograms. The computational efficiency of our frequency-domain approach is assessed against a recent time-domain FWI case study performed in a similar geological environment. This analysis highlights the efficiency of the frequency-domain approach to process a large number of sources and receivers with limited computational resources, thanks to the efficiency of the substitution step performed by the direct solver. This efficiency can be further improved by using a block-low rank version of the multifrontal solver and by exploiting the sparsity of the source vectors during the substitution step. Future work will aim to update attenuation and density at the same time of the vertical wave speed.

Shot- and angle-domain wave-equation traveltime inversion of reflection data: Synthetic and field data examples

Full-waveform inversion requires the accurate simulation of the dynamics and kinematics of wave propagation. This is difficult in practice because the amplitudes cannot be precisely reproduced for seismic waves in the earth. Wave-equation reflection traveltime tomography (WT) is proposed to avoid this problem by directly inverting the reflection-traveltime residuals without the use of the high-frequency approximation. We inverted synthetic traces and recorded seismic data for the velocity model by WT. Our results demonstrated that the wave-equation solution overcame the high-frequency approximation of ray-based tomography, was largely insensitive to the accurate modeling of amplitudes, and mitigated problems with ambiguous event identification. The synthetic examples illustrated the effectiveness of the WT method in providing a highly resolved estimate of the velocity model. A real data example from the Gulf of Mexico demonstrated these benefits of WT, but also found the limitations in traveltime residual estimation for complex models.

2D seismic reflection tomography in strongly anisotropic media

Seismic traveltime tomography is an effective method to reconstruct underground anisotropic parameters. Currently, most anisotropic tomographic methods were developed under the assumption of weak anisotropy. The tomographic method proposed here can be implemented for imaging subsurface targets in strongly anisotropic media with a known tilted symmetry axis, since the adopted ray tracing method is suitable for anisotropic media with arbitrary degree. There are three kinds of reflection waves (qP, qSV and qSH waves) that were separately used to invert the blocky abnormal body model. The reflection traveltime tomographiy is developed here because a surface observation system is the most economical and practical way compared with crosswell and VSP. The numerical examples show that the traveltimes of qP reflection wave have inverted parameters successfully. Traveltimes of qSV reflection wave have inverted parameters successfully, with the exception of the since it is less sensitive than other parameters. Traveltimes of qSH reflection wave also have inverted parameters successfully. In addition, we find that the velocity sensitivity functions (derivatives of phase velocity with respect to elastic moduli parameters) and raypath illuminating angles have a great influence on the qualities of tomograms according to the inversion of theoretical models. Finally, the numerical examples confirm that the reflection traveltime tomography can be applied to invert strongly anisotropic models.

Joint transmission and reflection traveltime tomography using the fast sweeping method and the adjoint-state technique

SUMMARY We present a joint transmission and reflection traveltime tomography algorithm based on the Fast Sweeping Method and the adjoint-state technique. In contrast to classical ray based tomography, this algorithm utilizes a grid-based Eikonal equation solver to circumvent the non-linearity of conventional ray shooting and bending approaches in complex media. The adjoint-state technique is used to obtain the gradient of the objective function without the explicit estimation of the Frechet derivative matrix, which is usually computationally prohibitive for large-scale problems. When combined with Huygens’ Principle, the tomographic inversion can simultaneously use direct and reflected arrivals to optimize a final velocity model, further mitigate the ambiguity of the inverse problem and reveal deeper structures not visible to transmission tomography alone. In this paper, we describe the theoretical basis of our algorithm, evaluate its performance on synthetic models, and then apply it to a 20 km long 2-D seismic survey acquired in the Mackenzie Delta, Northwest Territories of Canada. The subsurface at that location is characterized by a thick permafrost (600 m) comprising high- and low-velocity areas associated with thermokarst lakes. Our results show the potential of the joint tomography in characterizing multi-scale heterogeneous velocity structures within the permafrost.