This paper is concerned with a $C^1$-conforming Gauss collocation approximation to the solution of a model two-dimensional elliptic boundary problem. Superconvergence phenomena for the numerical solution at mesh nodes, at roots of a special Jacobi polynomial, and at the Lobatto and Gauss lines are identified with rigorous mathematical proof, when tensor products of $C^1$ piecewise polynomials of degree not more than $k, k ≥ 3$ are used. This method is shown to be superconvergent with $(2k−2)$-th order accuracy in both the function value and its gradient at mesh nodes, $(k+2)$-th order accuracy at all interior roots of a special Jacobi polynomial, $(k+1)$-th order accuracy in the gradient along the Lobatto lines, and $k$-th order accuracy in the second-order derivative along the Gauss lines. Numerical experiments are presented to indicate that all the superconvergence rates are sharp.