We present a new three-dimensional thermomechanically coupled ice sheet model of the northern hemisphere to reconstruct the Quaternary ice sheets during the last glacial cycle. The model includes basal sliding, internally calculated surface mass balance, glacial isostasy, and a treatment for marine calving. The time dependent forcing consists of temperature and precipitation anomalies from the United Kingdom Meteorological Office (UKMO) General Circulation Model scaled to the Greenland Ice Core Project (GRIP) ice core δ18O record. Model parameters were chosen to best match geomorphological inferences on Last Glacial Maximum extent and global eustatic sea level change. For our standard run we find a maximum ice volume of 57 ¿ 106 km3 at 18.5 ka cal BP. This corresponds to a eustatic sea level lowering of 110 m after correction for hydro-isostatic displacement and anomalous ice resulting from defects in the specified boundary conditions of the Paleoclimate Model Intercomparison Project (PMIP) for which the UKMO GCM results were generated. Of this 110 m, 82 m was stored in the North American ice sheet and 25 m in the Eurasian ice sheet. We determine the qualitative and quantitative response of the model from a comprehensive sensitivity study in which 11 important parameters were varied over their respective ranges of uncertainty. Model outputs comparable to the observational record were explored in detail as a linear function along the axes of parameter space of the reference model. The method reveals the dominance of climate uncertainty when modeling the Last Glacial Maximum configuration of the northern hemisphere ice sheets, but also highlights the role of ice rheology and basal processes for ice sheet thickness, and glacial isostasy and calving for the timing of maximum ice volume.