We propose an irrotationality-preserving total variation algorithm to solve the two dimensional (2D) phase unwrapping problem, which occurs in Interferometric Synthetic Aperture Radar (InSAR) imaging and other problems. Total variation methods aim at denoising the phase derivatives to reconstruct the absolute phase. We supplement these methods by adding an additional constraint driving the curl of the gradient of the 2D phase map to zero, i.e. imposing the irrotationality of the gradient map by suitably constructing a cost function which we then minimize. We test our method and compare with existing methods on several synthetic surfaces specific to the problem of InSAR imaging for different noise levels. We report better estimates of unwrapped phase maps for the terrains simulated and for all noise levels with a two-fold improvement in terms of root mean square (RMS) error in high noise level scenarios. © 2018 IEEE.