  The second order partial difference equation of two variables $ \CD u:=
A_{1,1}(x) \Delta_1 \nabla_1 u + A_{1,2}(x) \Delta_1 \nabla_2 u + A_{2,1}(x)
\Delta_2 \nabla_1 u + A_{2,2}(x) \Delta_2 \nabla_2 u
  & \qquad \qquad \qquad \qquad + B_1(x) \Delta_1 u + B_2(x) \Delta_2 u =
\lambda u, $ is studied to determine when it has orthogonal polynomials as
solutions. We derive conditions on $\CD$ so that a weight function $W$ exists
for which $W \CD u$ is self-adjoint and the difference equation has polynomial
solutions which are orthogonal with respect to $W$. The solutions are
essentially the classical discrete orthogonal polynomials of two variables.
