The ice flux divergence of a glacier is an important quantity to examine because it determines the rate of temporal change of its thickness. Here, we combine high-resolution ice surface velocity observations of Nioghalvfjerdsfjorden (79north) Glacier, a major outlet glacier in north Greenland, with a dense grid of ice thickness data collected with an airborne radar sounder in 1998, to examine its ice flux divergence. We detect large variations, up to 100 m/yr, in flux divergence on grounded ice that are incompatible with what we know of the glacier surface mass balance, basal mass balance and thinning rate. We examine the hypothesis that these anomalies are due to the three-dimensional flow of ice around and atop bumps and hollows in basal topography by comparing the flux divergence of three-dimensional numerical models with its surface equivalent. We find that three-dimensional effects have only a small contribution to the observed anomalies. On the other hand, if we degrade the spatial resolution of the data to 10 km the anomalies disappear. Further analysis shows that the source of the anomalies is not the ice velocity data but the interpolation of multiple tracks of ice thickness data onto a regular grid using a scheme (here block kriging) that does not conserve mass or ice flux. This problem is not unique to 79north Glacier but is common to all conventional ice thickness surveys of glaciers and ice sheets; and fundamentally limits the application of ice thickness grids to high-resolution numerical modeling of glacier flow.