The acquisition of magnetic gradient measurements in airborne surveys has led to the improved gridding of total-field data. The approach described in this paper takes the additional information from the observed gradients one step further and directly incorporates them into 3D inversion. These data spatially constrain the recovered model of which reproduces both observed total-field and gradient datasets. In this paper, we describe the methodology to invert for a common 3D susceptibility distribution. The method is demonstrated through synthetic and field examples. The end result is more accurate recovered susceptibilities and compact anomalies compared with recovered models from total-field data alone.