Abstract Improving the reliability and accuracy in the determination of the depth to velocity discontinuities in the Earth's mantle is essential for a better understanding of mantle dynamics, as well as for addressing such fundamental questions as the origin of mantle plumes and fate of subducted slabs. Most existing techniques that utilize receiver function stacking, which is perhaps the most-commonly used method to image mantle discontinuities, suffer from strong trade-offs between the depth and velocity anomalies above the discontinuities. Here we propose and test a procedure that utilizes both the P-to-S converted phase (Pds) and the multiply reflected and converted phase (Ppds) at the discontinuities to simultaneously determine the depth of mantle discontinuities and velocity anomalies in the overlying layer. The procedure includes masking the strong PP arrivals prior to computing receiver functions, computing non-plane wave travel-times for Pds and Ppds for accurate moveout corrections, and utilizing the discrepancies in the apparent discontinuity depths from Pds and Ppds to simultaneously estimate velocity anomalies and discontinuity depths that are independent of velocity anomalies. Application of the procedure to data recorded by stations in seven radius = 1.5° circles along a 780 km N-S profile centered at the Yellowstone hotspot reveals lower-than-normal upper mantle and MTZ velocities beneath Yellowstone, a positive temperature anomaly of 40-190 °C in the vicinity of the 410 km discontinuity (d410), and a lower-than-normal temperature in the vicinity of the 660 km discontinuity (d660). The perceived hotspot is not associated with a localized depression of the d410 or an uplift of the d660, suggesting that the proposed mantle plume does not traverse the MTZ beneath the profile.