Abstract A widely accepted global model of glacial isostatic adjustment (GIA) is based on ICE‐6G ice loading and 1D‐viscosity model VM5a, a viscosity profile which remains at odds with mantle convection constraints such as from geoid modeling. We explore 3D‐viscosity variations in models of GIA and show that a model with lateral viscosity variations (LVVs) scaled from seismic tomography and lower mantle viscosity compatible with mantle convection constraints (∼30 times upper mantle viscosity) significantly reduces the misfit to relative sea level (∼40%), J2˙ $dot{ {J}_{2} }$ (60%), and polar motion (75%) while incurring only a small misfit penalty (∼5%) for GNSS‐uplift compared to VM5a, which has lower mantle viscosity ∼6 times upper mantle viscosity. Regional results exhibit additional tradeoffs but indicate that comparable solutions might be found with additional model refinement. The origin of a longstanding disagreement on lower mantle viscosity may thus be resolved by accounting for LVVs in models of GIA.