We propose an application of joint tomographic inversion using both transmitted and reflected arrivals to map velocity structures in thick permafrost areas of the Mackenzie Delta, Northwest Territories of Canada. Our tomography algorithm combines a grid-based solver of the eikonal equation, Huygens' principle, and the adjoint method. The grid-based solver assigns a traveltime to each grid point and avoids the shadow-zone problem in classical ray tracing and is well-adapted for parallelization. The adjoint method used in the inversion provides the gradient of a given objective function without explicitly estimating the Fréchet derivative matrix. When combined with Huygens' principle, the tomographic inversion can simultaneously use first and later arrivals to optimize a final velocity model. We first demonstrate the performance of the joint tomography algorithm on a two-dimensional synthetic model with velocity variations typical of permafrost. The method is then applied to a 2D seismic survey covering over 20 km with offsets up to 4 km, acquired in the Mackenzie Delta. The subsurface at that location is characterized by a thick permafrost (600 m) comprising high-velocity and low-velocity areas associated with thermokarst lakes. Our results show the potential of the joint tomography in characterizing multiscale heterogeneous velocity structures within the permafrost.