We have carried out a tomographic inversion for seismic velocity in the vicinity of Uturuncu volcano (Bolivia) based on a 33-station temporary seismic network deployment. We combine travel times from earthquakes in the shallow crust with those from earthquakes on the subducting Nazca plate to broadly constrain velocities throughout the crust using the LOTOS tomography algorithm. The reliability and resolution of the tomography is verified using a series of tests on real and synthetic data. The resulting three-dimensional distributions of Vp, Vs, and Vp/Vs reveal a large tooth-shaped anomaly rooted in the deep crust and stopping abruptly 6 km below the surface. This feature exhibits very high values of Vp/Vs (up to 2.0) extending to ~80 km depth. To explain the relationship of this anomaly with the surface uplift observed in interferometric synthetic aperture radar (InSAR) data, we propose two scenarios. In the first, the feature is a pathway for liquid volatiles that convert to gas, due to decompression, at ~6 km depth, causing a volume increase. This expansion drives seismicity in the overlying crust. In the second model, this anomaly is a buoyant pulse of magma within the batholith, ascending due to gravitational instability. We propose a simplified numerical simulation to demonstrate how this second model generally supports many of the observations. We conclude that both of these scenarios might be valid and complement each other for the Uturuncu case. Based on joint analysis of the tomography results and available geochemical and petrological information, we have constructed a model of the Uturuncu magma system that illustrates the main stages of phase transitions and melting.