Uturuncu volcano in southern Bolivia has been steadily inflating since at least 1965, and combined with shallow seismicity and near-summit active fumaroles represents a volcanic system showing significant signs of unrest. This study focuses on the mechanism driving the 70 km wide region of ground uplift. Using Finite Element Analysis we test for first-order parameters that constrain a viable model for the observed maximum InSAR line-of-sight uplift rate of 1 – 2 cm/yr between 1992 and 2006. Stresses and strains from pressurised finite sources are solved numerically using COMSOL Multiphysics, accounting for homogeneous and heterogeneous subsurface structure in elastic and viscoelastic rheologies. Crustal heterogeneity is constrained from seismic velocity data, which indicates a pervasive large low-velocity zone approximately 17 km below the surface. This is deduced to represent one of the world's largest known regions of partial-melt: the Altiplano-Puna Magma Body (APMB). Contrasting crustal heterogeneity and homogeneity highlights the significant effect of a mechanically weak source-depth layer. This alters surface deformation patterns by absorbing relatively more of the subsurface strain than its surrounding layers, thereby acting as a mechanical buffer. Monotonic time-dependent deformation and an anomalously high crustal heat-flux preclude elastic conditions so we induce a viscoelastic crustal rheology using the standard linear solid model. The elastic models can also only account for the spatial component of the observed uplift so their results are used solely to guide the parameters tested in the viscoelastic models. We explore a range of possible source geometries but reject spherical and oblate shapes on the grounds of their depth below the APMB and likely unsustainable pressurisation given the expected crustal mechanics. Our final preferred model suggests that pressurisation of a magma source extending upward from the APMB is causing the observed surface uplift and requires a continued increase in this pressure to explain both the observed spatial and temporal patterns. We thus also demonstrate how a pressure-time function plays a first-order role in explaining the observed temporal deformation pattern.