A new two-dimensional (2-D) speckle tracking method for displacement estimation based on the gradients of the magnitude and phase of 2-D complex correlation in a search region is presented. The novelty of this approach is that it couples the phase and magnitude gradients near the correlation peak to determine its coordinates with sub-sample accuracy in both axial and lateral directions. This is achieved with a minimum level of lateral interpolation determined from the angles between the magnitude and phase gradient vectors on the sampled (laterally interpolated) 2-D cross-correlation grid. The key result behind this algorithm is that the magnitude gradient vectors' final approach to the true peak is orthogonal to the zero-phase contour. This leads to a 2-D robust projection on the zero-phase contour that results in subsample accuracy at interpolation levels well below those needed using previously proposed methods. A full description of the 2-D, phase-coupled approach is given, including two implementations based on a geometric projection and constrained optimization. In addition, a robust fast search algorithm that allows the localization of the true peak without the need for exhaustive search is given. Experimental validation on three data sets from speckle-generating phantoms undergoing uniform diagonal motion, uniform axial deformation, and nonuniform lateral flow is given. It is shown that estimated 2-D displacement fields obtained using the phase-coupled technique display a full range of values covering the dynamic range without evidence of quantization. In comparison, a previously published method using 1-D phase-projection after lateral interpolation produces severely quantized lateral displacement fields (at the same levels of interpolation as the 2-D, phase-coupled method).