An ice sheet model was constrained to reconstruct the evolution of the Greenland Ice Sheet (GrIS) from the Last Glacial Maximum (LGM) to present to improve our understanding of its response to climate change. The study involved applying a glaciological model in series with a glacial isostatic adjustment and relative sea-level (RSL) model. The model reconstruction builds upon the work of Simpson et al. (2009) through four main extensions: (1) a larger constraint database consisting of RSL and ice extent data; model improvements to the (2) climate and (3) sea-level forcing components; (4) accounting for uncertainties in non-Greenland ice. The research was conducted primarily to address data-model misfits and to quantify inherent model uncertainties with the Earth structure and non-Greenland ice. Our new model (termed Huy3) fits the majority of observations and is characterised by a number of defining features. During the LGM, the ice sheet had an excess of 4.7 m ice-equivalent sea-level (IESL), which reached a maximum volume of 5.1 m IESL at 16.5 cal ka BP. Modelled retreat of ice from the continental shelf progressed at different rates and timings in different sectors. Southwest and Southeast Greenland began to retreat from the continental shelf by ∼16 to 14 cal ka BP, thus responding in part to the Bølling-Allerød warm event (c. 14.5 cal ka BP); subsequently ice at the southern tip of Greenland readvanced during the Younger Dryas cold event. In northern Greenland the ice retreated rapidly from the continental shelf upon the climatic recovery out of the Younger Dryas to present-day conditions. Upon entering the Holocene (11.7 cal ka BP), the ice sheet soon became land-based. During the Holocene Thermal Maximum (HTM; 9-5 cal ka BP), air temperatures across Greenland were marginally higher than those at present and the GrIS margin retreated inland of its present-day southwest position by 40–60 km at 4 cal ka BP which produced a deficit volume of 0.16 m IESL relative to present. In response to the HTM warmth, our optimal model reconstruction lost mass at a maximum centennial rate of c. 103.4 Gt/yr. Our results suggest that remaining data-model discrepancies are affiliated with missing physics and sub-grid processes of the glaciological model, uncertainties in the climate forcing, lateral Earth structure, and non-Greenland ice (particularly the North American component). Finally, applying the Huy3 Greenland reconstruction with our optimal Earth model we generate present-day uplift rates across Greenland due to past changes in the ocean and ice loads with explicit error bars due to uncertainties in the Earth structure. Present-day uplift rates due to past changes are spatially variable and range from 3.5 to −7 mm/a (including Earth model uncertainty).