The long-wavelength geoid is sensitive to Earth’s mantle density structure as well as radial variations in mantle viscosity. We present a suite of inversions for the radial viscosity profile using whole-mantle models that jointly constrain the variations in density, shear- and compressional-wavespeeds using full-spectrum tomography. We use a Bayesian approach to identify a collection of viscosity profiles compatible with the geoid, while enabling uncertainties to be quantified. Depending on tomographic model parameterization and data weighting, it is possible to obtain models with either positive- or negative-buoyancy in the large low shear velocity provinces (LLSVPs). We demonstrate that whole-mantle density models in which density and $V_S$ variations are correlated imply an increase in viscosity below the transition zone, often near ~1000~km. Many solutions also contain a low-viscosity channel below 650~km. Alternatively, models in which density is less-correlated with $V_S$ – which better fit normal mode data – require a reduced viscosity region in the lower mantle. This feature appears in solutions because it reduces the sensitivity of the geoid to buoyancy variations in the lowermost mantle. The variability among the viscosity profiles obtained using different density models is indicative of the strong non-linearities in modeling the geoid and the limited resolving power of the geoid kernels. We demonstrate that linearized analyses of model resolution do not adequately capture the posterior uncertainty on viscosity. Joint and iterative inversions of viscosity, wavespeeds, and density using seismic and geodynamic observations are required to reduce bias from prior assumptions on viscosity variation and scalings between material properties.