Estimation of ice sheet mass balance from satellite altimetry requires interpolation of point-scale elevation change (dH/dt) data over the area of interest. The largest dH/dt values occur over narrow, fast-flowing outlet glaciers, where data coverage of current satellite altimetry is poorest. In those areas, straightforward interpolation of data is unlikely to reflect the true patterns of dH/dt. Here, four interpolation methods are compared and evaluated over Jakobshavn Isbrae, an outlet glacier for which widespread airborne validation data are available from NASA's Airborne Topographic Mapper (ATM). The four methods are ordinary kriging (OK), kriging with external drift (KED), where the spatial pattern of surface velocity is used as a proxy for that of dH/dt, and their spatiotemporal equivalents (ST-OK and ST-KED). KED assumes a linear relationship between spatial gradients of velocity and dH/dt, which is confirmed for both negative (Pearson's correlation rho < -0.85) and, to a lesser degree, positive (rho = 0.73) dH/dt values. When compared to ATM data, KED and ST-KED yield more realistic spatial patterns and higher thinning rates (over 20 m yr(-1) as opposed to 7 m yr(-1) for OK). Spatiotemporal kriging smooths inter-annual variability and improves interpolation in periods with sparse data coverage and we conclude, therefore, that ST-KED produces the best results. Using this method increases volume loss estimates from Jakobshavn Isbrae by up to 20% compared to those obtained by OK. The proposed interpolation method will improve ice sheet mass balance reconstructions from existing and past satellite altimeter data sets, with generally poor sampling of outlet glaciers.