In this paper we extend the Shepard-Bernoulli operators introduced in [6] to the bivariate case. These new interpolation operators are realized by using local support basis functions introduced in [23] instead of classical Shepard basis functions and the bivariate three point extension [13] of the generalized Taylor polynomial introduced by F. Costabile in [11]. The new operators do not require either the use of special partitions of the node convex hull or special structured data as in [8]. We deeply study their approximation properties and provide an application to the scattered data interpolation problem; the numerical results show that this new approach is comparable with the other well known bivariate schemes QSHEP2D and CSHEP2D by Renka [34, 35].