We present a new 2D speckle tracking method for displacement estimation based on the gradients of the magnitude and phase of 2D complex correlation in a search region in two consecutive frames. The novelty of this approach is that it couples the phase and magnitude gradients near the correlation peak to determine the true maximum with sub-sample accuracy in both axial and lateral directions. It is shown that, for wide range of speckle SNR values, both magnitude and phase gradients of the 2D cross correlation near the true correlation peak are well behaved. Furthermore, the magnitude gradient vectors final approach to the true peak is either tangential or orthogonal to the zero-phase contour. This leads to an efficient and robust, directed gradient search algorithm for determination of the true peak resulting in a very accurate estimate of the 2D displacement vector. This algorithm was tested using computer simulations (Field II) as well as tissue mimicking phantoms undergoing deformation and in vivo images of human liver under a variety of breathing conditions. Comparisons with standard interpolation based approach and 2D cross spectrum based method in  were made for accuracy, speed, and robustness. Compared to standard interpolation, the phase-coupled approach produces superior results using much smaller interpolation factors. Compared to the 2D cross spectrum method, it is more robust, more computationally efficient, and produces finer lateral displacement estimation with lower variance.