We perform joint nonlinear inversions of glacial isostatic adjustment (GIA) data, including the following: postglacial decay times in Canada and Scandinavia, the Fennoscandian relaxation spectrum (FRS), late-Holocene differential sea level (DSL) highstands (based on recent compilations of Australian sea level histories), and the rate of change of the degree 2 zonal harmonic of the geopotential, $J2$. Resolving power analyses demonstrate the following: (1) the FRS constrains mean upper mantle viscosity to be ∼3 × 10$\mathrm{<mrow\; data-mjx-texclass="ORD">20</mrow>}$ Pa s, (2) postglacial decay time data require the average viscosity in the top ∼1500 km of the mantle to be 10$\mathrm{<mrow\; data-mjx-texclass="ORD">21</mrow>}$ Pa s, and (3) the $J2$ datum constrains mean lower mantle viscosity to be ∼5 × 10$\mathrm{<mrow\; data-mjx-texclass="ORD">21</mrow>}$ Pa s. To reconcile (2) and (3), viscosity must increase to 10$\mathrm{<mrow\; data-mjx-texclass="ORD">22</mrow>}$-10$\mathrm{<mrow\; data-mjx-texclass="ORD">23</mrow>}$ Pa s in the deep mantle. Our analysis highlights the importance of accurately correcting the $J2$ observation for modern glacier melting in order to robustly infer deep mantle viscosity. We also perform a large series of forward calculations to investigate the compatibility of the GIA data sets with a viscosity jump within the lower mantle, as suggested by geodynamic and seismic studies, and conclude that the GIA data may accommodate a sharp jump of 1-2 orders of magnitude in viscosity across a boundary placed in a depth range of 1000-1700 km but does not require such a feature. Finally, we find that no 1-D viscosity profile appears capable of simultaneously reconciling the DSL highstand data and suggest that this discord is likely due to laterally heterogeneous mantle viscosity, an issue we explore in a companion study.