Uncertainty in present-day glacial isostatic adjustment (GIA) rates represents at least 44% of the total gravity-based ice mass balance signal over Antarctica. Meanwhile, physical couplings between solid Earth, sea level and ice dynamics enhance the dependency of the spatiotemporally varying GIA signal on three-dimensional variations in mantle rheology. Improved knowledge of thermomechanical mantle structure is therefore required to refine estimates of current and projected ice mass balance. Here, we present a Bayesian inverse method for self-consistently mapping shear-wave velocities from high-resolution adjoint tomography into thermomechanical structure using calibrated parameterisations of anelasticity at seismic frequency. We constrain the model using regional geophysical data sets containing information on upper mantle temperature, attenuation and viscosity structure. Our treatment allows formal quantification of parameter covariances, and naturally permits propagation of material parameter uncertainties into thermomechanical structure estimates. We find that uncertainty in steady-state viscosity structure at 150 km depth can be reduced by 4–5 orders of magnitude compared with a forward-modeling approach neglecting covariance between viscoelastic parameters. By accounting for the dependence of apparent viscosity on loading timescale, we find good agreement between our estimates of mantle viscosity beneath West Antarctica, and those derived from satellite GPS. Direct access to temperature structure allows us to estimate lateral variations in lithosphere-asthenosphere boundary (LAB) depth, geothermal heat flow (GHF), and associated uncertainties. We find evidence for shallow LAB depths (63 ± 13 km), and high GHF (76 ± 7 mW m−2) beneath West Antarctica that, combined with low asthenospheric viscosities, indicate a highly dynamic response to ice mass loss.