Mantle seismic tomography has historically relied on radially symmetric ID velocity models to trace ray paths through the mantle. The resulting travel time residuals are used to invert for seismic velocity perturbations around this 1D model. However, we know the Earth deviates from such ID velocity models; for example there are global variations in crustal thickness; in the age of oceanic lithosphere and presence of subducting oceanic lithosphere. In light of this, an a priori model which incorporated the three types of surface observable heterogeneity outlined above was constructed as part of this thesis. Tracing ray paths through this more heterogeneous starting model resulted in new travel time residuals which were subsequently employed in a simultaneous tomographic inversion solving for earthquake relocation parameters and slowness perturbations. This inversion method allows us to investigate whether tomography using a priori models results in improved images of mantle velocity perturbations and systematic earthquake relocations. A graphical earthquake browser was specifically written to establish, in a consistent manner, the shape of subducting oceanic lithosphere for all the major subduction zones. The resulting population of earthquakes, which best represent the shape of Wadati-Benioff zones, were subsequently interpolated into profiles following the path of oceanic lithosphere as it subducts. The temperature field in and around each profile was generated using a new analytic solution of the heat equation for subducting lithosphere, adapted to incorporate slab shape. The upper mantle a priori model was constructed on an equal area tomographic grid by combining the thermal models of the subducting lithosphere, plate cooling models of oceanic lithosphere and variations in crustal thickness away from that prescribed in a ID velocity model. Efficient 20 ray tracing through the a priori model was achieved via the adaptation of a ID ray tracer by perturbing the reference ID model, iasp91, using the a priori velocities in the cells connecting the event to the recording station for each ray. A new travel time residual was calculated and subsequently used in the simultaneous solution for slowness perturbation and earthquake relocation. So as not to bias the earthquake relocation procedure, phases were selected so as to maximise the azimuth and epicentral distance coverage, while minimising the number of duplicated ray paths which would be redundant in the inversion. The data selection resulted in some 3,450 events emitting 785,000 teleseismic P phases (bottoming in the lower mantle). The cell based SIRT inversion procedure, used to solve the standard system of linear tomographic equations, was augmented by explicit damping and smoothing matrices so as to control both poorly resolved cells and the relative importance between earthquake relocation parameters and slowness perturbations. For comparison, the ray population was also traced through the 3SMAC upper mantle model before undertaking a similar inversion. The 5° x 5° equal area, 100 km thick, cell inversions resulted in systematic earthquake relocations with an average relocation distance of= 5 km. In the upper mantle, the inversion procedure adjusts the a priori subducting slab velocity contrast, revealing images of subducting oceanic lithosphere. In the lower mantle, there is little difference between inversions produced in this thesis and those available digitally. Some of the main features are the pronounced lineations interpreted as the Farallon slab (beneath North and South America) and the Tethys (beneath Eurasia) clearly imaged between 1200 and 1500 km depth. All inversions undertaken in this thesis image hotspots throughout the upper mantle, and in places these pronounced slow features are observed passing through the upperllower mantle transition. A section through the South Pacific superswell images slow material as a continuous body, to at least 1300 km. Synthetic recovery tests indicate these hotspot features are well resolved.