As temperature increases with depth and the creep resistance of rock decreases exponentially, a high-viscosity sub-lithospheric layer, just beneath the 'elastic' lithosphere is expected to exist. Depending on the temperature profile, a low-viscosity asthenosphere may also exist if the temperature deeper down gets high enough. Since the temperature profile is expected to change laterally – especially from below the oceans to cratonic areas underneath continents, rock properties of the lithosphere, high-viscosity sub-lithosphere and low-viscosity asthenosphere are expected to change laterally. Our aim is to constrain sub-lithospheric properties (depth, thickness and viscosity), lateral lithospheric thickness variations and asthenospheric properties using observed GIA data. A Coupled Laplace-Finite Element Method is used to compute gravitationally self-consistent sea level with time-dependent coastline and rotational feedback in addition to changes in deformation, gravity and the state of stress. We start with the VM5a-ICE-6G_C model combination and then modify the lithospheric, sub-lithospheric and asthenospheric properties (including lateral thickness variation) while keeping the mantle viscosities the same as VM5a. Through this study, we confirm that the sub-lithospheric and asthenospheric properties can significantly affect the predicted global relative sea level (RSL), present-day gravity rate-of-change (g-dot) and uplift rate (u-dot) in Laurentia and Fennoscandia. In addition, incorporating the elastic lithosphere with lateral thickness variation, sub-lithosphere and asthenosphere can improve the fit to global RSL, but the predicted peak values of g-dot and u-dot in Laurentia may decrease slightly but not significant enough to affect the fit to the observed data. Our results prefer an elastic lithosphere that has maximum thickness of 140 km under continental cratons but reduces to 60 km underneath the oceans. The results preferred depth of the asthenospheric bottom is around 190–200 km with asthenospheric viscosity around 1020Pa s. Finally, we show that the best laterally heterogeneous mantle model we found in previous publication when combined with the lithosphere with lateral thickness variaion gives the best fit to global RSL and peak g-dot and u-dot in Laurentia simultaneously.