We use topographically corrected interferograms, repeated global positioning system observations of locations of stakes and time series of elevation data to produce time series of high-resolution three-dimensional (3-D) ice surface motion maps for the infilling of the ice depression created by the 1996 subglacial eruption at the Gjálp volcano in Vatnajökull, Iceland. The ice inflow generated uplift in the central parts of the depression. During the first months, the uplift was much reduced by basal melting as the subglacial volcano cooled. For those motions surface-parallel ice flow cannot be assumed. The 3-D motion maps are created by an optimization process that combines the complementary datasets. The optimization is based on a Markov random-field regularization and a simulated annealing algorithm. The 3-D motion maps show the pattern of gradually diminishing ice flow into the depression. They provide a consistent picture of the 3-D motion field, both spatially and with time, which cannot be seen by separate interpretation of the complementary observations. The 3-D motion maps were used to calculate the cooling rate of the subglacial volcano for the first year after the eruption. First an uplift rate resulting solely from the inflow of ice was calculated from inferred horizontal motions. Basal melting was then estimated as the difference between the calculated uplift generated by the inflow of ice, and the observed uplift that was the combined result of ice inflow and basal melting. The basal melting was found to decline from 55 m3 s–1 (due to power of 18 GW) in January 1997 to 5 m3 s–1 (2GW) in October 1997.