The derivation and implementation of a three-dimensional model used to investigate the Last Glacial Maximum ice sheet across Iceland by Hubbard et al. [2006. A modelling insight into the Icelandic Last Glacial Maximum ice sheet. Quaternary Science Reviews] is described. It is applied at 2 km resolution and requires boundary distributions of topography, geothermal heat flux, surface air temperature and mass balance calculated using a temperature-index approach based on reference distributions of annual temperature and precipitation. The model enables the variables of ice thickness, stress, strain and temperature to evolve freely through time and caters for the coupling of thermally triggered basal sliding with non-local dynamics through the computation of longitudinal stresses. It is driven through perturbations in sea-level and annual precipitation and temperature. A series of contemporary experiments are initiated to validate the model against the present ice cover across Iceland. Forcing the model from ice-free conditions with the 1961–1990 (reference) climatology yields a good simulation of all the ice masses except for Vatnajökull, where the model falls well short of its present margins. However, an experiment forced from ice-free conditions for 1000 years with a 2 °C cooling perturbation, followed by 100 years of reference climatology yields a good simulation of Vatnajökull (in addition to other ice masses), implying that it is a remnant icecap, inherited from the Little Ice Age and perpetuated through strong ice elevation/mass balance coupling. An ensemble of experiments are initiated to investigate the sensitivity of the optimum LGM model isolated in Hubbard et al. (2006). Ice sheet volume and aspect ratio (but not area) are found to be sensitive to basal boundary conditions, in particular to the choice of sliding parameter and the applied geothermal conditions. Due to strong topographic control, in particular the configuration of offshore bathymetry and shelf break, ice sheet volume and area is sensitive to the calving parameter and sea-level change. However, an asymmetrical response indicates that the ice sheet is effectively decoupled from further climatic deterioration once it advances to the continental shelf break. These experiments imply that there is little latitude in the selection of model parameters which yields an ice sheet compatible with the available evidence and that the optimum LGM experiment represents a sound result. By inference, at least 63% of the optimum LGM ice sheet was grounded below sea-level implying potential instability with the onset of deglaciation.